Keywords

• Evolutionary computations;
• Evolutionary multiple objective optimisation;
• Performance assessment;
• Hypervolume measure;
• OR in aerodynamic industries

1. Introduction

Dealing with preferences can be a difficult task in multiobjective optimisation, especially if the knowledge about the structure of the set of optimal solutions is incomplete. By adopting the a posteriori approach, the decision maker postpones the final decision on preferences between conflicting objectives until a set of compromise solutions is known. Such sets can be obtained by Pareto optimisation [1] and [2].

In Pareto optimisation, the purpose is to optimise a problem of m   objective functions f = (f  ₁, … , f  _m) with f_i:S→R,i∈{1,…,m}$f_i:\mathbb{S}\to \mathbb{R}\text{,}i\in \{1\text{,}\dots \text{,}m\}$, and S$\mathbb{S}$ denoting a search (or: decision) space. Here, only minimisation is considered. In the examined examples, the decision space is a subset of Rⁿ${\mathbb{R}}^n$, though the proposed algorithms can easily be generalised to other spaces. Since objectives can be conflicting, instead of searching for a single best solution, we seek for a set of good compromise solutions. Among solution vectors y,y^′∈R^m$\mathbf{y}\text{,}{\mathbf{y}}^{\prime }\in {\mathbb{R}}^m$ a partial order is defined by the dominance relation “≺”:

equation(1)

$\mathbf{y}\prec {\mathbf{y}}^{\prime }(\mathbf{y}\text{dominates}{\mathbf{y}}^{\prime }):\iff \forall i\in \{1\text{,}\dots \text{,}m\}:y_i⩽y_i^{\prime }\text{and}\exists j\in \{1\text{,}\dots \text{,}m\}:y_j<y_j^{\prime }\text{.}$

Turn MathJax on

For x,x^′∈S$\mathbf{x}\text{,}{\mathbf{x}}^{\prime }\in \mathbb{S}$, we write x≺x^′:⇔f(x)≺f(x^′)$\mathbf{x}\prec {\mathbf{x}}^{\prime }:\iff \mathbf{f}(\mathbf{x})\prec \mathbf{f}({\mathbf{x}}^{\prime })$. The minimal elements of the dominance relation are called Pareto-optimal   and the set of Pareto-optimal decision vectors is called Pareto  (-optimal  ) set  : {x∈S|∄x^′∈S:x^′≺x}$\{\mathbf{x}\in \mathbb{S}|\nexists {\mathbf{x}}^{\prime }\in \mathbb{S}:{\mathbf{x}}^{\prime }\prec \mathbf{x}\}$. The corresponding image under f in the solution space R^m${\mathbb{R}}^m$ is called Pareto front  . A point in the search space is called non-dominated   within a set A⊆S$A\subseteq \mathbb{S}$ if there is no point in A dominating it.

In practice, it is rarely possible to compute the true Pareto set with any algorithmic method, because of the black-box character of many objective functions and/or the complexity of the search space. Instead, algorithms are intended to yield approximations of the Pareto front and thereby the Pareto set.

Among Pareto optimisation techniques, EMOA established themselves in a leading position concerning research interests and progress. EMOA are applicable to a wide range of problems and their population based approach assures an inherent set of best qualified solutions. Therefore, they are well-suited methods, whenever analytical methods fail to obtain the precise Pareto set.

A variety of quality indicators is used in order to measure the quality of Pareto front approximations. Among them, the hypervolume measure   or S$\mathcal{S}$metric   is of outstanding importance. It is a quality indicator that rewards the convergence towards the Pareto front as well as the representative distribution of points along the front. The hypervolume measure was originally proposed by Zitzler and Thiele [3], who called it the size of dominated space   [4]. Let Λ   denote the Lebesgue measure, then the S$\mathcal{S}$ metric is defined as

equation(2)

$\mathcal{S}(B\text{,}{\mathbf{y}}_{\text{ref}})≔\Lambda \left(\bigcup _{\mathbf{y}\in B}\{{\mathbf{y}}^{\prime }|\mathbf{y}\prec {\mathbf{y}}^{\prime }\prec {\mathbf{y}}_{\text{ref}}\}\right)\text{,}B\subseteq {\mathbb{R}}^m\text{.}$

Turn MathJax on

Here, y_ref∈R^m${\mathbf{y}}_{\text{ref}}\in {\mathbb{R}}^m$ denotes a reference point that should be dominated by all Pareto-optimal solutions.

Fleischer [5] proved that, given a finite search space and a reference point, maximisation of the hypervolume measure is equivalent to finding the Pareto set. Therefore, a set of bounded size which obtains the maximal S$\mathcal{S}$ metric value that is possible for the given size, only consists of Pareto-optimal solutions. It has been empirically observed that for a fixed number of points, the maximisation of the hypervolume metric yields subsets of the Pareto front which are well-distributed (cf. [6] and [7]).

The absolute value of the S$\mathcal{S}$ metric depends on the choice of the reference point. This dependence can influence the ranking of Pareto front approximations (cf. [8]). Despite this disadvantage, the S$\mathcal{S}$ metric was used in several comparative studies of EMOA, e.g. [4], [9] and [10].

A rather nice characteristic of the hypervolume measure is, that the ranking of the solutions is invariant to linear scaling of the objective functions. Furthermore, the reference point only influences the relative contribution of the points containing worst values (boundary points). There are further theoretical qualities of the S$\mathcal{S}$ metric which make it a rather fair measure, and preferable to other performance metrics [11].

