Abstract The use of surrogate based optimization (SBO) is widely spread in engineering design to reduce the number of computational expensive simulations. However, "real-world" problems often consist of multiple, conflicting objectives leading to a set of competitive solutions (the Pareto front). The objectives are often aggregated into a single cost function to reduce the computational cost, though a better approach is to use multiobjective optimization methods to directly identify a set of Pareto-optimal solutions, which can be used by the designer to make more efficient design decisions (instead of weighting and aggregating the costs upfront). Most of the work in multiobjective optimization is focused on multiobjective evolutionary algorithms (MOEAs). While MOEAs are well-suited to handle large, intractable design spaces, they typically require thousands of expensive simulations, which is prohibitively expensive for the problems under study. Therefore, the use of surrogate models in multiobjective optimization, denoted as multiobjective surrogate-based optimization, may prove to be even more worthwhile than SBO methods to expedite the optimization of computational expensive systems. In this paper, the authors propose the efficient multiobjective optimization (EMO) algorithm which uses Kriging models and multiobjective versions of the probability of improvement and expected improvement criteria to identify the Pareto front with a minimal number of expensive simulations. The EMO algorithm is applied on multiple standard benchmark problems and compared against the well-known NSGA-II, SPEA2 and SMS-EMOA multiobjective optimization methods. Keywords Multiobjective optimization · Expected improvement · Probability of improvement · Hypervolume · Kriging · Gaussian process Surrogate modeling techniques, also known as metamodeling, are becoming rapidly popular in the engineering community to speed up complex, computational expensive design I. Couckuyt (B) · D. Deschrijver · T. Dhaene Department of Information Technology (INTEC), Ghent University-iMinds, Ghent, Belgium e-mail: ivo.couckuyt@ugent.be 123 576 J Glob Optim (2014) 60:575-594 problems [22, 37] . Surrogate models, or metamodels, are mathematical approximation models that mimic the behavior of computational expensive simulation codes such as mechanical or electrical finite element simulations, or computational fluid dynamic simulations. This paper deals with the use of surrogate models for expediting the optimization of time-consuming (black-box) problems of a deterministic nature, in contrast to stochastic simulation. While several types of surrogate modeling uses can be distinguished, this work is concerned with the integration of surrogate models into the optimization process, often denoted by surrogate based optimization (SBO) or metamodel-assisted optimization (MAO). SBO methods typically generate surrogate models on the fly that are only accurate in certain regions of the input space, e.g., around potentially optimal regions. The generated surrogate models can then be used to intelligently guide the optimization process to the global optimum. The focus of this work is the global SBO method based on the probability of improvement (PoI) and expected improvement (EI), popularized by Jones et al. [25] . These "statistical criteria" guide the selection of new data points in such a way that the objective function is optimized, while minimizing the number of expensive simulations. The advantage of EI and PoI is that, besides the prediction (mean), the uncertainty (variance) of the surrogate model is taken into account as well, providing a balance between exploration 1 and exploitation. 2 Most often EI or PoI is used in conjunction with the Kriging surrogate model (Gaussian Processes) [27] which provides by construction a prediction of the mean as well as the variance, but other surrogate models are also possible, such as radial basis functions (RBF), support vector regression (SVR) [13] , etc. The single-objective SBO problem is well described in literature, however, most (if not all) "real-world" problems actually consists of multiple, conflicting objectives leading to a set of Pareto-optimal solutions. Often the objectives are aggregated into a single cost function, e.g., using a weighted sum, that can be optimized by standard optimization techniques. Subsequently, by repeating this process many times using varying starting conditions, e.g., different set of weights, several solutions on the Pareto front can be found. On the other hand, a multiobjective optimization method can optimize the different objective functions simultaneously, and try to find the Pareto front in just a single run. Examples of such methods are primarily the multiobjective evolutionary algorithms (MOEAs), e.g., the "Non-dominated Sorting Genetic Algorithm II" (NSGA-II; [14]), the "Strength Pareto Evolutionary Algorithm 2" (SPEA2; [44] ) and the "S-Metric Selection Evolutionary MultiObjective Algorithm" (SMS-EMOA; [5] ). Unfortunately, MOEAs typically require a massive amount of function evaluations, which is infeasible for computational expensive simulators. Hence, it is vital to economize on the number of function evaluations, e.g., by using surrogate models. multiobjective surrogatebased optimization (MOSBO) methods only appeared quite recently in literature. Most work is focused on integrating surrogate models in MOEAs [41] . use neural networks to either approximate the fitness function or as a local approximation technique to generate search points more efficiently. Voutchkov and Keane [35] apply the NSGA-II algorithm to Kriging models instead of the expensive simulator. For an overview of available techniques and approaches the reader is referred to [30, 42] . While the PoI and EI approach is well-developed and used for single-objective SBO, its use in MOSBO is not well spread. Single-objective versions of EI and PoI are utilized by Knowles [28] and Knowles et al. [29] to solve MOSBO problems. This approach, known as ParEGO, uses Kriging and EI to optimize a weighted sum of objective functions. By 1 Improving the overall accuracy of the surrogate model (space-filling). 2 Enhancing the accuracy of the surrogate model solely in the region of the (current) optimum. J Glob Optim (2014) 60:575-594 577 randomizing the weights every iteration several solutions along the Pareto front can be identified. More recently, Keane [26] proposed multiobjective versions of PoI and Euclidean distance-based EI. At the same time Emmerich et al. [17] proposed the hypervolume-based EI criterion. Similarly to a weighted sum, the multiobjective versions of EI and PoI aggregate information from the surrogate models into a single cost function, balancing between exploration (see footnote 1) and exploitation. 3 Unfortunately, only formulae for two objective functions are given by Keane as the statistical criteria become rather cumbersome and complex for a higher number of objective functions. Similarly, while Emmerich et al. [16] describe formulae for an arbitrary number of dimensions for the hypervolume-based EI, the computation cost increases at least exponentially with the number of objectives and, hence, has only been applied to two objectives. The key contribution of this paper is the efficient multiobjective optimization (EMO) algorithm which is a much more efficient method of evaluating multiobjective versions of the PoI and EI criteria for multiobjective optimization problems. In fact, the problem at hand is similar to calculating the hypervolume (a Pareto set quality estimator) [45] as will be shown below and, hence, hypervolume algorithms can be adapted to aid in the evaluation of the statistical criteria. Moreover, a new statistical criterion is proposed, based on the hypervolume-based EI, which is significantly cheaper to compute while still delivering promising results. In Sect. 2 the Kriging surrogate model is briefly discussed. In Sect. 3, an overview of the EMO algorithm is given, including general expressions for the PoI and several variants of EI. Subsequently, a fundamental part needed for the calculation of the statistical criteria is discussed in Sect. 3.4. Afterwards, in Sect. 4 the EMO algorithm is tested on several functions from the DTLZ benchmark suite [15] . Lastly, in Sect. 5 conclusions and future work are discussed. Kriging is a popular surrogate model to approximate deterministic noise-free data, and has proven to be very useful for tasks such as optimization [25] , design space exploration, visualization, prototyping, and sensitivity analysis [37] . A thorough mathematically treatment of Kriging is given in [19, 33] . Basically, Kriging is a two-step process: first a regression function h(x) is constructed, and, subsequently, a centered Gaussian process Z with variance σ 2 and a correlation matrix Ψ is constructed through the residuals. Consider a set of n samples, ( 2) and associated function values, y = (y 1 , . . . , y n ) , where (·) is the transpose of a vector or matrix. Essentially, the regression part is encoded in the n× p model matrix F using basis functions  while the stochastic process is mainly defined by the n × n correlation matrix Ψ , where ψ(·, ·) is the correlation function. ψ(·, ·) is parameterized by a set of hyperparameters θ . The choice of correlation function is crucial to obtain good accuracy. This paper focuses on using the Matérn correlation function [34] , with ν = 3 2 , Matérn ν= In addition, the popular Gaussian correlation function is also used, The hyperparameters θ are identified by Maximum Likelihood Estimation (MLE). In particular, the negative concentrated log-likelihood is minimized, where . Subsequently, the prediction mean and prediction variance of Kriging are derived, respectively, as, where is the model matrix of the predicting point x, α is a p × 1 vector denoting the coefficients of the regression function, determined by generalized least squares (GLS), and r (x) is an 1 × n vector of correlations between the point x and the samples X . Efficient multiobjective optimization (EMO) A flow chart of the EMO algorithm is shown in Fig. 1 . First an initial set of points X is generated and evaluated on the expensive objective functions f j (x), for j = Kriging models are updated with this new information and this process is repeated in an iterative fashion until some stopping criterion is met. Of particular interest are the PoI and EI statistical criteria which are widely used for single-objective optimization [10, 24] . Hence, it may be useful to extend the concept of the PoI and EI directly to multiobjective optimization. Multiobjective versions of the PoI and EI are defined for an arbitrary number of objective functions in Sects. 3.2 and 3.3. For ease of notation in the forthcoming sections, the output of all the Kriging models can be considered as mutually independent Gaussian random variables Y j (x), The associated probability density function φ j and cumulative distribution function Φ j of Y j (x) are compactly denoted as, Given a set of n points X as in (2), a Pareto set P can be constructed that comprises v ≤ n Pareto-optimal (non-dominated) solutions, Each solution f(x * i ) is a vector that contains the objective function values for an associated input point 3.2 Probability of improvement (PoI) Evaluating the probability that the objective function values of a new input point x are located inside a well-defined region A in the objective space requires a multidimensional integration over that region. Naturally, several variants of the multiobjective PoI can be constructed as the concept of improvement is ambiguously defined in the context of multiobjective optimization. This is reflected in the selection of the integration region A, e.g., A can be the non-dominated part of the objective space or A can be the region in the objective space that solely extends (and does not dominate) the Pareto set, etc. In general, the probability that a new input point x yields improvement over the Pareto set P is denoted by the PoI P[I ], To evaluate Eq. (10), the integration area A can be decomposed into q (hyper-)rectangular cells, which yields a finite summation of contributing terms, see Fig. 2a . The lower and upper bound [l k , u k ] of each cell, for k = 1 . . . q, will be computed in Sect. 3.4. While the cells can be chosen to disjointedly cover the integration area A, the algorithm described in Sect. 3.4 decomposes the region A in overlapping cells. In this case, cells may negate the overlapping contribution of other cells by subtraction, denoted by the ± symbol in Eq. (11). While the PoI criterion is already quite useful and insensitive to the scaling of the objective functions, it does not, necessarily, encourage the generation of a uniform Pareto set. The EI quantifies the amount of improvement using an improvement function I (y, P) and, thus, prefers solutions that are lying farther from existing members of the Pareto set. The EI integral is defined as, J Glob Optim (2014) 60:575-594 581 In contrast to the various types of the PoI criteria (e.g., dominating 4 Pareto points or extending the Pareto set as in Sect. 3.2), it arguably makes more sense to only integrate the EI criteria over the region A corresponding to the non-dominated part of the objective space. The improvement function will automatically prefer new points that dominate the most points within the Pareto set P (the largest improvement). When no such points are found, the improvement function encourages the selection of points that extend the Pareto set P in an uniform way. Consequently, the design of the improvement function for the EI is crucial in identifying an optimal and uniform Pareto set. A good theoretical overview of different types of EI is given by [36] , including work on scalar improvement functions [16, 26] as well as using the single-objective EI in a multiobjective setting [23, 28] . Below we focus on evaluating the Euclidean distance-based EI [26] as well as the hypervolume-based EI [16, 17] efficiently for many objectives. In addition, a simplified version of the hypervolume-based EI is proposed that is significantly cheaper to compute. The hypervolume metric (or S-metric) [45] is widely used in multiobjective optimization to assess the quality of a Pareto set or to drive multiobjective optimization algorithms [5] . The hypervolume indicator H(P) denotes the volume of the region dominated by the Pareto set P, bounded by a reference point r which needs to be dominated by all points of the Pareto set, see Fig. 2a . Larger values of the hypervolume indicates better Pareto sets. Moreover, the exclusive hypervolume (or hypervolume contribution, see Fig. 2b ) of a Pareto set P relative to a point p is defined as, H exc measures the contribution (or improvement) of the point p to the Pareto set P and, hence, can also be used to define a scalar improvement function, namely, Subsequently, the integral of the hypervolume-based EI is, where A is the non-dominated region (bounded by the reference point r). Initially it was suggested to approximate the hypervolume-based EI using Monte Carlo techniques [17] . Recently, Emmerich et al. [16] proposed a method to calculate it exactly for an arbitrary number of dimensions by decomposing the non-dominated region into a set of cells as is also done in this work. Unfortunately, in [16] the proposed mathematical expressions assume that the non-dominated region is decomposed into an uniform grid of cells based on the Pareto set, see the cells bounded by the dashed lines in Fig. 2a . Hence, the number of cells required to evaluate the criterion scales at least exponentially with the number of Pareto points and objectives. Moreover, for each cell a separate hypervolume calculation needs to be done and, hence, it is infeasible to apply the method for three objectives or higher. This 123 582 J Glob Optim (2014) 60:575-594 work develops a new mathematical expression for the hypervolume-based EI that alleviates some of its computational complexities by decomposing the non-dominated region into a much smaller set of cells as well as removing the separate hypervolume calculations. In contrast to other statistical criteria, the proposed expressions for the hypervolumebased EI requires the non-dominated region to be covered by a set of disjoint cells, then the hypervolume-based EI can be written in closed form as, where I C k denotes the improvement contribution of cell k, namely, For each cell k = 1 . . . q contributions are calculated per dimension as follows, There are several ways a cell k can contribute to the I C k depending on its position relative to cell k, see Fig. 3 . Cell k is represented by the dark shaded cell. The pairs (·, ·) inside each cell k = 1 . . . q denote the type of contribution (A, B or C) per dimension. A, B and C refer to the Equations of the piecewise function G, see the previous Table. If a cell is completely covered by cell k in any one dimension, its length in that dimension is always included (type B). If a cell is only partially covered by the cell k, the contribution is divided into two parts: the integration when the cell is possibly fully covered (first line of A; zero if a j ≤ l j ) and when the cell is partially covered (second and third line of A). If a cell is not covered in any The correctness of this algorithm has been verified by an extensive numerical comparison against the publicly available code of Emmerich et al. [16] for the two objective case. In addition, the two and three objective cases have been verified using Monte Carlo methods. Regardless of the fact that the new procedure is already significantly cheaper than the method proposed by Emmerich et al. [16] , this hypervolume-based EI is still more expensive to evaluate than other statistical criteria. This is due to the computation time being more sensitive to the number of cells as well as the reliance on a binary partitioning of the nondominated region into disjoint cells, which requires more cells to cover the integration area than the Walking Fish Group (WFG) algorithm explained in Sect. 3.4. Hence, a simplification of the hypervolume-based EI is proposed in the next section which can be evaluated using the WFG algorithm. Inspired by the definition of the hypervolume-based EI [16], the hypervolume-based PoI can be written as the product of the improvement function I (μ, P) and the PoI P[I ], and so the advantages of using the hypervolume contribution can be preserved while significantly reducing the overall computational complexity, where I (μ, P) is defined as in (14) and μ = (μ 1 (x) , . . . , μ m (x)) is a vector that contains the prediction of the Kriging models of each objective function for a point x. In effect, the prediction variance is not taken into account anymore for the improvement function, in contrast to E I hv , as it is moved outside of the integral. The integration area A of P[I ] corresponds to the non-dominated region and, hence, a closed-form expression of the hypervolume-based PoI can be derived from the same set of cells used to evaluate P[I ], see Fig. 2b , namely, 123 584 J Glob Optim (2014) 60:575-594 where, Similarly to the hypervolume-based PoI, Keane [26] defines the EI as the product of the PoI P[I ] and an Euclidean distance-based improvement function. Let f c be the solution in P that is located closest to the centroidŷ(x) of the P[I ] integral, where the weight vector w = (w 1 , . . . , w m ) is used to scale the objective functions into the same range and,ŷ with the centroidŷ over arbitrary integral bounds [l, b] defined by, . Lastly, E euclid [I ] for an input vector x is defined as, Like all other EI criteria the integration area A of P[I ] andŷ(x) is the non-dominated region. In order to evaluate these statistical criteria efficiently, one or more integrals need to be evaluated over an integration area A. As A is non-rectangular and often irregularly shaped, especially for a higher number of objective functions, the integral must first be decomposed into a sum of k integrals over rectangular cells. While these cells can be identified analytically upfront for two objectives [26] or one can use the most fine-grained cells possible (for a total of q ≈ (v + 1) m cells; [16] ), this becomes rather prohibitively complex and cumbersome for a higher number of objective functions (>2). J Glob Optim (2014) 60:575-594 585 Instead, the authors propose to decompose the integration area in as few cells as possible using an efficient computer algorithm, i.e., each cell encompasses a large part of the integration area. A straightforward approach to determine the required bounds of the cells for the evaluation of the criteria is to use binary partitioning [11], see Fig. 2a . While this approach is quite flexible as it allows to identify different kind of integration areas (e.g., leading to several variants of statistical criteria), it becomes prohibitively expensive as the number of objectives exceeds four. Nonetheless, by terminating the binary partitioning early the statistical criteria can still be approximated fairly well for a higher number of objectives. However, the focus of this work is to improve the performance of the exact evaluation of the criteria. To that end, it makes sense to take advantage of the numerous algorithms for calculating the hypervolume. Formally, the hypervolume is the Lebesgue integral, where A is the region dominated by the Pareto set P and bounded by some reference point r. As the statistical criteria are integrals evaluated over a similar area, but using a different integrand, the idea is to adapt a hypervolume routine and retrieve the integration area A as a set of cells instead of immediately calculating its (hyper)volume. Exact algorithms [3, 4, 31, 43] for calculating the hypervolume as well as approximations [1, 2, 7, 8, 20, 28] and alternative versions of the hypervolume problem, e.g., finding the Pareto point(s) that contributes least to the hypervolume [9], have been suggested in literature. While the algorithm proposed by Beume et al. [3] has the best worst case complexity, the WFG algorithm [39] is actually faster on most practical optimization problems and, hence, is adapted in this work to evaluate the statistical criteria. The basic WFG algorithm operates by defining the hypervolume as a sum of exclusive hypervolumes, originally introduced by [18] , where each exclusive hypervolume in the summation corresponds to one of the slices bounded by dashed lines in Fig. 4 . At first sight this may look expensive as by definition the exclusive hypervolume (13) itself requires two separate hypervolume calculations. Fortunately, several optimizations can be made based on the following main ideas, see Fig. 4 . -Slices to the left of the contributing point p (i.e., with smaller values than p in the current objective) can be discarded as they contain no hypervolume that is dominated by p. -If p is dominated in the remaining (unprocessed) objectives it dominates no more exclusive hypervolume. Advantage of these insights can be taken by rewriting the exclusive hypervolume from (13) as H exc (p; P) = H({p}) − H(P ), where P = {limit(p, q)|q ∈ P} and limit (p, q) = (max{ p 1 , q 1 }, . . . , max{ p m , q m }). The first term H({p}) is simply the volume of the cell bounded by p and the reference point r, the second term H(P ) is a recursion where the hypervolume is calculated for P . P is obtained by taking the maximum (limit) of the contributing point and for each point in P, i.e., the points are projected on p, see Fig. 4 . Obviously, many points will be introduced that are dominated by some other point in P . These points do not contribute anything to the hypervolume and can be removed by taking the non-dominated subset of P before continuing calculation. It is this last step that significantly  One step in the WFG algorithm to calculate H exc (f 2 = p; {f 3 , f 4 }), denoted by the shaded region (f 1 has already been processed). The exclusive hypervolume slice is efficiently calculated as the volume of the cell bounded by p and r minus the hypervolume of a reduced Pareto set, represented by the squares, where all points are limited by the contributing point p. This creates many dominating points which can be removed before continuing calculation. In this illustration the volume of the cell bounded by f 3 and r is subtracted to calculate the exclusive hypervolume improves the performance as it has been shown that most datasets already lose over 50 % of their points after only one recursion. Note that if the contributing point p has generally higher objective values then likely more points can be pruned. Hence, it makes sense to process the worse Pareto points first by sorting the Pareto set P descending on the current objective before each iteration. Moreover, now the hypervolume calculation can also be sliced on the sorted objective each recursion, hence the use of the terms current objective and unprocessed objectives in the previous paragraph. Each recursion the current objective is sorted and is not considered in subsequent calculations of H excl , for more information on this optimization please see [38] . This is similar to the Hypervolume by Slicing Objectives algorithm (HSO; [40] ). A final optimization is the base case for two objectives, for which we can easily calculate the hypervolume in O(v) assuming the Pareto set is sorted on the current objective. The base case for two objectives is similar to the algorithm proposed by [26] to calculate the Euclidean distance-based expected improvement and, hence, the performance of the WFG algorithm is at least equal to it. The WFG algorithm is easily adapted to keep a record of a cell's lower-and upperbound instead of calculating its hypervolume. The adapted WFG algorithm, in case of minimization, identifies the cells that are dominated by the Pareto set. Naturally, the algorithm can be modified to find the dominated region by temporarily viewing it as a maximization problem. However, for the evaluation of the statistical criteria, especially the EI, it is much more useful to identify the non-dominated region. This is achieved by subtracting the cells obtained from the adapted WFG algorithm from the cell that covers m (possibly bounded by a reference point r when using one of the hypervolume-based criteria), see Fig. 2a . Hence, the adapted WFG algorithm can be used for statistical criteria where the integration area A corresponds to either the dominated or non-dominated region. Should other integration areas be required, e.g., for calculating the probability that a new point dominates at least two Pareto points, more flexible but slower methods such as binary partitioning [11] can be used. for an arbitrary input point x is also included for the 2D case. In contrast to the proposed algorithm, the exi2D algorithm needs to be run multiple times during optimization of the criteria After q << (v + 1) m sets of cells (= integral bounds) have been identified the actual PoI and EI statistical criteria can be evaluated using Eqs. (10), (20), (16) or (18). While evaluating the criteria the point f max + is replaced by (∞, . . . , ∞) or r depending on the criterion. A plot with the practical computation time and the number of the cells is shown in Fig. 5a , b, applying the adapted WFG algorithm to sets of Pareto points randomly drawn from the first quadrant of a unit sphere (taking the mean values of 1,000 repetitions). The computation time of the cells poses no problem for the evaluation of the criteria (well within the seconds range). The limited factor of the EMO algorithm is the number of cells the integration area A is decomposed into, as for each cell the corresponding PoI or EI equations needs to be evaluated which can become prohibitively expensive when many cells are required to identify the integration area A. The decomposition of the integration area A into cells and the actual evaluation of the PoI and EI criteria is separate in the sense that the cells only need to be identified once every sampling iteration, and then only if the (intermediate) Pareto set has changed with respect to the previous iteration. Afterwards, the PoI and EI criteria can be evaluated multiple times 123 588 J Glob Optim (2014) 60:575-594 for a point x, e.g., during optimization, using the same set of cells. Lastly, the hypervolume indicator is also easily obtained by summation of the volume of the cells. Examples A good set of configurable multiobjective benchmark problems has been proposed by Deb et al. [15], of which four benchmark functions are chosen and adapted slightly to benchmark the EMO algorithm. A summary of the selected benchmark functions is found in Table 1 . For a complete description of the benchmark functions the authors refer to [15] . All benchmark functions are configured to have six input parameters. Specifically, the first example is the DTLZ1 function with three objective functions where the Pareto front lies on the plane y 1 + y 2 + y 3 = 1. The second example is the DTLZ2 function with three objective functions where the Pareto front is the first quadrant of an unit sphere centered on the origin. The third example is the DTLZ7 function with four objective functions which has 2 m−1 = 2 4−1 = 8 disconnected Pareto-optimal regions in the objective space. The last example, the DTLZ5 function configured to have six objective functions, is similar to DTLZ2 except that the Pareto front is just one slice of the unit hypersphere, i.e., the Pareto front is a (densely populated) curve in a m = 6 dimensional objective space. An initial set of 65 samples is generated by a near-optimal maximin Latin Hypercube Design (LHD; [12] ). Subsequently, a statistical criterion is optimized for each iteration to select the next point to evaluate. The criterion is optimized using a combination of Monte Carlo sampling and a local search. Specifically, 20 × n Monte Carlo candidate points are generated and evaluated on the criterion. The best Monte Carlo candidate is further refined using Matlab's fmincon optimizer. Various configurations of the EMO algorithm are applied on the benchmark functions. In particular, EMO is configured with the E I euclid and P hv criterion together with Kriging models using the Matérn correlation function [32] with ν = 3 2 and a constant regression function (M = 1 and F = 1). The hyperparameters of the Kriging models are optimized using SQPLab [6] (http://www-rocq.inria.fr/ gilbert/modulopt/optimization-routines/ sqplab/sqplab.html), utilizing likelihood derivative information. The optional weight vector w of the E I euclid criterion is set to 1 for all benchmark functions, except for DTLZ7 where w = (1, 1, 1, 0.02) to scale the last objective function into the same range as the other objective functions. Furthermore, the EMO runs of the E I euclid criterion are repeated with Kriging models using the Gaussian correlation function, these runs are denoted by E I gauss euclid in the results. J Glob Optim (2014) 60:575-594 589 Lastly, extra EMO runs are configured for the DTLZ1 and DTLZ2 functions using the expensive E I hv criterion with Kriging models using the Matérn correlation. The E I hv criterion is not applied to the DTLZ5 and DTLZ7 problems. Each of the in total 14 EMO runs is repeated 10 times for statistical robustness and halts when the sample budget is met, namely, 250 samples. The EMO runs are compared against the NSGA-II, SPEA2 and SMS-EMOA evolutionary algorithms with a varying population size and maximum number of generations. The first run is configured with a population size of 25 and a maximum number of generations of 10 (total sample budget 250) and the second run is configured with a population size of 50 and a maximum number of generations of 50 (total sample budget 2,500). The remaining parameters have been left to their default values. Similarly to the EMO runs, the evolutionary algorithm runs are repeated 10 times. Beside assessing the performance of the algorithms using the hypervolume metric, the convergence measure is used too. The convergence measure is the mean distance of every point of the Pareto set to the closest Pareto point of the known Pareto front. In this work the known Pareto fronts are represented by 100,000 Monte Carlo points. Results for the benchmark functions have been summarized in Table 2 . Note that the differences on the hypervolume metric are more significant than they appear because of the conservative choice of the reference point r (needed to accommodate the results of all test configurations). In general, it is seen that the EMO runs have better performance than the MOEAs in terms of hypervolume score for most functions except for DTLZ1. After a closer examination it is observed that the accuracy of the Kriging models of DTLZ1 for most statistical criteria is sub-optimal. In particular, the first objective function is difficult to approximate using the Kriging models, see Fig. 6a . A plot of the final Pareto sets generated of the DTLZ2 problem is shown in Fig. 7 . It is seen that the hypervolume-based criteria emphasizes the edges of Pareto front more while leaving a small gap between the edge and the inner portion of the Pareto front. This is not unlike the DTLZ2 results as reported in [5] and is due to the nature of the hypervolume indicator. Logically, the farther away the reference point is located, the larger the exclusive hypervolume will be for points lying on the edge of the current Pareto set (as the exclusive hypervolume is then solely bounded by the reference point). Further research is needed to determine the influence of the choice of reference point r on the statistical criteria [1] . While the EMO algorithm outperforms the MOEAs on the hypervolume indicator on most problems, there are some limitations. The EMO algorithm, and other MOSBO techniques, rely on the quality of the surrogate model to guide the selection of new expensive data points. The Kriging models do not have to be accurate at the start of the algorithm when using the EI and PoI criteria, but they should be able to capture the behavior of the objective functions sufficiently well when enough samples become available, which might not always be the case (see Fig. 6 and the DTLZ1 results). Furthermore, the construction of the Kriging models and the evaluation of the statistical criteria comes at a computational cost, similar to the computational cost of MOEAs that rely on the hypervolume (i.e., SMS-EMOA), which might limit the practical usage of the EMO algorithm for some (less expensive) optimization problems. Specifically, the construction of the Kriging models and the thorough optimization of the statistical criteria make the EMO algorithm more expensive than SMS-EMOA. Moreover, 123 590 J Glob Optim (2014) 60:575-594 the EMO algorithm requires the storage of the integral bounds in memory, which can be prohibitively expensive for a higher number of objectives. The authors presented the EMO algorithm, which uses multiobjective versions of the PoI and EI to identify the Pareto front with a limited sample budget. Different configurations of the EMO algorithm are compared against the well-known SPEA2, NSGA-II and SMS-EMOA evolutionary methods with promising results. In theory an arbitrary number of objective functions can be handled. However, in practice due to the nature of the multiobjective EI and PoI statistical criteria the EMO algorithm also does not escape the curse of dimensionality (no-free-lunch theorem) with respect to the number of objective functions and number of Pareto points. Nevertheless, the EMO algorithm can be applied to problems up to eight objectives. Furthermore, it should be noted that in practice a high number of objectives are often aggregated to a smaller number of objectives (≤4) using, e.g., weighted sums, as such high-dimensional spaces are not easily interpretable for designers. Future work will focus more on exploring the key benefits of the EMO algorithm on various industrial applications and benchmark problems. Preliminary results show that the benefits of the EMO algorithm are more pronounced when limited to a very low number of function evaluations (<100). Moreover, the benefits of using statistical criteria, which balance exploration and exploitation automatically, may manifest more for problems with a complex optimization landscape and multiple local Pareto fronts. In addition, future work will focus on minimizing the number of cells and on an iterative update scheme for the cells, which will be considerable more efficient than recalculating the cells almost each iteration. Indirectly, a speedup can also be achieved by selecting multiple update points at a time. Finally, it may also be worthwhile to investigate the use of approximated statistical criteria, namely, adapting hypervolume approximation routines for decomposing the integration area into cells.
Recently, Couckuyt et al. REF published an algorithm that is faster based on empirical tests.
Fast calculation of multiobjective probability of improvement and expected improvement criteria for Pareto optimization
Journal of Global Optimization
Mathematics Computer Science