A tutorial on multiobjective optimization: fundamentals and evolutionary methods

In almost no other field of computer science, the idea of using bio-inspired search paradigms has been so useful as in solving multiobjective optimization problems. The idea of using a population of search agents that collectively approximate the Pareto front resonates well with processes in natural evolution, immune systems, and swarm intelligence. Methods such as NSGA-II, SPEA2, SMS-EMOA, MOPSO, and MOEA/D became standard solvers when it comes to solving multiobjective optimization problems. This tutorial will review some of the most important fundamentals in multiobjective optimization and then introduce representative algorithms, illustrate their working principles, and discuss their application scope. In addition, the tutorial will discuss statistical performance assessment. Finally, it highlights recent important trends and closely related research fields. The tutorial is intended for readers, who want to acquire basic knowledge on the mathematical foundations of multiobjective optimization and state-of-the-art methods in evolutionary multiobjective optimization. The aim is to provide a starting point for researching in this active area, and it should also help the advanced reader to identify open research topics.


Introduction
Consider making investment choices for an industrial process. On the one hand the profit should be maximized and on the other hand environmental emissions should be minimized. Another goal is to improve safety and quality of life of employees. Even in the light of mere economical decision making, just following the legal constraints and minimizing production costs can take a turn for the worse.
Another application of multiobjective optimization can be found in the medical field. When searching for new therapeutic drugs, obviously the potency of the drug is to be maximized. But also the minimization of synthesis costs and the minimization of unwanted side effects are much-needed objectives (van der Horst et al. 2012;Rosenthal and Borschbach 2017).
There are countless other examples where multiobjective optimization has been applied or is recently considered as a promising field of study. Think, for instance, of the minimization of different types of error rates in machine learning (false positives, false negatives) (Yevseyeva et al. 2013;Wang et al. 2015), the optimization of delivery costs and inventory costs in logistics(Geiger and Sevaux 2011), the optimization of building designs with respect to health, energy efficiency, and cost criteria (Hopfe et al. 2012).
In the following, we consider a scenario where given the solutions in some space of possible solutions, the so-called decision space which can be evaluated using the so-called objective functions. These are typically based on computable equations but might also be the results of physical experiments. Ultimately, the goal is to find a solution on which the decision maker can agree, and that is optimal in some sense.
When searching for such solutions, it can be interesting to pre-compute or approximate a set of interesting solutions that reveal the essential trade-offs between the objectives. This strategy implies to avoid so-called Pareto dominated solutions, that is solutions that can improve in one objective without deteriorating the performance in any other objective. The Pareto dominance is named after Vilfredo Pareto, an Italian economist. As it was earlier mentioned by Francis Y.Edgeworth, it is also sometimes called Edgeworth-Pareto dominance (see Ehrgott 2012 for some historical background). To find or to approximate the set of non-dominated solutions and make a selection among them is the main topic of multiobjective optimization and multicriterion decision making. Moreover, in case the set of nondominated solutions is known in advance, to aid the decision maker in selecting solutions from this set is the realm of decision analysis (aka decision aiding) which is also part of multi-criterion decision making.
Definition 1 Multiobjective Optimization. Given m objective functions f 1 : X ! R; . . .; f m : X ! R which map a decision space X into R, a multiobjective optimization problem (MOP) is given by the following problem statement: minimize f 1 ðxÞ; . . .; minimize f m ðxÞ; x 2 X ð1Þ Remark 1 In general, we would demand m [ 1 when we talk about multiobjective optimization problems. Moreover, there is the convention to call problems with large m, not multiobjective optimization problems but many-objective optimization problems (see Fleming et al. 2005;Li et al. 2015). The latter problems form a special, albeit important case of multiobjective optimization problems.
Remark 2 Definition 1 does not explicitly state constraint functions. However, in practical applications constraints have to be handled. Mathematical programming techniques often use linear or quadratic approximations of the feasible space to deal with constraints, whereas in evolutionary multiobjective optimization constraints are handled by penalties that increase the objective function values in proportion to the constraint violation. Typically, penalized objective function values are always higher than objective function values of feasible solutions. As it distracts the attention from particular techniques in MOP solving, we will only consider unconstrained problems. For strategies to handle constraints, see Coello Coello (2013).
Considering the point(s) in time when the decision maker interacts or provides additional preference information, one distinguishes three general approaches to multiobjective optimization (Miettinen 2012): 1. A priori: A total order is defined on the objective space, for instance by defining a utility function R m ! R and the optimization algorithm finds a minimal point (that is a point in X ) and minimum value concerning this order. The decision maker has to state additional preferences, e.g., weights of the objectives, prior to the optimization. 2. A posteriori: A partial order is defined on the objective space R m , typically the Pareto order, and the algorithm searches for the minimal set concerning this partial order over the set of all feasible solutions. The user has to state his/her preferences a posteriori, that is after being informed about the trade-offs among nondominated solutions. 3. Interactive (aka Progressive): The objective functions and constraints and their prioritization are refined by requesting user feedback on preferences at multiple points in time during the execution of an algorithm.
In the sequel, the focus will be on a posteriori approaches to multiobjective optimization. The a priori approach is often supported by classical single-objective optimization algorithms, and we refer to the large body of the literature that exists for such methods. The a posteriori approach, however, requires interesting modifications of theorems and optimization algorithms-in essence due to the use of partial orders and the desire to compute a set of solutions rather than a single solution. Interactive methods are highly interesting in real-world applications, but they typically rely upon algorithmic techniques used in a priori and a posteriori approaches and combine them with intermediate steps of preference elicitation. We will discuss this topic briefly at the end of the tutorial.

Related work
There is a multiple of introductory articles that preceded this tutorial: • In  a tutorial on state-of-the-art evolutionary computation methods in 2004 is provided including Strength Pareto Evolutionary Algorithm Version 2 (SPEA2) (Zitzler et al. 2001), Non-dominated Sorting Genetic Algorithm II (NSGA-II) (Deb et al. 2002), Multiobjective Genetic Algorithm (MOGA) (Fonseca and Fleming 1993) and Pareto-Archived Evolution Strategy (PAES) (Knowles and Corne 2000) method. Indicator-based methods and modern variants of decomposition based methods, that our tutorial includes, were not available at that time. • In Deb (2008) an introduction to earlier multiobjective optimization methods is provided, and also in the form of a tutorial. The article contains references to early books in this field and key articles and also discusses applications.
• Derivative-free methods for multiobjective optimization, including evolutionary and direct search methods, are discussed in Custódio et al. (2012). • On conferences such as GECCO, PPSN, and EMO there have been regularly tutorials and for some of these slides are available. A very extensive tutorial based on slides is the citable tutorial by Brockhoff (2017).
Our tutorial is based on teaching material and a reader for a course on Multiobjective Optimization and Decision Analysis at Leiden University, The Netherlands (http:// moda.liacs.nl). Besides going into details of algorithm design methodology, it also discusses foundations of multiobjective optimization and order theory. In the light of recent developments on hybrid algorithms and links to computational geometry, we considered it valuable to not only cover evolutionary methods but also include the basic principles from deterministic multiobjective optimization and scalarization-based methods in our tutorial.

Order and dominance
For the notions we discuss in this section a good reference is Ehrgott (2005). The concept of Pareto dominance is of fundamental importance to multiobjective optimization, as it allows to compare two objective vectors in a precise sense. That is, they can be compared without adding any additional preference information to the problem definition as stated in Definition 1.
In this section, we first discuss partial orders, pre-orders, and cones. For partial orders on R m there is an important geometric way of specifying them with cones. We will define the Pareto order (aka Edgeworth-Pareto order) on R m . The concept of Pareto dominance is of fundamental importance for multiobjective optimization, as it allows to compare two objective vectors in a precise sense (see Definition 5 below). That is, comparisons do not require adding any additional preference information to the problem definition as stated in Definition 1. This way of comparison establishes a pre-order (to be defined below) on the set of possible solutions (i.e., the decision space), and it is possible to search for the set of its minimal elements-the efficient set.
As partial orders and pre-orders are special binary relations, we digress with a discussion on binary relations, orders, and pre-orders.
Remark 3 Sometimes we will also write xRy for ðx; yÞ 2 R.

Now we can define different types of orders:
Definition 3 Pre-order, Partial Order, Strict Partial Order. A binary relation R is said to be a -pre-order (aka quasi-order), if and only if it is transitive and reflexive, -partial order, if and only if it is an antisymmetric preorder, -strict partial order, if and only if it is irreflexive and transitive Remark 4 Note that a strict partial order is necessarily asymmetric (and therefore also anti-symmetric).
Proposition 1 Let X be a set and D ¼ fðx; xÞjx 2 Xg be the diagonal of X.
-If R is an anti-symmetric binary relation on X, then any subset of R is also an anti-symmetric binary relation. -If R is irreflexive, then (R is asymmetric if and only if R is antisymmetric). Or: the relation R is asymmetric if and only if R is anti-symmetric and irreflexive. -If R is a pre-order on X, then fðx; yÞ j ðx; yÞ 2 R and ðy; xÞ 6 2 Rg, denoted by R strict , is transitive and irreflexive. In other words, R strict is a strict partial order associated to the pre-order R.
-If R is a partial order on X, then R n D is irreflexive and transitive. In other words, R n D is a strict partial order. Moreover R n D is anti-symmetric (or asymmetric). -If R is a pre-order on X, then (R n D is a strict partial order if and only if R is asymmetric).
Remark 5 In general, if R is a pre-order, then R n D does not have to be transitive. Therefore, in general, R n D will not be a strict partial order.
Definition 4 Minimal Element. A minimal element x 2 X in a (strictly) partially ordered set (X, R) is an element for which there does not exist an x 0 2 X with x 0 Rx and x 0 6 ¼ x.
(In case, the order R is a strict partial order, x 0 Rx implies x 0 6 ¼ x). In words, in case that y ð1Þ 0 Pareto y ð2Þ the first vector is not worse in each of the objectives and better in at least one objective than the second vector.
Proposition 2 The Pareto order 0 Pareto on the objective space R m is a strict partial order. Moreover ð0 Pareto [ DÞ is a partial order. We denote this by " Pareto or also by " if the context provides enough clarity.
In multiobjective optimization we have to deal with two spaces: The decision space, which comprises all candidate solutions, and the objective space which is identical to R m and it is the space in which the objective function vectors are represented. The vector-valued function f ¼ ðf 1 ; . . .; f m Þ > maps the decision space X to the objective space R m . This mapping and the Pareto order on R m as defined in Definition 5 can be used to define a pre-order on the decision space X as follows.
Definition 6 Pre-order on Search Space. Let x 1 ; x 2 2 X.
The solution x 1 is said to Pareto dominate the solution x 2 if and only if fðx 1 Þ 0 Pareto fðx 2 Þ. Notation: x 1 Pareto dominates x 2 is denoted by x 1 0 f x 2 .

Remark 6
The binary relation 0 f on X is a strict partial order on X and ð0 f [fðx; xÞ j x 2 XgÞ is a partial order on X . Note that the pre-order R associated to " Pareto via f ( i.e., x 1 Rx 2 if and only if fðx 1 Þ " Pareto fðx 2 Þ ) is, in general, not asymmetric and therefore, in general, 0 f 6 ¼ R n fðx; xÞ j x 2 Xg.
Sometimes we need the notion of the so called strict component order on R m and its accompanying notion of weak non-dominance.

Cone orders
The Pareto order is a special case of a cone order, which are orders defined on vector spaces. Defining the Pareto order as a cone order gives rise to geometrical interpretations. We will introduce definitions for R m , although cones can be defined on more general vector spaces, too. The binary relations in this subsection are subsets of R m Â R m and the cones are subsets of R m .
Definition 9 Non-trivial Cone. A set C & R m with ; 6 ¼ C 6 ¼ R m is called a non-trivial cone, if and only if 8a 2 R; a [ 0; 8c 2 C : ac 2 C.

Remark 8
In the sequel when we say cone we mean nontrivial cone.
Definition 10 Minkowski Sum. The Minkowski sum (aka algebraic sum) of two sets A 2 R m and B 2 R m is defined as A È B :¼ fa þ b j a 2 A^b 2 Bg. Moreover we define aA ¼ faaj a 2 Ag. Definition 11 The binary relation, R C , associated to the cone C. Given a cone C the binary relation associated to this cone, notation R C , is defined as follows: 8x 2 R m : 8y 2 R m : ðx; yÞ 2 R C if and only if y 2 fxg È C.
Remark 10 It is clear that for any cone C the associated binary relation is translation invariant (i.e, if 8u 2 R m : ðx; yÞ 2 R C ) ðx þ u; y þ uÞ 2 R C ) and also multiplication invariant by any positive real (i.e., 8a [ 0 : ðx; yÞ 2 R C ) ðax; ayÞ 2 R C ). Conversely, given a binary relation R which is translation invariant and multiplication invariant by any positive real, the set C R :¼ fy À x j ðx; yÞ 2 Rg is a cone. The above two operations are inverses of each other, i.e., to a cone C one associates a binary relation R C which is translation invariant and multiplication invariant by any positive real, and the associated cone of R C is C, and conversely starting from a binary relation R which is translation invariant and multiplication invariant by any positive real one obtains the cone C R and the binary relation associated to this cone is R. In short, there is a natural one to one correspondence between cones and translation invariant and multiplication-invariant-bypositive-reals binary relations on R m .
Note that for a positive multiplication invariant relation R the set C R ¼ fy À x j xRy g is a cone. We restrict our attention to relations which are translation invariant as well in order to get the above mentioned bijection between cones and relations.
Also note if a positive multiplication invariant and translation invariant relation R is such that ; 6 ¼ R 6 ¼ R m Â R m , then the associated cone C R is nontrivial. Relations associated to non-trivial cones are nonempty and not equal to all of R m Â R m .

Remark 11
In general the binary relation R C associated to a cone is not reflexive nor transitive nor anti-symmetric. For instance, the binary relation R C is reflexive if and only if 0 2 C. The following definitions are needed in order to state for which cones the associated binary relation is antisymmetric and/or transitive.
Definition 12 Pointed cone and convex cone. A cone C is pointed if and only if C \ ÀC f0g where ÀC ¼ fÀc j c 2 Cg and C is convex if and only if 8c 1 2 C; c 2 2 C; 8a such that 0 a 1 : ac 1 þ ð1 À aÞc 2 2 C.
With these definitions we can specify for which cones the associated relation is transitive and/or anti-symmetric: Proposition 3 Let C be a cone and R C its associated binary relation (i.e., R C ¼ fðx; yÞ j y À x 2 Cg) . Then the following statements hold.
-R C is translation and positive multiplication invariant, -R C is anti-symmetric if and only if C is pointed, -R C is transitive if and only if C is convex, and moreover, -R C is reflexive if and only if 0 2 C.
A similar statement can be made if we go in the other direction, i.e.: Proposition 4 Let R be a translation and positive multiplication invariant binary relation and the C R the associated cone (i.e., C R ¼ fy À x j ðx; yÞ 2 Rg). Then the following statements hold.
-C R is a cone, -R is anti-symmetric if and only if C R is pointed, -R is transitive if and only if C R is convex, and moreover, -R is reflexive if and only if 0 2 C R .
In the following definition some important subsets in R m ; m ! 1 are introduced.
Definition 13 Let m be a natural number bigger or equal to 1. The non-negative orthant (aka hyperoctant) of R m , denoted by R m ! 0 is the set of all elements in R m whose coordinates are non-negative. Furthermore, the zero-dominated orthant, denoted by R m 10 , is the set R m ! 0 n f0g. Analogously we define the non-positive orthant of R m , denoted by R 0 , as the set of elements in R m the coordinates of which are non-positive. Furthermore, the set of elements in R m which dominate the zero vector 0, denoted by R m 00 , is the set R m 0 n f0g. The set of positive reals is denoted by R [ 0 and the set of non-negative reals is denoted by R ! 0 .
Remark 12 The sets defined in the previous definition are cones.
Proposition 5 The Pareto order 0 Pareto on R m is given by the cone order with cone R m 10 , also referred to as the Pareto cone.
Remark 13 As R m 10 is a pointed and convex cone, the associated binary relation is irreflexive, anti-symmetric and transitive (see Proposition 3). Of course, this can be verified more directly.
The reason to view the Pareto order as derived from a cone is that it gives the opportunity to study this order more geometrically. For instance, the definition and many of the properties of the very important hypervolume indicator (to be defined later) readily become intuitive. A reason for deviating from the Pareto cone could be to add constraints to the trade-off between solutions. Moreover, see later for a discussion, the more general cones turned out to be very useful in generalizing the hypervolume indicator and influence the distribution of points in the approximation set to the Pareto front.
Alternatives to the standard Pareto order on R m can be easily imagined and constructed by using pointed, convex cones. The alternatives can be used, for instance, in preference articulation.

Time complexity of basic operations on ordered sets
Partial orders do not have cycles. Let R be a partial order. It is easy to see that R does not have cycles. We show that the associated strict partial order does not have cycles. That is, there do not exist where D is the diagonal. For suppose such b i ; i ¼ 1; Á Á Á ; t À 1 can be found with this property. Then by transitivity of R n D (see Proposition 1), we get ðb 1 ; b tÀ1 Þ 2 R n D. By assumption, we have ðb tÀ1 ; b 1 Þ 2 R n D. Again by transitivity, we get ðb 1 ; b 1 Þ 2 R n D which is a contradiction. In other words, R does not have cycles. (The essence of the above argument is, that any strict partial order does not have cycles.) The absence of cycles for (strict) partial orders gives rise to the following proposition.
Proposition 6 Let S be a (strict) partially ordered set. Then any finite, non-empty subset of S has minimal elements (with respect to the partial order). In particular, any finite, non-empty subset Y & R m has minimal elements with respect to the Pareto order 0 Pareto . Also any, finite nonempty subset X & X has minimal elements with respect to 0 f . The question arises: How fast can the minimal elements be obtained?
Proposition 7 Given a finite partially ordered set ðX; "Þ, the set of minimal elements can be obtained in time Hðn 2 Þ.
Proof A double nested loop can check non-domination for each element. For the lower bound consider the case that all elements in X are incomparable. Only in this case is X the minimal set. It requires time Xðn 2 Þ to compare all pairs (Daskalakis et al. 2011). h Fortunately, in case of the Pareto ordered set ðX; 0 Pareto Þ, one can find the minimal set faster. The algorithm suggested by Kung et al. (1975) combines a dimension sweep algorithm with a divide and conquer algorithm and finds the minimal set in time Oðnðlog nÞÞ for d ¼ 2 and in time Oðnðlog nÞ dÀ2 Þ for d ! 3. Hence, in case of small finite decision spaces, efficient solutions can be identified without much effort. In the case of large combinatorial or continuous search spaces, however, optimization algorithms are needed to find them.

Scalarization techniques
Classically, multiobjective optimization problems are often solved using scalarization techniques (see, for instance, Miettinen 2012). Also in the theory and practice of evolutionary multiobjective optimization scalarization plays an important role, especially in the so-called decomposition based approaches.
In brief, scalarization means that the objective functions are aggregated (or reformulated as constraints), and then a constrained single-objective problem is solved. By using different parameters of the constraints and aggregation function, it is possible to obtain different points on the Pareto front. However, when using such techniques, certain caveats have to be considered. In fact, one should always ask the following two questions: 1. Does the optimization of scalarized problems result in efficient points? 2. Can we obtain all efficient points or vectors on the Pareto front by changing the parameters of the scalarization function or constraints?
We will provide four representative examples of scalarization approaches and analyze whether they have these properties.

Linear weighting
A simple means to scalarize a problem is to attach nonnegative weights (at least one of them positive) to each objective function and then to minimize the weighted sum of objective functions. Hence, the multiobjective optimization problem is reformulated to: Definition 14 Linear Scalarization Problem. The linear scalarization problem (LSP) of an MOP using a weight vector w 2 R m 10 , is given by Proposition 8 The solution of an LSP is on the Pareto front, no matter which weights in R m 10 are chosen. Proof We show that the solution of the LSP cannot be a dominated point, and therefore, if it exists, it must necessarily be a non-dominated point. Consider a solution of the LSP against some weights w 2 R m 10 , say x Ã and suppose this minimal point is dominated. Then there exists an objective vector y 2 fðX Þ with 8i 2 f1; . . .; mg y i f i ðx Ã Þ and for some index j 2 f1; . . .; mg it holds that y j \f j ðx Ã Þ.
Hence, it must also hold that P m i¼1 w i y i \ P m i¼1 w i f i ðx Ã Þ, which contradicts the assumption that x Ã is minimal. h In the literature the notion of convexity (concavity) of Pareto fronts is for the most part not defined formally. Possible formal definitions for convexity and concavity are as follows.
Proposition 9 In case of a convex Pareto front, for each solution in Y N there is a solution of a linear scalarization problem for some weight vector w 2 R m 10 .
If the Pareto front is non-convex, then, in general, there can be points on the Pareto front which are the solutions of no LSP. Practically speaking, in the case of concave Pareto fronts, the LSP will tend to give only extremal solutions, that is, solutions that are optimal in one of the objectives. This phenomenon is illustrated in Fig. 2, where the tangential points of the dashed lines indicate the solution obtained by minimizing an LSP for different weight choices (colors). In the case of the non-convex Pareto front (Fig. 2, right), even equal weights (dark green) cannot lead to a solution in the middle part of the Pareto front. Also, by solving a series of LSPs with minimizing different weighted aggregation functions, it is not possible to obtain this interesting part of the Pareto front.

Chebychev scalarization
Another means of scalarization, that will also uncover points in concave parts of the Pareto front, is to formulate the weighted Chebychev distance to a reference point as an objective function.
Definition 17 Chebychev Scalarization Problem. The Chebychev scalarization problem (CSP) of an MOP using a weight vector k 2 R m 10 , is given by minimize max i2f1;...;mg where z Ã is a reference point, i.e., the ideal point defined as Proposition 10 Let us assume a given set of mutually nondominated solutions in R m (e.g., a Pareto front). Then for every non-dominated point p there exists a set of weights for a CSP, that makes this point a minimizer of the CSP provided the reference point z Ã is properly chosen (i.e., the vector p À z Ã either lies in the positive or negative orthant).
Practically speaking, Proposition 10 ensures that by changing the weights, all points of the Pareto front can, in principle, occur as minimizers of CSP. For the two example Pareto fronts, the minimizers of the Chebychev scalarization function are points on the iso-height lines of the smallest CSP function value which still intersect with the Pareto front. Clearly, such points are potentially found in convex parts of Pareto fronts as illustrated in Fig. 3 (left) as well as in concave parts (right).However, it is easy to construct examples where a CSP obtains minimizers in weakly dominated points (see Definition 8). Think for instance of the case fðX Þ ¼ ½0; 1 2 . In this case all points on the line segment ð0; 0Þ > ; ð0; 1Þ > and on the line segment ð0; 0Þ > ð1; 0Þ > are solutions of some Chebychev scalarization. (The ideal point is 0 ¼ ð0; 0Þ > , one can take as weights (0, 1) for the first scalarization, and (1, 0) for the second scalarization; the Pareto front is equal to fð0; 0Þ > g).
In order to prevent this, the augmented Chebychev scalarization provides a solution. It reads: where is a sufficiently small, positive constant.

-constraint method
A rather straightforward approach to turn a multiobjective optimization problem into a constraint single-objective optimization problem is the -constraint method.
The method is illustrated in Fig. 4

Boundary intersection methods
Another often suggested way to find an optimizer is to search for intersection points of rays with the attained subset fðX Þ (Jaszkiewicz and Słowiński 1999). For this method, one needs to choose a reference point in R m , say r, which, if possible, dominates all points in the Pareto front. Alternatively, in the Normal Boundary Intersection method (Das and Dennis 1998) the rays can emanate from a line (in the bi-objective case) or an m À 1 dimensional hyperplane, in which case lines originate from different evenly spaced reference points (Das and Dennis 1998). Then the following problem is solved: Definition 19 Boundary Intersection Problem. Let d 2 R m 10 denote a direction vector and r 2 R m denote the reference vector. Then the boundary intersection problem is formulated as: subject to ðaÞ r þ td À fðxÞ ¼ 0; ðbÞ x 2 X; and ðcÞ t 2 R ! 0 Constraints (a) and (b) in the above problem formulation enforce that the point is on the ray and also that there exists a pre-image of the point in the decision space. Because t is minimized, we obtain the point that is closest to the reference point in the direction of d. This method allows some intuitive control on the position of resulting Pareto front points. Excepting rare degenerate cases, it will obtain points on the boundary of the attainable set fðX Þ. However, it also requires an approximate knowledge of the position of the Pareto front. Moreover, it might result in dominated points if the Pareto front is not convex. The method is illustrated in Fig. 4 (left) for a single direction and reference point.

Numerical algorithms
Many of the numerical algorithms for solving multiobjective optimization problems make use of scalarization with varying parameters. It is then possible to use single-objective numerical optimization methods for finding different points on the Pareto front.
Besides these, there are methods that focus on solving the Karush-Kuhn-Tucker conditions. These methods aim for covering all solutions to the typically underdetermined nonlinear equation system given by these condition. Again, for the sake of clarity and brevity, in the following treatment, we will focus on the unconstrained case, noting that the full Karush-Kuhn-Tucker and Fritz-John conditions also feature equality and inequality constraints (Kuhn and Tucker 1951).
Theorem 1 Fritz-John Conditions. A neccessary condition for x 2 X to be locally efficient is given by Theorem 2  Remark 14 The equation in the Fritz-John Condition typically does not result in a unique solution. For an ndimensional decision space X we have n þ 1 equations and we have m þ n unknowns (including the k multipliers). Hence, in a non-degenerate case, the solution set is of dimension m À 1.
It is possible to use continuation and homotopy methods to obtain all the solutions. The main idea of continuation methods is to find a single solution of the equation system and then to expand the solution set in the neighborhood of this solution. To decide in which direction to expand, it is necessary to maintain an archive, say A, of points that have already been obtained. To obtain a new point x new in the neighborhood of a given point from the archive x 2 A the homotopy method conducts the following steps: 1. Using the implicit function theorem a tangent space at the current point is obtained. It yielded an m À 1 dimensional hyperplane that is tangential to fðxÞ and used to obtain a predictor. See for the implicit function theorem, for instance, Krantz and Parks (2003 See Hillermeier (2001) and Schütze et al. (2005) for examples and more detailed descriptions. The continuation and homotopy methods require the efficient set to be connected. Moreover, they require points to satisfy certain regularity conditions (local convexity and differentiability). Global multiobjective optimization research is still a very active field of research. There are some promising directions, such as subdivision techniques , Bayesian global optimization , and Lipschitz optimization (Ž ilinskas 2013). However, these require the decision space to be of low dimension.
Moreover, there is active research on derivative-free methods for numerical multiobjective optimization. Direct search techniques have been devised, for instance, Custódio et al. (2011), and by Audet et al. (2010). For a summary of derivative-free methods, see Custódio et al. (2012).

Evolutionary multiobjective optimization
Evolutionary algorithms are a major branch of bio-inspired search heuristics, which originated in the 1960ties and are widely applied to solve combinatorial and non-convex numerical optimization problems. In short, they use paradigms from natural evolution, such as selection, recombination, and mutation to steer a population (set) of individuals (decision vectors) towards optimal or near-optimal solutions (Bäck 1996).
Multiobjective evolutionary algorithms (MOEAs) generalize this idea, and typically they are designed to gradually approach sets of Pareto optimal solutions that are well-distributed across the Pareto front. As there are-in general-no single-best solutions in multiobjective optimization, the selection schemes of such algorithms differ from those used in single-objective optimization. First MOEAs were developed in the 1990ties-see, e.g., Kursawe (1990) and Fonseca and Fleming (1993), but since around the year 2001, after the first book devoted exclusively to this topic was published by Deb (2001), the number of methods and results in this field grew rapidly.
With some exceptions, the distinction between different classes of evolutionary multiobjective optimization algorithms is mainly due to the differences in the paradigms used to define the selection operators, whereas the choice of the variation operators is generic and dependent on the problem. As an example, one might consider NSGA-II (see Deb et al. 2002) as a typical evolutionary multiobjective optimization algorithm; NSGA-II can be applied to continuous search spaces as well as to combinatorial search spaces. Whereas the selection operators stay the same, the variation operators (mutation, recombination) must be adapted to the representations of solutions in the decision space.
There are currently three main paradigms for MOEA designs. These are: 1. Pareto based MOEAs: The Pareto based MOEAs use a two-level ranking scheme. The Pareto dominance relation governs the first ranking and contributions of points to diversity is the principle of the second level ranking. The second level ranking applies to points that share the same position in the first ranking. NSGA-II (see Deb et al. 2002) and SPEA2 (see Zitzler and Thiele 1999) are two popular algorithms that fall into this category. In this tutorial, we will introduce typical algorithms for each of these paradigms: NSGA-II, SMS-EMOA, and Algorithm 1 NSGA-II Algorithm end for 11: {End variate} 12: {Selection step, select μ-"best" out of (Pt ∪ Qt) by a two step procedure:} 13: (R 1 , ..., R ) ← non-dom sort(f , Pt ∪ Qt) 14: Find the element of the partition, R iµ , for which the sum of the cardinalities |R 1 | + {End of selection step.} 16: t ← t + 1 17: end while 18: return Pt MOEA/D. We will discuss important design choices, and how and why other, similar algorithms deviate in these choices.

Pareto based algorithms: NSGA-II
The basic loop of NSGA-II (Deb et al. 2002) is given by Algorithm 1.
Firstly, a population of points is initialized. Then the following generational loop is repeated. This loop consists of two parts. In the first, the population undergoes a variation. In the second part, a selection takes place which results in the new generation-population. The generational loop repeats until it meets some termination criterion, which could be convergence detection criterion (cf. Wagner et al. 2009) or the exceedance of a maximal computational budget.
In the variation part of the loop k offspring are generated. For each offspring, two parents are selected. Each one of them is selected using binary tournament selection, that is drawing randomly two individuals from P t and selecting the better one concerning its rank in the population. The parents are then recombined using a standard recombination operator. For real-valued problems simulated binary crossover (SBX) is used (see Deb and Argawal 1995). Then the resulting individual is mutated. For real-valued problem polynomial mutation (PM) is used (see Mateo and Alberto 2012). This way, k individuals are created, which are all combinations or modifications of individuals in P t . Then the parent and the offspring populations are merged into P t [ Q t .
In the second part, the selection part, the l best individuals of P t [ Q t with respect to a multiobjective ranking are selected as the new population P tþ1 .
Next we digress in order to explain the multiobjective ranking which is used in NSGA-II. The key ingredient of NSGA-II that distinguishes it from genetic algorithms for single-objective optimization, is the way the individuals are ranked. The ranking procedure of NSGA-II consists of two levels. First, non-dominated sorting is performed. This ranking solely depends on the Pareto order and does not depend on diversity. Secondly, individuals which share the same rank after the first ranking are then ranked according to the crowding distance criterion which is a strong reflection of the diversity.
Let NDðPÞ denote the non-dominated solutions in some population. Non-dominated sorting partitions the populations into subsets (layers) based on Pareto non-dominance and it can be specified through recursion as follows.
As in each step of the recursion at least one solution is removed from the population, the maximal number of layers is |P|. We will use the index ' to denote the highest non-empty layer. The rank of the solution after non-dominated sorting is given by the subindex k of R k . It is clear that solutions in the same layer are mutually incomparable. The non-dominated sorting procedure is illustrated in Fig. 5 (upper left). The solutions are ranked as follows R 1 ¼ fy ð1Þ ; y ð2Þ ; y ð3Þ ; y ð4Þ g, R 2 ¼ fy ð5Þ ; y ð6Þ ; y ð7Þ g, R 3 ¼ fy ð8Þ ; y ð9Þ g. Now, if there is more than one solution in a layer, say R, a secondary ranking procedure is used to rank solutions within that layer. This procedure applies the crowding distance criterion. The crowding distance of a solution x 2 R is computed by a sum over contributions c i of the i-th objective function: The crowding distance is now given as: For m ¼ 2 the crowding distances of a set of mutually non-dominated points are illustrated in Fig. 5 (upper right). In this particular case, they are proportional to the perimeter of a rectangle that just is intersecting with the neighboring points (up to a factor of 1 4 ). Practically speaking, the value of l i is determined by the nearest neighbor of x to the left according to the i-coordinate, and l i is equal to the i-th coordinate of this nearest neighbor, similarly the value of u i is determined by the nearest neighbor of x to the right according to the i-coordinate, and u i is equal to the i-th coordinate of this right nearest neighbor. The more space there is around a solution, the higher is the crowding distance. Therefore, solutions with a high crowding distance should be ranked better than those with a low crowding distance in order to maintain diversity in the population. This way we establish a second order ranking. If the crowding distance is the same for two points, then it is randomly decided which point is ranked higher.
Now we explain the non-dom_sort procedure in line 13 of Algorithm 1 the role of P is taken over by P t \ Q t : In order to select the l best members of P t [ Q t according to the above described two level ranking, we proceed as follows. Create the partition R 1 ; R 2 ; Á Á Á ; R ' of P t [ Q t as described above. For this partition one finds the first index i l for which the sum of the cardinalities jR 1 j þ Á Á Á þ jR i l j is for the first time ! l. If jR 1 j þ Á Á Á þ jR i l j ¼ l, then set P tþ1 to [ i l i¼1 R i , otherwise determine the set H containing l À ðjR 1 j þ Á Á Á þ jR i l À1 jÞ elements from R i l with the highest crowding distance and set the next generationpopulation, P tþ1 , to ð[ Pareto-based Algorithms are probably the largest class of MOEAs. They have in common that they combine a ranking criterion based on Pareto dominance with a diversity based secondary ranking. Other common algorithms that belong to this class are as follows. The Multiobjective Genetic Algorithm (MOGA) (Fonseca and Fleming 1993), which was one of the first MOEAs. The PAES (Knowles and Corne 2000), which uses a grid partitioning of the objective space in order to make sure that certain regions of the objective space do not get too crowded. Within a single grid cell, only one solution is selected. The Strength Pareto Evolutionary Algorithm (SPEA) (Zitzler and Thiele 1999) uses a different criterion for ranking based on Pareto dominance. The strength of an individual depends on how many other individuals it dominates and by how many other individuals dominate it. Moreover, clustering serves as a secondary ranking criterion. Both operators have been refined in SPEA2 (Zitzler et al. 2001), and also it features a strategy to maintain an archive of non-dominated solutions. The Multiobjective Micro GA. Coello and Pulido (2001) is an algorithm that uses a very small population size in conjunction with an archive. Finally, the Differential Evolution Multiobjective Optimization (DEMO) (Robic and Filipic 2005) algorithm combines concepts from Pareto-based MOEAs with a variation operator from differential evolution, which leads to improved efficiency and more precise results in particular for continuous problems.

Indicator-based algorithms: SMS-EMOA
A second algorithm that we will discuss is a classical algorithm following the paradigm of indicator-based multiobjective optimization. In the context of MOEAs, by a performance indicator (or just indicator), we denote a scalar measure of the quality of a Pareto front approximation. Indicators can be unary, meaning that they yield an absolute measure of the quality of a Pareto front approximation. They are called binary, whenever they measure how much better one Pareto front approximation is concerning another Pareto front approximation.
The SMS-EMOA ) uses the hypervolume indicator as a performance indicator. Theoretical analysis attests that this indicator has some favorable properties, as the maximization of it yields approximations of the Pareto front with points located on the Pareto front and well distributed across the Pareto front. The hypervolume indicator measures the size of the dominated space, bound from above by a reference point.
For an approximation set A & R m it is defined as follows: HIðAÞ ¼ Volðfy 2 R m : y " Pareto r^9a 2 A : a " Pareto ygÞ Here, Vol ð:Þ denotes the Lebesgue measure of a set in dimension m. This is length for m ¼ 1, area for m ¼ 2, volume for m ¼ 3, and hypervolume for m ! 4. Practically speaking, the hypervolume indicator of A measures the size of the space that is dominated by A. The closer points move to the Pareto front, and the more they distribute along the Pareto front, the more space gets dominated. As the size of the dominated space is infinite, it is necessary to bound it. For this reason, the reference point r is introduced.
The SMS-EMOA seeks to maximize the hypervolume indicator of a population which serves as an approximation set. This is achieved by considering the contribution of points to the hypervolume indicator in the selection procedure. Algorithm 2 describes the basic loop of the standard implementation of the SMS-EMOA.
The algorithm starts with the initialization of a population in the search space. Then it creates only one offspring individual by recombination and mutation. This new individual enters the population, which has now size l þ 1. To reduce the population size again to the size of l, a subset of size l with maximal hypervolume is selected. This way as long as the reference point for computing the hypervolume remains unchanged, the hypervolume indicator of P t can only grow or stay equal with an increasing number of iterations t.
Next, the details of the selection procedure will be discussed. If all solutions in P t are non-dominated, the selection of a subset of maximal hypervolume is equivalent to deleting the point with the smallest (exclusive) hypervolume contribution. The hypervolume contribution is defined as: DHIðy; YÞ ¼ HIðYÞ À HIðY n fygÞ An illustration of the hypervolume indicator and hypervolume contributions for m ¼ 2 and, respectively, m ¼ 3 is given in Fig. 6. Efficient computation of all hypervolume contributions of a population can be achieved in time Hðl log lÞ for m ¼ 2 and m ¼ 3 (Emmerich and Fonseca 2011). For m ¼ 3 or 4, fast implementations are described in Guerreiro and Fonseca (2017). Moreover, for fast logarithmic-time incremental updates for 2-D an algorithm is described in Hupkens and Emmerich (2013). For achieving logarithmic time updates in SMS-EMOA, the non-dominated sorting procedure was replaced by a procedure, that sorts dominated solutions based on age. For m [ 2, fast incremental updates of the hypervolume indicator and its contributions were proposed in for more than two dimensions (Guerreiro and Fonseca 2017).
In case dominated solutions appear the standard implementation of SMS-EMOA partitions the population into layers of equal dominance ranks, just like in NSGA-II. Subsequently, the solution with the smallest hypervolume contribution on the worst ranked layer gets discarded.
SMS-EMOA typically converges to regularly spaced Pareto front approximations. The density of these approximations depends on the local curvature of the Pareto front. For biobjective problems, it is highest at points where the slope is equal to À45 (Auger et al. 2009). It is possible to influence the distribution of the points in the approximation set by using a generalized cone-based hypervolume indicator. These indicators measure the hypervolume dominated by a cone-order of a given cone, and the resulting optimal distribution gets more uniform if the cones are acute, and more concentrated when using obtuse cones (see Emmerich et al. 2013).
Besides the SMS-EMOA, there are various other indicator-based MOEAs. Some of them also use the hypervolume indicator. The original idea to use the hypervolume indicator in an MOEA was proposed in the context of archiving methods for non-dominated points. Here the hypervolume indicator was used for keeping a boundedsize archive (Knowles et al. 2003). Besides, in an early work hypervolume-based selection which also introduced a novel mutation scheme, which was the focus of the paper (Huband et al. 2003). The term Indicator-based Evolutionary Algorithms (IBEA) (Zitzler and Künzli 2004) was introduced in a paper that proposed an algorithm design, in which the choice of indicators is generic. The hypervolume-based IBEA was discussed as one instance of this class. Its design is however different to SMS-EMOA and makes no specific use of the characteristics of the hypervolume indicator. The Hypervolume Estimation Algorithm (HypE) (Bader and Zitzler 2011) uses a Monte Carlo Estimation for the hypervolume in high dimensions and thus it can be used for optimization with a high number of objectives (so-called many-objective optimization problems). MO-CMA-ES (Igel et al. 2006) is another hypervolume-based MOEA. It uses the covariance-matrix adaptation in its mutation operator, which enables it to adapt its mutation distribution to the local curvature and scaling of the objective functions. Although the hypervolume indicator has been very prominent in IBEAs, there are some algorithms using other indicators, notably this is the R2 indicator (Trautmann et al. 2013), which features an ideal point as a reference point, and the averaged Hausdorff distance (D p indicator) (Rudolph et al. 2016), which requires an aspiration set or estimation of the Pareto front which is dynamically updated and used as a reference. The idea of aspiration sets for indicators that require knowledge of the 'true' Pareto front also occurred in conjunction with the a-indicator (Wagner et al. 2015), which generalizes the approximation ratio in numerical single-objective optimization. The Portfolio Selection Multiobjective Optimization Algorithm (POSEA) (Yevseyeva et al. 2014) uses the Sharpe Index from financial portfolio theory as an indicator, which applies the hypervolume indicator of singletons as a utility function and a definition of the covariances based on their overlap. The Sharpe index combines the cumulated performance of single individuals with the covariance information (related to diversity), and it has interesting theoretical properties.

Decomposition-based algorithm: MOEA/D
Decomposition-based algorithms divide the problem into subproblems using scalarizations based on different weights. Each scalarization defines a subproblem. The subproblems are then solved simultaneously by dynamically assigning and re-assigning points to subproblems and exchanging information from solutions to neighboring subproblems.
The method defines neighborhoods on the set of these subproblems based on the distances between their aggregation coefficient vectors. When optimizing a subproblem, information from neighboring subproblems can be exchanged, thereby increasing the efficiency of the search as compared to methods that independently solve the subproblems.
MOEA/D (Zhang and Li 2007) is a very commonly used decomposition based method, that succeeded a couple of preceding algorithms based on the idea of combining decomposition, scalarization and local search (Ishibuchi and Murata 1996;Jin et al. 2001;Jaszkiewicz 2002). Note that even the early proposed algorithms VEGA (Schaffer 1985) and the vector optimization approach of Kursawe (see Kursawe 1990) can be considered as rudimentary decomposition based approaches, where these algorithms obtain a problem decomposition by assigning different members of a population to different objective functions. These early algorithmic designs used subpopulations to solve different scalarized problems. In contrast, in MOEA/ D one population with interacting neighboring individuals is applied, which reduces the complexity of the algorithm.
Typically, MOEA/D works with Chebychev scalarizations, but the authors also suggest other scalarization methods, namely scalarization based on linear weightingwhich however has problems with approximating nonconvex Pareto fronts-and scalarization based on boundary intersection methods-which requires additional parameters and might also obtain strictly dominated points.
MOEA/D evolves a population of individuals, each individual x ðiÞ 2 P t being associated with a weight vector k ðiÞ . The weight vectors k ðiÞ are evenly distributed in the search space, e.g., for two objectives a possible choice is: . . .; l. The i-th subproblem gðxjk i ; z Ã Þ is defined by the Chebychev scalarization function (see also Eq. 2): The main idea is that in the creation of a new candidate solution for the i-th individual the neighbors of this individual are considered. A neighbor is an incumbent solution of a subproblem with similar weight vectors. The neighborhood of i-th individual is the set of k subproblems, for so predefined constant k, that is closest to k ðiÞ in the Euclidean distance, including the i-th subproblem itself. It is denoted with B(i). Given these preliminaries, the MOEA/ D algorithm-using Chebychev scalarization-reads as described in Algorithm 3. Note the following two remarks about MOEA/D: (1) Many parts of the algorithm are kept generic. Here, generic options are recombination, typically instantiated by standard recombination operators from genetic algorithms, and local improvement heuristic. The local improvement heuristic should find a solution in the vicinity of a given solution that does not violate constraints and has a relatively good performance concerning the objective function values.
(2) MOEA/D has additional statements to collect all non-dominated solutions it generates during a run in an external archive. Because this external archive is only used in the final output and does not influence the search dynamics, it can be seen as a generic feature of the algorithm. In principle, an external archive can be used in all EMOAs and could therefore also be done in SMS-EMOA and NSGA-II. To make comparisons to NSGA-II and SMS-EMOA easier, we omitted the archiving strategy in the description.
Recently, decomposition-based MOEAs became very popular, also because they scale well to problems with many objective functions. The NSGA-III (Deb and Jain 2014) algorithm is specially designed for many-objective optimization and uses a set of reference points that is dynamically updated in its decomposition. Another decomposition based technique is called Generalized Decomposition (Giagkiozis et al. 2014). It uses a mathematical programming solver to compute updates, and it was shown to perform well on continuous problems. The combination of mathematical programming and decomposition techniques is also explored in other, more novel, hybrid techniques, such as Directed Search (Schütze et al.

Algorithm 3 MOEA/D
input: Λ = {λ (1) , ..., λ (µ) } {weight vectors} input: z * : reference point for Chebychev distance initialize P 0 ⊂ X µ initialize neighborhoods B(i) by collecting k nearest weight vectors in Λ for each λ (i) while not terminate do for all i ∈ {1, ..., μ} do Select randomly two solutions x (1) , x (2) in the neighborhood B(i). y ← Recombine x (1) , x (2) by a problem specific recombination operator. y ← Local problem specific, heuristic improvement of y, e.g. local search, based on the scalarized objective function g(x|λ (i) , z * ) . if g(y |λ (i) , z * ) < g(x (i) |λ (i) , z * ) then x (i) ← y end if Update z * , if neccessary, i.e, one of its component is larger than one of the corresponding components of f (x (i) ). end for t ← t + 1 end while return Pt 2016), which utilizes the Jacobian matrix of the vectorvalued objective function (or approximations to it) to find promising directions in the search space, based on desired directions in the objective space.

Performance assessment
In order to make a judgement (that is, gain insight into the advantages and disadvantages) of multiobjective evolutionary (or for that matter also deterministic) optimizers we need to take into account (1) the computational resources used, and (2) the quality of the produced result(s).
The current state of the art of multiobjective optimization approaches are mainly compared empirically though theoretical analyses exist (see, for instance, the convergence results described in Rudolph and Agapie (2000), Beume et al. (2011) albeit for rather simple problems as more realistic problems elude mathematical analysis.
The most commonly computational resource which is taken into account is the computation time which is very often measured implicitly by counting fitness function evaluations-in this respect, there is no difference with single-objective optimization. In contrast to single-objective optimization, in multiobjective optimization, a close distance to a (Pareto) optimal solution is not the only thing required but also good coverage of the entire Pareto front. As the results of multiobjective optimization algorithms are (finite) approximation sets to the Pareto front we need to be able to say when one Pareto front approximation is better than another. One good way to define when one approximation set is better than another is as in Definition 22 (see Zitzler et al. 2003).
Definition 21 Approximation Set of a Pareto Front. A finite subset A of R m is an approximation set of a Pareto front if and only if A consists of mutually (Pareto) nondominated points.
Definition 22 Comparing Approximation Sets of a Pareto Front. Let A and B be approximation sets of a Pareto front in R m . We say that A is better than B if and only if every b 2 B is weakly dominated by at least one element a 2 A and A 6 ¼ B. Notation: A.B.
In Fig. 7 examples are given of the case of one set being better than another and in Fig. 8 examples are given of the case that a set is not better than another.
This relation on sets has been used in Zitzler et al. (2003) to classify performance indicators for Pareto fronts. To do so, they introduced the notion of completeness and compatibility of these indicators with respect to the set relation 'is better than'.
Definition 23 Unary Set Indicator. A unary set indicator is a mapping from finite subsets of the objective space to the set of real numbers. It is used to compare (finite) approximations to the Pareto front.
Definition 24 Compatibility of Unary Set Indicators concerning the 'is better than' order on Approximation Sets. A unary set indicator I is compatible concerning the 'is better than' or . The hypervolume indicator and some of its variations are complete. Other indicators compared in the paper ) are weakly-complete or not even weakly-complete. It has been proven in the same paper that no unary indicator exists that is complete and compatible at the same time. Moreover for the hypervolume indicator HI it has be shown that HI ðAÞ [ HI ðBÞ ) :ðB.AÞ. The latter we call weakly-compatible.
In all the discussions of the hypervolume indicator, the assumption is that all points of the approximation sets under consideration dominate the reference point. Recently, variations of the hypervolume indicator have been proposed-the so-called free hypervolume indicators-that do not require the definition of a reference point  Fig. 7 In the left picture, the set of points denoted by blue squares is better than (.) the set consisting of the red-circle points. Also in the right picture the set consisting of blue squares is better than the set of red-circle points-in this case the intersection of the two sets is non-empty and are complete and weakly-compatible for all approximation sets .
Besides unary indicators, one has introduced binary indicators (see Riquelme et al. 2015). The most used ones are unary indicators followed up by binary indicators. For binary indicators, the input is a pair of approximation sets and the output is again a real number. Here the goal is to determine which of the two approximation sets is the better one (and how much better it is) 1 . Binary indicators can also be used, if the true Pareto front is known, e.g., in benchmarking on test problems. A common binary indicator is the binary -indicator. In order to introduce this indicator we first define for each d 2 R a binary relation on the points in R m .
Definition 25 d-domination. Let d 2 R and let a 2 R m and b 2 R m . We say that a d-dominates b (notation: Next, we can define the binary indicator I .

Definition 26
The Binary Indicator I . Given two approximation sets A and B, then I ðA; BÞ :¼ inf d2R f8b 2 B 9a 2 A such that a " d bg.
Clearly for a fixed B the smaller I ðA; BÞ is the better the approximation set A is relative to B. The following properties hold: A.B ) I ðB; AÞ [ 0, the second notable property is as follows: I ðA; BÞ 0 and I ðB; AÞ [ 0 ) A.B. These two properties show that based on the binary -indicator it is possible to decide whether or not A is better than B. In contrast, the knowledge of the hypervolume indicator on the sets A and B does not allow to decide whether or not A is better than B.
Some of indicators are useful in case there is knowledge or complete knowledge about the Pareto front. For instance (see Rudolph et al. 2016), it has been suggested to compute the Hausdorff distance (or variations of it) of an approximation set to the Pareto front. Moreover, the binaryindicator can be transformed into a complete unary indicator in case the second input is the known Pareto frontnote that this indicator needs to be minimized.
Another useful way to learn about the behavior of evolutionary multiobjective algorithms is the attainment curve approach (see da Fonseca et al. 2001). The idea is to generalize the cumulative distribution function and for the study of algorithms it is approximated empirically. The distribution is defined on the set of (finite) approximation sets of the Pareto front. For each point in the objective space R m it is the probability that the Pareto front approximation attains this point (that is, it is either one point in the approximation set, or it is dominated by some point in the approximation set). Formally, it reads  Fig. 8 In each of the pictures, the set consisting of the blue square points is not better than the set consisting of the red circle points. Clearly, in each of the two pictures on the right the set consisting of the red circle points is better than the set consisting of the blue square points Pða ð1Þ z _ a ð2Þ z _ . . . _ a ðkÞ zÞ; where A ¼ fa ð1Þ ; a ð2Þ ; . . .; a ðkÞ g is the approximation set and z 2 R m . In other words P is the probability of an algorithm to find at least one solution which attains the goal z in a single run. We define the attainment function a A : R m ! ½0; 1 associated to the approximation set A as follows: where A 1 ; . . .; A s are the outcome approximation sets of an algorithm in s runs of the algorithm and I : ð set of approximation sets Þ Â R m ! f0; 1g is a function which associates to an approximation set and a vector in R m the value 1 in case the vector is a member of the approximation set or if some element of the approximation set dominates it, otherwise the value is 0. For m ¼ 2 or 3 we can plot the boundaries where this function changes its value. These are the attainment curves (m ¼ 2Þ and attainment surfaces (m ¼ 3). In particular the median attainment curve/surface gives a good summary of the behavior of the optimizer. It is the boundary where the function changes from a level below 0.5 to a level higher than 0.5. Alternatively one can look at lower and higher levels than 0.5 in order to get an optimistic or respectively a pessimistic assessment of the performance. In Fig. 9 an example of the median attainment curve is shown. We assume that the four approximation sets are provided by some algorithm.

Recent topics in multiobjective optimization
Recently, there are many new developments in the field of multiobjective optimization. Next we will list some of the most important trends.

Many-objective optimization
Optimization with more than 3 objectives is currently termed many-objective optimization [see, for instance, the survey ]. This is to stress the challenges one meets when dealing with more than 3 objectives. The main reasons are: 1. problems with many objectives have a Pareto front which cannot be visualized in conventional 2D or 3D plots instead other approaches to deal with this are needed; 2. the computation time for many indicators and selection schemes become computationally hard, for instance, time complexity of the hypervolume indicator computation grows super-polynomially with the number of objectives, under the assumption that P 6 ¼ NP; 3. last but not least the ratio of non-dominated points tends to increase rapidly with the number of objectives. For instance, the probability that a point is nondominated in a uniformly distributed set of sample points grows exponentially fast towards 1 with the number of objectives.
In the field of many-objective optimization different techniques are used to deal with these challenges. For the first challenge, various visualization techniques are used such as projection to lower dimensional spaces or parallel coordinate diagrams. In practice, one can, if the dimension is only slightly bigger than 3, express the coordinate values by colors and shape in 3D plots. Naturally, in many-objective optimization indicators which scale well with the number of objectives (say polynomially) are very much desired. Moreover, decomposition based approaches are typically preferred to indicator based approaches.
The last problem requires, however, more radical deviations from standard approaches. In many cases, the reduction of the search space achieved by reducing it to the efficient set is not sufficiently adequate to allow for  Fig. 9 The median attainment curve for the case of four approximation sets; one approximation set consists of the blue squares, the second set consists of points denoted by brown triangles, the third consists of the red circles, and the fourth consists of points denoted by black crosses; the darker gray the region is the more approximation sets dominate it. The median attainment curve is the black polygonal line subsequent decision making because too many alternatives remain. In such cases, a stricter order than the Pareto order is required which requires additional preference knowledge. To elicit preference knowledge, interactive methods often come to the rescue. Moreover, in some cases, objectives are correlated which allows for grouping of objectives, and in turn, such groups can be aggregated to a single objective. Dimensionality reduction and community detection techniques have been proposed for identifying meaningful aggregation of objective functions.

Preference modeling
The Pareto order is the most applied order in multiobjective optimization. However, different ranking schemes and partial orders have been proposed in the literature for various reasons. Often additional knowledge of user preferences is integrated. For instance, One distinguishes preference modeling according to at what stage of the optimization the preference information is collected (a priori, interactively, and a posteriori). Secondly one can distinguish the type of information requested from the decision maker, for instance, constraints on the trade-offs, relative importance of the objectives, and preferred regions on the Pareto front. Another way to elicit preference information is by ordinal regression; here the user is asked for pairwise comparisons of some of the solutions. From this data, the weights of utility functions are learned (Branke et al. 2015). The interested reader is also referred to interesting work on non-monotonic utility functions, using the Choquet integral (Branke et al. 2016). Notably, the topic of preference elicitation is one of the main topics in Multiple Criteria Decision Analysis (MCDA). In recent years collaboration between MCDA and multiobjective optimization (MOO) brought forward many new useful approaches. A recommended reference for MCDA is Belton and Stewart (2002). For a comprehensive overview of preference modelling in multiobjective optimization we refer to Li et al. (2017) and Hakanen et al. (2016). Moreover Greco et al. (2016) contains an updated collection of state of the art surveys on MCDA. A good reference discussing the integration of MCDA into MOO is Branke et al. (2008).

Optimization with costly function evaluations
In industrial optimization very often the evaluation of (an) objective function(s) is achieved by simulation or experiments. Such evaluations are typically time-consuming and expensive. Examples of such costly evaluations occur in the optimization based on crash tests of automobiles, chemical experiments, computational fluid dynamics simulations, and finite element computations of mechanical structures. To deal with such problems techniques that need only a limited number of function evaluations have been devised. A common approach is to learn a surrogate model of the objective functions by using all available past evaluations. This is called surrogate model assisted optimization. One common approach is to optimize on the surrogate model to predict promising locations for evaluation and use these evaluations to further improve the model. In such methods, it is also important to add points for developing the model in under-explored regions of the search space. Some criteria such as expected improvement take both exploitation and exploration into account. Secondly, surrogate models can be used in pre-processing in the selection phase of evolutionary algorithms. To save time, less interesting points can be discarded before they would be evaluated by the costly and precise evaluator. Typically regression methods are used to construct surrogate models; Gaussian processes and neural networks are standard choices. Surrogate modeling has in the last decade been generalized to multiobjective optimization in various ways. Some important early work in this field was on surrogate assisted MOEAs (Emmerich et al. 2006) and ParEGO algorithm (Knowles 2006). A state of the art review can be found in Allmendinger et al. (2017).

New bio-inspired paradigms
Inspiration by nature has been a creative force for dealing with optimization algorithm design. Apart from biological evolution, many other natural phenomena have been considered. While many of these algorithmic ideas have so far remained in a somewhat experimental and immature state, some non-evolutionary bio-inspired optimization algorithms have gained maturation and competitive performance. Among others, this seems to hold for particle swarm optimization (Reyes-Sierra and Coello Coello 2006), ant colony optimization (Barán and Schaerer 2003), and artificial immune systems Coello Coello and Cortés (2005). As with evolutionary algorithms, also these algorithms have first been developed for single-objective optimization, and subsequently, they have been generalized to multiobjective optimization. Moreover, there is some recent research on bio-inspired techniques that are specifically developed for multiobjective optimization. An example of such a development is the Predator-Prey Evolutionary Algorithm, where different objectives are represented by different types of predators to which the prey (solutions) have to adapt (Laumanns et al. 1998;Grimme and Schmitt 2006). In the field of natural computing, it is also investigated whether algorithms can serve as models for nature. It is an interesting new research direction to view aspects of natural evolution as a multiobjective optimization process, and first such models have been explored in Rueffler (2006) and Sterck et al. (2011).

Set-oriented numerical optimization
Traditionally, numerical techniques for multiobjective optimization are single point techniques: They construct a Pareto front by formulating a series of single-objective optimization problems (with different weights or constraints) or by expanding a Pareto front from a single solution point by point using continuation. In contrast, setoriented numerical multiobjective optimization operates on the level of solution sets, the points of which are simultaneously improved, and that converge to a well-distributed set on the Pareto front. Only very recently such methods have been developed, and techniques that originated from evolutionary multiobjective optimization have been transferred into deterministic methods. A notable example is the hypervolume indicator gradient ascent method for multiobjective optimization (HIGA-MO) (Wang et al. 2017). In this method a set of say l points is represented as a single vector of dimension ln, where n is the dimension of the search space. In real-valued decision space the mapping HI: R ld ! R from the such population vectors to the hypervolume indicator has a well-defined derivative in almost all points. The computation of such derivatives has been described in Emmerich and Deutz (2014). Viewing the vector of partial derivatives as a gradient, conventional gradient methods can be used. It requires, however, some specific adaptations in order to construct robust and practically usable methods for local optimization. On convex problems, fast linear convergence can be achieved. By using second-order derivatives in a hypervolume-based Newton-Raphson method, even quadratic convergence speed has been demonstrated empirically on a set of convex bi-objective problems. The theory of such secondorder methods is subject to ongoing research (Hernández et al. 2014).

Advanced performance assessment
Despite significant progress on the topic of performance assessment in recent years, there are still many unanswered questions. A bigger field of research is on performance assessment of interactive and many objective optimization. Moreover, the dependency of performance measures on parameters, such as the reference point of the hypervolume indicator requires further investigation. Some promising work in that direction was recently provided in Ishibuchi et al. (2017).
9 How to get started?
In the following, we list some of the main resources for the field of (Evolutionary) Multiobjective Optimization. Aside from the resources mentioned above, there are many research groups and labs which maintain a repository of software accompanying their published research, e.g., the MODA group at Leiden University http://moda.liacs.nl and the research group of Carlos Fonseca at Coimbra University eden.dei.uc.pt/cmfonsec/software.html.

Summary and outlook
In this tutorial, we gave an introduction to the field of multiobjective optimization. We covered the topics of order-theoretical foundations, scalarization approaches, and optimality conditions. As solution methods, we discussed homotopy and evolutionary methods. In the context of Evolutionary methods, we discussed three state-of-theart techniques in detail, namely NSGA-II, SMS-EMOA, and MOEA/D, each representing a key paradigm in evolutionary multiobjective algorithm design. NSGA-II served as a representative of Pareto based approaches, SMS-EMOA as an example of Indicator-based approaches, and MOEA/D as an example of decomposition based approaches. These algorithms have some advantages and disadvantages: -Pareto-based approaches follow a straightforward design principle, that is directly based on Pareto dominance and diversity preservation (using, for instance, crowding distance). Usually, these algorithms require only a few parameters, and larger numbers of objective functions do not cause problems. However, it might be difficult to guarantee and measure convergence and achieve a very regular spacing of solutions. -Indicator-based approaches use an indicator for the performance of an approximation set to guide the search. It is possible to assess their convergence behavior online, and they hold the promise to be more amenable to theoretical analysis. However, the computation time often increases rapidly with the number of dimensions and the distribution of points in the approximation sets might depend critically on the settings of the reference point or other parameters of the indicator. -Decomposition-based methods provide a very flexible framework for algorithm design, as they can incorporate various scalarization methods. A disadvantage is that they require some a priori knowledge of the position of the Pareto front in the objective space and the number of weight vectors might grow exponentially with the objective space size, even if the Pareto front is of low dimension.
According to the above, choosing the right method depends much on the dimension of the objective space, the number of solutions one wants to output, the desired distribution of the solutions (knee-point focused or uniformly spread) and the a priori knowledge on the location and shape of the Pareto front. Due to space constraints, many advanced topics in multiobjective optimization are not covered in depth. We refer for these topics to the literature. For instance, constraint handling (Coello Coello 2013), multimodality (Kerschke et al. 2016), non-convex global optimization (Ž ilinskas 2013), and combinatorial optimization (Ehrgott and Gandibleux 2000).
Multiobjective Optimization is a very active field of research. There are still many open, challenging problems in the area. For future development of the research field it will be essential to provide EMO algorithms that are built around a robust notion of performance and, ideally, also can be analyzed more rigorously. Major topics for current research are also uncertainty handling and robustness, many-objective optimization, theoretical foundations and computational complexity, generalizations, for instance, level set approximation, diversity optimization, and setoriented optimization, customization and integration into multidisciplinary workflows, and scalability to big data, or expensive evaluations.