Although it seems straightforward to consider the quality indicator directly in the selection operator, with a few exceptions [12], [13] and [7], this option has hardly been considered for the design of EMOA, yet. This article presents an approach proposed recently by the authors [7], namely the selection of solution candidates according to their contribution to the dominated hypervolume. It deepens the discussion by emphasising on problems with three objectives.

The article is structured as follows. In Section 2, we provide a comprehensive description of the SMS-EMOA and relate it to similar state-of-the-art approaches. A detailed discussion on the algorithm’s behaviour and qualities is given, including implementation as well as theoretical aspects. Besides, an efficient modification of the SMS-EMOA is presented which has an improved runtime and in some cases even higher quality than the basic SMS-EMOA. In addition, efficient algorithms for computing hypervolume contributions of non-dominated solutions are discussed, with a focus on the two- and three-dimensional case. In Section 3, we compare the performance of the SMS-EMOA to conventional EMOA on the most common benchmark problems for evolutionary multiobjective optimisation. The studies of design optimisation tasks from aeronautic industry, featuring two and three objectives, confirm the applicability of the new approach to real-world problems (Section 4). Only two- and three-dimensional problems are covered, since these cases are most frequently addressed in real-world applications. Finally, Section 5 outlines the main conclusions of this work.

2. The S$\mathcal{S}$ metric selection (SMS-) EMOA

The SMS-EMOA has been designed to cover a maximal hypervolume with a finite number of points. Furthermore, we wanted to diminish the problem of choosing the right reference point and aimed at a simple algorithm that can be parallelised in an asynchronous manner and extended by knowledge integrating techniques, such as approximate function evaluations. The SMS-EMOA combines ideas borrowed from other EMOA, like the well-established NSGA-II [14] and archiving strategies presented by Knowles et al. [6] and [8]. It is a steady-state algorithm founded on two pillars: (1) non-dominated sorting is used as a ranking criterion and (2) the hypervolume is applied as selection criterion to discard that individual, which contributes the least hypervolume to the worst-ranked front.

2.1. Details of the SMS-EMOA

The basic algorithm is described in Algorithm 1. Starting with an initial population of μ   individuals, a new individual is generated by means of randomised variation operators. ¹ The new individual will become a member of the next population, if replacing another individual leads to a higher quality of the population with respect to the S$\mathcal{S}$ metric.

The procedure Reduce described in Algorithm 2 selects the μ   individuals of the subsequent population. The algorithm fast-nondominated-sort used in NSGA-II [14] is applied to partition the population in v   sets R₁,…,R_v${\mathcal{R}}_1\text{,}\dots \text{,}{\mathcal{R}}_v$ according to the non-dominated sorting defined by Goldberg [15]. The subsets are called fronts and are provided with an index representing a hierarchical order (the level of domination) whereas the solutions within each front are mutually non-dominated. The first subset contains all non-dominated solutions of the original set Q  . The second front consists of individuals that are non-dominated in the set (Q⧹R₁)$(Q⧹{\mathcal{R}}_1)$, e.g. each member of R₂${\mathcal{R}}_2$ is dominated by at least one member of R₁${\mathcal{R}}_1$. More general, the ith front consists of individuals that are non-dominated if the individuals of the fronts j with j < i were removed from Q.

Afterwards, one individual is discarded from the worst ranked front. Whenever this front comprises |R_v|>1$|{\mathcal{R}}_v|>1$ individuals, the individual s∈R_v$s\in {\mathcal{R}}_v$ is eliminated that minimises

equation(3)

Δ_S(s,R_v)≔S(R_v)-S(R_v⧹{s}).${\Delta }_{\mathcal{S}}(s\text{,}{\mathcal{R}}_v)≔\mathcal{S}({\mathcal{R}}_v)-\mathcal{S}({\mathcal{R}}_v⧹\{s\})\text{.}$

Turn MathJax on

The value of Δ_S(s,R_v)${\Delta }_{\mathcal{S}}(s\text{,}{\mathcal{R}}_v)$ can be interpreted as the exclusive contribution of s   to the S$\mathcal{S}$ metric value of its appropriate front. By definition of Δ_S(s,R_v)${\Delta }_{\mathcal{S}}(s\text{,}{\mathcal{R}}_v)$, an individual, which dominates another is always kept and no non-dominated individual is replaced by a dominated one. This measure keeps those individuals which maximise the population’s S$\mathcal{S}$ metric value, which implies that the covered hypervolume of a population cannot decrease by application of the Reduce operator. Thus, for Algorithm 1 the following invariant holds:

equation(4)

S(P_t)⩽S(P_t+1).$\mathcal{S}(P_t)⩽\mathcal{S}(P_{t+1})\text{.}$

Turn MathJax on

2.1.1. Steady-state selection

Due to the high computational effort of the hypervolume calculation (see Section 2.3), a steady state selection scheme is used. Since only one individual is created, only one has to be deleted from the population in each generation. Thus, the selection operator has to compute at most μ   + 1 values of the S$\mathcal{S}$ metric (exactly μ   + 1 values in case all solutions are non-dominated). These are the values of the subsets of the worst ranked front, in which one point of the front is left out, respectively. A (μ   + λ  ) selection scheme would require the calculation of $\left(\begin{array}{c}\hfill \mu +\lambda \hfill \\ \hfill \mu \hfill \end{array}\right)$ possible S$\mathcal{S}$ metric values to identify an optimally composed population, maximising the S$\mathcal{S}$ metric net value. Besides, a steady-state scheme can easily be parallelised by means of asynchronously distributed function evaluations.

2.1.2. Population size

In contrast to other strategies that store non-dominated individuals in an archive, the SMS-EMOA keeps a population of non-dominated and dominated individuals at constant size. Keeping only non-dominated individuals, might lead to small or even single-membered populations and thus to a crucial loss of diversity. Modifications of SMS-EMOA with dynamic population size have been considered in [16]. These variants are interesting whenever an acceleration of the evolution process in its beginning is desired. To avoid the loss of diversity, we recommend to define a lower bound for the population size.

2.1.3. Handling of boundary solutions

In a two-dimensional solution space, the reference point y_ref of the S$\mathcal{S}$ metric is only needed to calculate the hypervolume of the two extremal points that consist of one best and one worst objective value, respectively. We decided to omit y_ref and always keep these extremal solutions.

In the case of more objective functions, there can be far more points, whose hypervolume contributions depend on the choice of the reference point. These exactly are the boundary points which contain atleast one worst objective value. It is not intended to keep all boundary solutions, thus contributions are calculated for them as well. We used a dynamic handling of the reference point, where y_ref is recalculated in each generation. It is chosen as the vector of the currently worst objective values increased by 1.0. Thus, the contribution of a point that contains a worst objective value depends on the remaining objectives as the difference to the reference point behaves neutral in the product. Thus, good non-extremal solutions might outperform extremal ones concerning the hypervolume. First observations advised, that the dynamic handling of the reference point is preferable to a static one.

2.2. Selection variants of SMS-EMOA

A competitive study of several selection modi within the SMS-EMOA framework was presented in [16]. The concepts all feature the contributing hypervolume in combination with non-dominated sorting, the common dominance criterion, or the number of dominating points. The latter will be described in more detail here.

The number of dominating points d(s, P(t)) is the number of points from set P(t) that dominate point s, formally

equation(5)

d(s,P(t))≔|{y∈P(t)|y≺s}|.$d(s\text{,}P(t))≔|\{y\in P(t)|y\prec s\}|\text{.}$

Turn MathJax on

A modified Reduce procedure is given in Algorithm 3. The number of dominating points d  (s  , P  (t  )) is used as selection criterion in case of dominated solutions and the contributing hypervolume Δ_S${\Delta }_{\mathcal{S}}$ is applied whenever all individuals of P  (t  ) are non-dominated. If the population consists of more than one front, the individual with the highest d  (s  , P  (t  )) value among the solutions of the worst ranked front is discarded. Otherwise, the d  (s  , P  (t  )) values of all individuals equal zero and the Δ_S${\Delta }_{\mathcal{S}}$ selection is used instead.

One motivation for developing such a measure is the smaller runtime complexity compared to the hypervolume measure. Moreover, the motivation was to achieve a different ranking of dominated solutions in order to emphasise on sparsely filled regions of the solution space. Here, individuals are kept for proceeding generations to fill gaps in the Pareto front approximation. The contributing hypervolume is applied to distribute solutions well on the front they are lying on. Though the final purpose is to distribute solutions well on the first non-dominated front (R₁${\mathcal{R}}_1$), this is not an end in itself on the other fronts. The d(s, P(t)) measure favours solutions located in those areas, where the better fronts are sparsely populated. The idea is that emanating offspring solutions may rise in rank to better fronts and fill those vacancies. In areas where the non-dominated front is densely populated, it is of no use to keep dominated individuals.

The measure is visualised by means of an example in Fig. 1. The contributing hypervolume measure favours y₈ over y₉ but the d(s, P(t)) measure favours y₉ since there is only one dominating point. Solution y₉ is obviously much more interesting because of the big vacancy between y₅ and y₆ on the first front.

The number of dominating points can be calculated efficiently: the selection candidates are compared with all solutions to check the dominance. This is possible in time O(mμ²), with m denoting the number of objectives and μ the population size. The expected runtime complexity is lower than the one of the basic SMS-EMOA. The worst case complexity of the SMS-EMOA is the same, because the population may consist of non-dominated solutions only. In this case, the d(s, P(t)) selection scheme is never applied. The average case complexity is assumed to be better.

2.3. Calculation of contributing hypervolume Δ_S${\Delta }_{\mathcal{S}}$

The best known algorithms calculate the hypervolume in a runtime, that is polynomial in the number of points, but grows exponentially with the number of objectives. An overview of such algorithms is given next. Additionally, dedicated algorithms to directly calculate all values of Δ_S${\Delta }_{\mathcal{S}}$ at once in case of two and three objectives are presented. These algorithms are significantly faster than repetitive calls of procedures to compute the hypervolume as a whole.

In case of two objectives, we take the points of the worst-ranked non-dominated front and sort them in ascending order according to the values of the first objective function f  ₁. We get a sequence that is additionally sorted in descending order concerning the f  ₂ values, because the points are mutually non-dominated. Given a sorted front R_v={s₁,…,s_|Rv|}${\mathcal{R}}_v=\{s_1\text{,}\dots \text{,}{s}_{|{\mathcal{R}}_v|}\}$, Δ_S${\Delta }_{\mathcal{S}}$ is calculated as follows (i=2,…,|R_v|-1$i=2\text{,}\dots \text{,}|{\mathcal{R}}_v|-1$):

equation(6)

Δ_S(s_i,R_v)=(f₁(s_i+1)-f₁(s_i))·(f₂(s_i-1)-f₂(s_i)).${\Delta }_{\mathcal{S}}(s_i\text{,}{\mathcal{R}}_v)=(f_1(s_{i+1})-f_1(s_i))·(f_2(s_{i-1})-f_2(s_i))\text{.}$

Turn MathJax on

A fast algorithm for computing all contributions in a three objective space is presented in [16]. It has a runtime of O(μ³) and is based on a two-dimensional projection of the set of points. The f₁–f₂ plane is divided into a grid according to the coordinates of all μ + 1 solution vectors and the reference point. For each grid cell c_ij (i = 1, … ,μ + 1, j = 1, … , μ + 1), the μ + 1 values of f₃ are computed at the corner point of the cell, which dominates it in the f₁–f₂ plane. Among these values, the best value h₁(c_ij) and second best value h₂(c_ij) are determined. Given this information, which can be stored in two matrices of dimension (μ + 1) × (μ + 1), the exclusive contribution of the ith solution vector (i = 1, … , μ + 1) can be calculated as follows. For any grid cell c_ij with h₁(c_ij) value determined by the ith solution, we determine the volume v_ij as the product of the surface area of c_ij and the difference ∣h₁(c_ij) − h₂(c_ij)∣. Accumulating all these values of v_ij yields to the exclusive contribution of the ith solution to the dominated hypervolume. The runtime is O(μ³) because the number of cells is quadratic and the calculations for each cell value are done in linear time. It is likely possible to further tune this algorithm, i.e. by using incremental updates. However, for the computations performed in this paper the algorithm performed reasonably fast. Another potential advantage is that the algorithm is non-recursive and uses simple data structures. Thus it is relatively easy to implement/debug.

An algorithm that computes the whole hypervolume of a set of points is the LebMeasure algorithms by Fleischer [5]. While [17] recently showed, that the runtime complexity is exponential in the number of objectives in the worst case, contrary to the runtime complexity mistakenly stated in [5].

The HSO algorithm (Hypervolume by Slicing Objectives) described by Knowles [18] and Zitzler [19] is exponential in the number of objectives as well. It is easier to implement, uses simpler data structures, and needs less space. Recently While et al. [20] presented heuristics for accelerating HSO. For a set of l points, HSO processes objectives one after the other and thereby divides the problem into l sets with a dimension reduced by one in each iteration. The runtime differs with the order in which the objectives are processed. The heuristics are applied as a pre-processing of HSO and calculate a promising order of objectives which reliably avoids the worst possible case (but not reliably the worst case complexity) of HSO. While et al. [20] demonstrated for some DTLZ functions the significant runtime reduction of 25–98% compared to the basic HSO.

2.4. Comparison to other EMOA

2.4.1. ESP

Huband et al. [12] incorporated an approximation of the hypervolume contribution into a selection operator resembling SPEA2 [21], replacing the common nearest neighbour criterion. They involved the “product of the one-dimensional lengths to the next worse objective function value in the front for each objective” [12] as selection criterion in case of equal fitness values. In the two-objective case, this exactly corresponds to the S$\mathcal{S}$ metric contribution (cf. Eq. (6)). In case of more than two objectives, this procedure provides only a lower bound of the hypervolume value. The great advantage is the reduced runtime of O(k) with k being the number of solutions on the front. The selection is performed using a (μ + λ)-ES selection scheme, setting μ = λ = 100. Selection with equal numbers of parents and offsprings is very common in evolutionary multiobjective optimisation; though evidence for the advantage of this setting is missing up to now. As variation operator, the ESP algorithm employs usual operators from evolution strategies coupled with a mutation probability p. Huband et al. choose p = 1/n as the probability to apply mutation for each decision variable.

ESP has been compared to several well-known EMOA like NSGA-II and SPEA2 on a set of test functions (ZDT1 to ZDT4 and ZDT6 [22]). In [12], the authors account the probabilistic mutation operator to be responsible for the obtained quality of ESP.

2.4.2. IBEA

The indicator-based EA (IBEA), proposed by Zitzler and Künzli [13], is a general framework to incorporate quality indicators in EMOA. Next to a hypervolume based indicator, they also incorporated a binary additive ϵ-indicator. The fitness of a solution is based on the sum of the indicator values resulting from pairwise comparisons to all other solutions.

They utilised ordinary EMOA variation operators for different test functions, i.e. SBX recombination and polynomial mutation on the real valued test functions under investigation (ZDT6, KUR, DTLZ2 and DTLZ6). Again, a (μ + μ) selection scheme and a specially scaled value of the objective and indicator values were incorporated in the improved adaptive IBEA. Both proposed variants yielded results “significantly better than SPEA2 and NSGA-II with respect to the performance indicators” [13].

2.4.3. The AA_reduce scheme

The SMS-EMOA presented here nearly fits into the generic algorithm scheme AA_reduce presented by Knowles and Corne [8] within the context of archiving strategies. The archiving strategy called AA_S uses the contributing hypervolume of the non-dominated points to determine the worst solution. Knowles and Corne showed that the AA_S strategy converges to a subset of the true Pareto front. Provided that the population size in the SMS-EMOA equals the archive size in AA_S and only non-dominated solutions are concerned, the AA_S strategy is equivalent to our selection method. If dominated solutions are considered as well, the SMS-EMOA population contains even more solutions than the AA_S archive. Thus, the proof of convergence holds for SMS-EMOA as well.

2.4.4. NSGA-II

The NSGA-II by Deb et al. [14] is an improved version of the preceding NSGA. It features non-dominated-sorting and a (μ   + μ  ) selection scheme. The secondary ranking criterion for solutions on the same front is called crowding distance. Its task corresponds to Δ_S${\Delta }_{\mathcal{S}}$ in the SMS-EMOA. For a solution, the crowding distance is defined as the sum of the side lengths of the cuboid through the neighbouring solutions on the front. For extremal solutions it is defined as infinity. The crowding distance value of a solution depends on its neighbours and not directly on the position of the point itself, in contrast to the Δ_S${\Delta }_{\mathcal{S}}$ measure.

The crowding distance is meant to distribute solutions uniformly on the fronts. In contrast to this, the hypervolume measure is meant to distribute them in a way that approximately maximises the covered hypervolume. This indicates that good compromise solutions, which are located in parts of the front with high curvature, are valued higher in SMS-EMOA than in NSGA-II.

3. Studies on artificial testfunctions

3.1. Distribution of solutions

A study was performed to get an impression of how the SMS-EMOA distributes solutions on Pareto fronts of different curvature, e.g. convex, concave and linear Pareto fronts. We conducted a study on simple two-objective test functions with a high dimensional decision space. The following family of generic functions, named EBN, is devised:

equation(7)

$f_1(\mathbf{x})≔{\left(\sum _{i=1}^n|x_i|\right)}^{\gamma }{n}^{-\gamma }\text{,}{f}_2(\mathbf{x})≔{\left(\sum _{i=1}^n|x_i-1|\right)}^{\gamma }{n}^{-\gamma }\text{,}\mathbf{x}\in [0\text{,}1{]}^n\text{,}$

Turn MathJax on

equation(8)

$y_2=(1-y_1^{(1/\gamma )}{)}^{\gamma }\text{.}$

Turn MathJax on

All points in [0, 1]ⁿ are Pareto-optimal solutions [23].

Fig. 2 shows that the solutions are not uniformly distributed on the Pareto front, but distributed in a way that maximises the covered hypervolume. The results demonstrate that the SMS-EMOA favours solutions in regions where the Pareto front has knee-points and captures the regions with fair trade-offs between different objectives. The regions with unbalanced trade-offs, located on the flanks of the Pareto front, are covered with less density, although extremal solutions are maintained. On the linear Pareto front, the points get uniformly distributed. In case of a concave Pareto front the regions with fair trade-offs are emphasised. These are located near the angular point of the Pareto front. This serves the practitioner who is mainly interested in a limited number of solutions on the Pareto front to get an overview of possible compromise solutions.

3.2. Settings for benchmarks

The SMS-EMOA was tested on several test problems from literature. We aimed at comparability to frequently cited papers of Deb et al. presenting the ϵ-MOEA approach [9] and [10]. Therefore, we also invoked the same variation operators, namely SBX recombination and polynomial mutation with the parametrisation given in [9] and [10]. The test problems named ZDT1 to ZDT4 and ZDT6 from [10] and [22] have been considered and are presented in Section 3.3, as well as the functions DTLZ1 to DTLZ4 [24] (see Section 3.4). This way, we compared our SMS-EMOA to NSGA-II, C-NSGA-II, SPEA2, and ϵ-MOEA.

The same parameter settings as in the benchmark from Deb et al. are invoked [9] and [10]. 20,000 function evaluations are calculated for the ZDT functions and 30,000 on all DTLZ functions except for DTLZ3, where 100,000 function evaluations are used. The population size is set to 100 for all functions, unlike in [10] where a size of 200 individuals was applied for DTLZ3. For DTLZ4, ten runs are executed and five for all other DTLZ and ZDT functions. Like in [9], only ‘successful’ runs are used to calculated the metric values. A run is termed successful, if the solution set is distributed over the whole true Pareto front. The quality of the algorithms is compared considering the S$\mathcal{S}$ metric and another performance measure, the convergence metric like in [9] and [10]. For every point in the solution set, the euclidean distance to the nearest point on the Pareto front is measured in the solution space. The convergence metric is the arithmetic mean of these distances. It is only applicable if the Pareto front is known a priori, as it is the case for many benchmark problems. The convergence metric does only account for the precision of the approximation to the Pareto front but not for the distribution of points. Thus, it is recommended to use it in combination with other measures [1].

For a clear overview, the best results for the S$\mathcal{S}$ metric and the convergence achieved in [10] and [9] have been copied to Table 1 and Table 2.

3.3. Two-dimensional functions (ZDT1 to ZDT4 and ZDT6)

The following results are depicted in Table 1. The results of the SMS-EMOA with the modified selection operator using the number of dominating points (cf. Section 2.2) are given in the lines named ‘SMS-EMOA dp’. The results of the two SMS-EMOA variants are very similar and they often share the same rank. This is a very positive result since the quality of the SMS-EMOA is maintained though algorithmic complexity and runtime have been reduced. In the following, we start addressing the achievements for each test functions separately, but do not further distinguish the results for the SMS-EMOA variants.

ZDT1 has a smooth convex Pareto front where the SMS-EMOA is ranked best on the S$\mathcal{S}$ metric and near to the best concerning the convergence measure. ZDT4 provides multiple parallel local fronts and the Pareto front is equivalent to that of ZDT1. On the basis of the given values from [9] and [10], we assume that all algorithms jumped across the second best front with most solutions and aimed at the best front, like SMS-EMOA. The worse ranked values of the other algorithms seem to stem from disadvantageous distributions. ZDT2 has a smooth concave front and the SMS-EMOA covers most hypervolume. ZDT3 has a discontinuous Pareto front that consists of five slightly convex parts. Here, the SMS-EMOA is a little bit better concerning the S$\mathcal{S}$ metric than the second ranked ϵ-MOEA and significantly better with respect to convergence. ZDT6 has a concave Pareto front that is equivalent to that of ZDT2, except for the differences that the front is truncated to a smaller range and that points are non-uniformly spaced. Here, the SMS-EMOA is ranked second on both measures, outperformed by SPEA2, which in contrast shows apparently bad results on the other more simple problems.

The SMS-EMOA is ranked best concerning the S$\mathcal{S}$ metric in all functions except for ZDT6. Concerning the convergence measure, it has two first ranks, two second ranks and one third rank. According to the sum of ranks of the two measures on each function, one can state that the SMS-EMOA provides best results on all considered functions, except for ZDT6, where it is outperformed by SPEA2. Comparing the sum of the achieved ranks of each measure shows that our algorithm obtains best results concerning both the convergence measure (with 9) and the S$\mathcal{S}$ metric (with 6). In order to test the significance of the results a Student t  -test has been applied. The results of the SMS-EMOA were tested for significance to all strategies that performed worse. The results with respect to the S$\mathcal{S}$ metric are significant with a significance level of more than 90%. Except for the ZDT6 test problem, even a higher significance level of 95% holds. Also for the convergence measure, the SMS-EMOA is better than the next best strategy with a significance level of more than 90%, except for ZDT4, where the comparison to CNSGA-2 and ϵ-EMOA does not lead to a sufficient significance level.

3.4. Three-dimensional functions (DTLZ1 to DTLZ4)

The results for the DTLZ functions are presented in Table 2. For easier comparison, the best values of those published in [9] are copied into the table. The abbreviation ‘NC’ stands for ‘not computed’, indicating missing values in the literature [9].

The results indicate that both SMS-EMOA variants clearly outperform the aforementioned EMOA concerning the S$\mathcal{S}$ metric value. Nearly the same holds for the convergence, where only the ϵ-MOEA applied to DTLZ1 is better. Just like on the ZDT functions, the comparison of the SMS-EMOA variants shows no significant differences. This confirms that the hypervolume-based selection is successful within different frameworks.

The SMS-EMOA solutions for functions DTLZ1 and DTLZ2 are visualised in Fig. 3. The median run concerning the S$\mathcal{S}$ metric values is shown, respectively. As there are no observably differences between the selection variants, only results received with the basic SMS-EMOA are displayed.

The Pareto front of DTLZ1 is a beveled, triangular plane in the first quadrant. The solutions of the SMS-EMOA variants lay exactly on this plane. They are distributed uniformly and extremal solutions at the edges of the plane are contained (cf. Fig. 3). The functions DTLZ2 to DTLZ4 share the same Pareto front, which is the eighth of the unit sphere, bounded to the first quadrant. There are no major differences between the results of the SMS-EMOA variants, again. Fig. 3 shows that the solutions for DTLZ2 are not equally distributed over the whole Pareto front. An adequate amount of boundary solutions are positioned along the eighth sphere. The remaining solutions are uniformly distributed over the inner part of the sphere surface, leaving an obvious gap to the boundary solutions. The function DTLZ3 has several local fronts parallel to the global one. All runs of the SMS-EMOA variants reached the global Pareto front and their convergence values are substantially better than the best EMOA studied in previous work (cf. Table 2). The distribution of solutions corresponds to the one on DTLZ2. The difficulty for function DTLZ4 is that the first two object variables are raised to the power of 100. Hence, there is a bias towards the special solutions resulting in sets located solely in the f₁–f₂-plane or the f₁–f₃-plane. Just like in [10], only the successful runs are considered for the metric values and their number out of 10 performed runs is given in brackets. The results of the successful runs are similar to those of the solutions of DTLZ2 and DTLZ3.

On almost all DTLZ functions the obtained results for the convergence metric and S$\mathcal{S}$ metric are significantly better for the SMS-EMOA than for the other approaches (with a significance level of 99% using a Student t  -test). Only on DTLZ1, the significance of the better S$\mathcal{S}$ metric values is not given and the ϵ-MOEA achieves better convergence values than the SMS-EMOA.

In conjunction, concerning this bundle of test problems, the SMS-EMOA can be regarded as the best evolutionary multiobjective algorithm. The outstanding performance concerning the S$\mathcal{S}$ metric is a very encouraging result even though good results seem to be natural because of the use of the S$\mathcal{S}$ metric as selection criterion and quality measure. One should admit that our approach is a rather simple one with only one population and only a few parameters. It is steady-state, resulting in a low selection pressure per generation. Neither there are any special variation operators fitted to the selection strategy, nor it is tuned for performance in any way.

The results in the convergence measure are even more surprising. Especially on the function that are supposed to be more difficult, the SMS-EMOA achieves very good results. A possible explanation might be that a population of well-distributed points is able to sample individuals with larger improvement. Further investigations are required to clarify this topic.

4. Application to airfoil design optimisation

Airfoil design applications are frequently addressed multiobjective real-world problems. The work on these applications is pushed by aeronautic industry because new solutions promise huge potential savings as well as better and more secure flight characteristics. Moreover, research in aeronautics receives governmental support in order to push up-to-date research in the related countries.

The multiobjective point of view is introduced by contradicting demands made to airfoils. On the one hand, it should provide a maximum lift to ensure save and stable flight qualities and a low energy consumption during the starting phase. On the other hand, a minimum drag is wanted during the cruising flight allowing for minimal energy consumption. Many more qualities influence the complex interplay of parameters to ensure best flight characteristics.

In this paper, we exemplarily deal with two airfoil design test cases. A computational fluid dynamics (CFD) tool based on the solutions of Navier–Stokes equations computes the properties of the airfoils proposed by the optimisation technique. The parameterisation of the airfoils is done using Bezier-splines. One computation of an objective function value by time consuming CFD tools typically takes several minutes or more. Hence, only a limited number of evaluations can be afforded. We allow 1000 evaluations to stay comparable to previous studies on these test problems.

An airfoil redesign test case is the application investigated firstly. The two objectives defined are the differences in pressure distribution between target airfoils that are nearly optimal for given flow conditions and the airfoil proposed by the optimisation technique. The target airfoils are NACA 0012 and NACA 4412, why we call this one the ‘NACA test case’. These flow conditions correspond to the take-off phase (high lift required) and the cruising phase (low drag desired) during an aircraft flight and can be taken from the left hand part of Table 3. For both flow conditions, the pressure distribution around the airfoil is calculated to derive the differences to the two target airfoils. The integrated differences serve as fitness function values for the test case. Due to the parameterisation in use, this test case features 18 degrees of freedom. A more detailed description can be found in [7].

Within the second test case, three drag coefficients are to be minimised. The considered airfoil type is the RAE 2822 airfoil, hence we will term the corresponding test case ‘RAE test case’. The different drag coefficients stem from three different flow conditions, each yielding different values for drag, lift and pitching moment. The conditions can be taken from the right hand part of Table 3. The task is to minimise the drag coefficients while not losing lift and keeping the pitching moment within a 2% range. Furthermore, geometrical constraints have to be met for the thickness of the airfoil measured at two particular points: the maximum thickness, its leading edge radius, and its trailing edge angle.

The geometrical information about a proposed airfoil can be received from the evaluation software just after the airfoil shape is generated. The whole time-consuming procedure of grid generation, solving the flow describing equations, and all post-processing tasks are not required to receive this information. Therefore, the geometrical constraints are treated differently from the aerodynamic ones, which require the costly evaluation. Variation operators are applied up to 1000 times unless a feasible solution subject to the geometrical constraints is obtained. This was done in order to achieve a higher ratio of feasible individuals. If no feasible solution is found after 1000 samples, the geometrical constraints will be treated like all other constraints. Due to the parameterisation in use here, the RAE test case features six degrees of freedom.

4.1. Fitness function approximations

We use Kriging metamodels as described by Sacks et al. [25] as fitness function approximation tools (FFA) to accelerate the SMS-EMOA in the presence of time consuming objective function evaluations. For a new design point x∈S$\mathbf{x}\in \mathbb{S}$, the Kriging metamodel allows for a prediction ${\hat{y}}_i(\mathbf{x})(i\in \{1\text{,}\dots \text{,}m\})$ of the objective function value f  _i(x) from the set of previously evaluated (neighbouring) points stored in a database. Moreover, it estimates the uncertainty of the prediction by means of a standard deviation ${\hat{s}}_i(\mathbf{x})$.

In each generation, the metamodel-assisted SMS-EMOA   [7] generates a surplus of individuals by mutations of one chosen individual. From these individuals, we retrieve the most promising solution by means of the metamodel. The chosen individual is then evaluated exactly using the precise evaluation function and considered in the (μ   + 1)-selection. An optimistic guess of the objective function values, ${\hat{y}}_i(\mathbf{x})-2{\hat{s}}_i(\mathbf{x})\text{,}(i=1\text{,}\dots \text{,}m)$ serves as a ranking criterion in the selection. The eligible individual is the one with the potentially highest value of Δ_S${\Delta }_{\mathcal{S}}$ [7] based on this optimistic guess. This selection criterion will be termed lower bound criterion in the following. Another variant is to apply simply the predicted function value for estimating the contribution in hypervolume, thereby neglecting its associated uncertainty measure. The latter strategy will be referred to as mean value criterion. ²

4.2. Results

The SMS-EMOA provided very good and encouraging results on the theoretical test problems as well as on the design optimisation tasks. For the current study on airfoils, we collected five runs for each setting.

4.2.1. NACA results

Here, we considered SMS-EMOA without fitness function approximation (FFA) as well as the metamodel-assisted SMS-EMOA with mean value and lower bound criterion as described above. The points of the sets achieved using fitness function approximations are much better distributed than the ones obtained without. While runs without FFA are always biased towards one special region in the solution space, runs utilising FFA do not show any special crowding. The solutions are more equally distributed all over the Pareto front, with a desired higher density in regions with fair trade-offs as discussed above. The reason is the huge amount of pre-evaluations that are used to find promising regions of the search space to locate individuals that are to be evaluated exactly. Compared to the results with FFA, the runs without FFA do not seem to tap into their full potential probably due to the fact of too few evaluations.

The resulting sets without FFA are the worst fronts of all presented ones except for one special region. In the other regions, the SMS-EMOA with lower bound criterion seems to be better than the other algorithms shortly followed by old results from NSGA-II with lower bound criterion. This is the best result from [26] where different EMOA with and without FFA were studied. We therefore conclude, that the SMS-EMOA with lower bound criterion yields better results than all EMOA tested in [26].

The SMS-EMOA with mean value criterion provided the worst results utilising FFA within this study. The observation that the lower bound approximation technique offers better results than the mean value approximation was made in [26] as well. This seems to be a general achievement that deserves further attention.

4.2.2. RAE results

The results presented here have been obtained by the basic SMS-EMOA and are compared with results from an earlier publication, where a standard NSGA-II implementation was considered [27]. NSGA-II (with SBX, polynomial mutation and a population size of 20 individuals) was the best algorithm in that publication. We therefore conclude, that algorithms outperforming the NSGA-II results, provide better results than all other algorithmic variants studied in [27].

The NSGA-II results (termed NSGA-II-SBP-20 in Fig. 4) are compared to results from SMS-EMOA featuring the same variation operators and the same population size (SMS-SBP-20). Furthermore, SMS-EMOA was tested with standard variation operators from evolution strategies.³ This setting was tested with a population size of 10 and 20 individuals (SMS-ES-10, SMS-ES-20), respectively.

Fig. 4 presents the non-dominated solutions from the union of all five received Pareto front approximations per setting in the left hand part. The right hand part gives a median attainment surface consisting of all points from a Pareto front approximation that are weakly dominated by at least two of the other Pareto front approximations obtained with the same strategy.

In the top row of Fig. 4, the projections of the three-dimensional solutions into the f₁–f₂-plane are shown. The results demonstrate, that the fitness function components f₁ and f₂ are highly correlated to each other. Nearly all solutions appear to improve the baseline design (marked with a circle in the upper right corner). Furthermore, solutions received from one strategy seem to dominate each other. This is an artifact caused by the projection of the three-dimensional points. The high correlation between the fitness function components f₁ and f₂ causes the projections in planes f₁–f₃ (middle row) and f₂–f₃ (bottom row) to look similar. In the f₁–f₃-plane (middle row in Fig. 4), the variance in the results received from NSGA-II is lower than the ones from SMS-EMOA. The results of SMS-EMOA show a wider distribution, which is strictly in accordance with the intentions followed within the development of the algorithm.

From the right hand part of Fig. 4 (median attainment surface) it is hardly possible to distinguish the different results. The left hand part (the best obtained solutions) allow for a much better distinction. Here, a clear best setting can be recognised, which is SMS-EMOA, featuring ES variation operators and a population size of 10 individuals. The second best setting is the aformentioned one with a population size of 20 individuals. Both other strategies can be ranked equally with an advantage for the SMS-EMOA variant in distribution of solutions and quality in the third fitness function component (f₃). The NSGA-II behaves a bit better with respect to the other fitness function components (f₁ and f₂). Note, that SMS-ES-10 is significantly better than all other strategies with respect to all characteristics mentioned above.

The superiority of these results may be explained by the use of the ES variation operators in contrast to the SBX crossover and the polynomial mutation as well as the smaller population size. In concordance with our experience the SBX and polynomial mutation perform well on the theoretical test cases under investigation. Real-world problems, however, seems to be better addressed with the ES variation operators. The benefit of the smaller population size is detectable in all plots. This is probably due to the higher reproduction probability of each individual in smaller populations and the resulting higher selection pressure.

From an aerodynamic point of view, it is important to improve the baseline design with respect to all given flight conditions; this means finding solutions dominating the baseline design. Only one such result has been calculated until now, using the metamodel assisted NSGA-II as discussed in [27]. The ability to produce such good results is due to the fact that a large number of pre-evaluations is comprised. Our results show that SMS-EMOA is always able to simultaneously improve all objectives of the baseline design. This can only be explained by the favourable hypervolume selection incorporated in SMS-EMOA.

5. Summary and prospects

Multiple results on standard benchmark problems indicate that the SMS-EMOA is well-suited for Pareto optimisation with two and three objectives. It clearly outperforms established techniques such as SPEA2, ϵ  -MOEA, and NSGA-II, both in convergence and in S$\mathcal{S}$ metric values received. The examples reveal that results are well-distributed on the Pareto front. Boundaries of the Pareto front and regions around knee points are favoured. A distinguishing feature of the SMS-EMOA is that it is well-suited for approximating Pareto sets by a small number of individuals. This is often desired in practice.

A new selection criterion has been developed, which replaces the hypervolume-based selection in the presence of dominated solutions in the population. Using this new concept, the rating of a solution depends on the number of solutions dominating it. Hence, it focuses the evolutionary search towards less explored regions near the Pareto front. Both variants of the SMS-EMOA (the basic algorithm and the coupling with the number of dominated points criterion) led to better results than obtained by conventional strategies. This indicates that the hypervolume selection mechanism, performs well, regardless of the particular details of the selection scheme. This indicates a kind of robustness with regard to the special selection method chosen.

Besides the comparison on state-of-the-art benchmarks, the SMS-EMOA has been tested on challenging real-world applications, namely the optimisation of airfoils within two or three different flight conditions. Not only the distributions of solutions could be improved, but also solutions that clearly dominate the baseline design have been found. Additionally, SMS-EMOA was coupled with a metamodelling tool. This tool provides fast fitness function approximations to save costly exact evaluation in unfavourable regions. The SMS-EMOA was the first algorithm which always produced solutions dominating the baseline one. Moreover this algorithm produces the solutions robustly.

Future research should envisage the study of further algorithmic variants and a deepened analysis of strategy parameters. By now, the computational effort to compute the S$\mathcal{S}$ metric values is very high for more than three objectives. Therefore, the SMS-EMOA is hardly applicable to higher dimensional problems that normally desire a large number of function evaluations. Nevertheless, many real-world applications allow only a very limited number of evaluations due to costly simulations that govern the runtime of the optimisation process. In these cases, the runtime of the operations performed within the EMOA can almost be neglected, and the SMS-EMOA is certainly a well-suited optimiser.

References

Copyright © 2006 Elsevier B.V. All rights reserved.