A survey on handling computationally expensive multiobjective optimization problems with evolutionary algorithms

Evolutionary algorithms are widely used for solving multiobjective optimization problems but are often criticized because of a large number of function evaluations needed. Approximations, especially function approximations, also referred to as surrogates or metamodels are commonly used in the literature to reduce the computation time. This paper presents a survey of 45 different recent algorithms proposed in the literature between 2008 and 2016 to handle computationally expensive multiobjective optimization problems. Several algorithms are discussed based on what kind of an approximation such as problem, function or fitness approximation they use. Most emphasis is given to function approximation-based algorithms. We also compare these algorithms based on different criteria such as metamodeling technique and evolutionary algorithm used, type and dimensions of the problem solved, handling constraints, training time and the type of evolution control. Furthermore, we identify and discuss some promising elements and major issues among algorithms in the literature related to using an approximation and numerical settings used. In addition, we discuss selecting an algorithm to solve a given computationally expensive multiobjective optimization problem based on the dimensions in both objective and decision spaces and the computation budget available.


Introduction
Many engineering problems have multiple objectives to be optimized, and these objectives are typically conflicting in nature, i.e., improvement in one objective is possible only by allowing deterioration of at least one of the other objectives. These kinds of problems are known as multiobjective optimization problems (MOPs). Because of the conflicting nature, there typically does not exist one optimal solution, Communicated  but multiple so-called Pareto optimal solutions. The set of all Pareto optimal solutions in the objective space is called a Pareto front. In many problems, explicit formulations of objective or constraint functions are not known and such functions are called black box functions. Usually, problems involving such functions need a long time to be solved, e.g., problems involving computational fluid dynamics simulations utilizing finite element algorithms take substantial time to obtain one solution. These are examples of problems that we refer to as computationally expensive multiobjective optimization problems.
In the last few decades, evolutionary algorithms (EAs) (Coello et al. 2007;Deb 2001) have been widely used for solving MOPs because of their advantages such as obtaining a set of nondominated solutions in one solution process, ability to handle problems with multiple local and nonconvex Pareto fronts, ability to easily deal with different kinds of variables (such as binary, real, integer or mixed) and no assumptions set on convexity and differentiability of objectives and constraints involved. Despite of these advantages, EAs do not guarantee convergence to optimal solutions. Moreover, they are often criticized as they consume many function evaluations which increases the computation time. This concern is particularly relevant when dealing with computationally expensive problems. It is therefore required to adapt EAs in a way that they can be used to obtain solutions in less computation time without too much reduction in the quality of solutions. In this paper, a survey of different algorithms proposed in the literature to handle MOPs with computationally expensive objective functions is presented. For other challenges related to solving MOPs, see Coello et al. (2007), Deb (2001) and Miettinen (1999).
Several algorithms have been proposed reduce the computational cost while solving MOPs using EAs. Different surveys exist in the literature (Jin 2005(Jin , 2011Santana-Quintero et al. 2010;Knowles and Nakayama 2008) on this topic. One of the popular approaches mentioned in these surveys is the use of approximation, especially function approximation or metamodeling to handle computationally expensive problems. The potential of metamodeling techniques has been widely studied in the machine learning community (Forrester and Keane 2009). Many algorithms exist in the literature which use a combination of these techniques and evolutionary algorithms.
In Jin (2005), different ways of using approximation such as problem, function and fitness approximation were discussed based on their use in the literature. Additionally, various issues such as how evolutionary algorithms can benefit from these approximations and what kind of approximation to be selected, etc., were also discussed. It is to be noted that fitness approximation is also used as evolutionary approximation in the literature (Jin 2005). In Santana-Quintero et al. (2010), the main focus was to discuss various algorithms proposed in the literature before the year 2008 based on the classification of Jin (2005). Some real-world applications were also mentioned, where different algorithms had been applied to reduce the computation time. In Jin (2011), the main focus was to discuss the recent developments for using function approximation in reducing the computation time. In that survey, different issues such as strategies for managing approximation, selection of individuals to be evaluated using approximation functions were discussed. In addition, the author also mentioned the potential of using function approximation in dynamic, robust and constrained optimization problems. A taxonomy of metamodel-based optimization algorithms was provided in Horn et al. (2015). In Knowles and Nakayama (2008), several algorithms based on the management of using function approximation (e.g., Kriging, radial basis function) were discussed. In addition, the importance of interactive algorithms with EAs was detailed. Two real-world applications were also mentioned, which were solved by two different algorithms, ParEGO (Knowles 2006) and the algorithm proposed in Nakayama et al. (2006).
In this survey, we extend these surveys and discuss 45 algorithms based on the classification of Jin (2005) from the year 2008-2016 published in different journals and conference proceedings in English. We also found that most of the algorithms use Kriging when compared to other function approximation techniques and therefore classify different function approximation-based algorithms into Kriging-and non-Kriging-based algorithms. This survey is also different from other surveys in following ways.
1. Algorithms using function approximation or surrogates are emphasized as they are more widely used. 2. Classification of function approximation-based algorithms further into Kriging-and non-Kriging-based algorithms shows the wide applicability of Kriging. 3. Algorithms are described in reference to a general function approximation framework which will provide the reader an understanding for using any function approximation in EAs. 4. The efficiency of different algorithms in terms of reducing computational cost or number of function evaluations is emphasized. 5. Various shortcomings in algorithms are observed and discussed, e.g., handling constraints, dimensions (in both decision and objective spaces) of the problem solved, training time 6. Some promising elements are also identified which can be helpful in overcoming the limitations of several issues, e.g., efficient management of function approximation technique to handle more than three objectives. 7. Some guidelines are also provided in selecting a particular algorithm based on the characteristic of the problem to be solved, dimensions in both objective and decision spaces and the budget available.
In the rest of this paper, in Sect. 2, some relevant concepts are described which are frequently used when reducing the computational burden. In Sect. 3, different algorithms are discussed based on the steps of general approximation framework. In addition, different algorithms are compared based on the type of approximation, evolutionary algorithm, evolution control and characteristics of the optimization problem solved. A brief discussion on various issues and using promising elements related to the use of different algorithms is presented in Sect. 4. Finally, conclusions are drawn in Sect. 5 along with future research directions.

Basic concepts and terminology
We consider multiobjective optimization problems of the form (Miettinen 1999): with k(≥ 2) objective functions f i (x) : S → . The vector of objective function values is denoted by f (x) = ( f 1 (x), . . . , f k (x)) T . For the simplicity of presentation, we assume that all the objective functions are to be minimized.
If some objective function f i is to be maximized, it is equivalent to minimize − f i . The (nonempty) feasible region S is a subset of the decision variable space n and consists of decision variable vectors x = (x 1 , . . . , x n ) T that satisfy all the constraints. The image of the feasible space S in the objective space k is called the feasible objective set denoted by Z . The elements of Z are called feasible objective vectors denoted by f (x) or z = (z 1 , . . . , z k ) T , where z i = f i (x), i = 1, . . . , k, are the objective function values. As discussed in the introduction, objective functions in a MOP are typically conflicting in nature, and, thus, there is no single well-defined optimal solution but a set of socalled Pareto optimal solutions exists. A decision vector x * ∈ S is Pareto optimal if there does not exist another decision vector x ∈ S such that f i (x) ≤ f i (x * ) for all i = 1, . . . , k and f j (x) < f j (x * ) for at least one index j. An objective vector is Pareto optimal if the corresponding decision vector is Pareto optimal. A Pareto optimal set consists of all Pareto optimal solutions in the decision space, and a Pareto front consists of all Pareto optimal solutions in the objective space.
A solution x 1 ∈ S is said to dominate another solution x 2 ∈ S, denoted by x 1 x 2 if f i (x 1 ) ≤ f i (x 2 ) for all i = 1, . . . , k and f i (x 1 ) < f i (x 2 ) for at least one i ∈ {1, . . . , k} .
Ideal and nadir objective vectors represent bounds for objective function values in the Pareto front. Components of an ideal objective vector z * ∈ k are determined by minimizing each objective function individually, that is, z * i = minimize x∈S f i (x). A nadir objective vector consists of upper bounds of objective functions in the Pareto front. It is usually difficult to obtain for problems with more than two objectives. Several ways of approximating it have been proposed in the literature (see, e.g., Deb et al. 2010;Miettinen 1999). In what follows, we define concepts relevant to the present survey.
A metamodel is an approximation of some computationally expensive element in a multiobjective optimization problem. A metamodel replaces the computationally expensive element by an element which consumes less computation time. There can be different computationally expensive elements while solving a MOP, for example, original functions or constraints, hypervolume. A surrogate is often used as a synonym for a metamodel and the process to build a metamodel is known as metamodeling or function approximation. Therefore, what to approximate using a metamodel is an important concern in reducing the computation time. Neural networks (Mitra and Majumder 2011;Yuan and Guangchen 2009), radial basis functions (Qasem et al. 2013), support vector regression (Aytug and Sayin 2009) and Kriging (Husain and Kim 2010) are some examples of commonly used metamodeling techniques.
An ensemble of metamodels means using more than one metamodel to reduce the computation time and usually there are two ways to use an ensemble of metamodels in the literature. In the first one, one metamodel having the highest accuracy among different metamodels is selected to evaluate individuals (Singh et al. 2010). In the second one, different metamodels are used to evaluate individuals with a different weight coefficient (ω j ) assigned to them. For instance, in Lim and Jin (2010), the predicted fitness value from an ensemble of m metamodels was defined as: using jth metamodel and w i , i = 1, . . . , k is the weight vector assigned to each objective function. A weighted sum method was used in Lim and Jin (2010); however, other scalarizing methods such as -constraint method, achievement scalarizing function (Miettinen 1999) can also be used to convert a multiobjective optimization problem into a single objective optimization problem. Fitness values, F ens (x) is approximated by m metamodels and ω j is the weight assigned to the fitness value of the jth metamodel. The fitness value of a metamodel is assigned a larger weight if it is found to be more accurate than other metamodels, and the accuracy can be obtained by statistical measurements such as root mean square error (Tenne and Armfield 2008).
An evolution control or model management (Jin et al. 2002) is a strategy to manage metamodels in an EA. There are two different ways to manage metamodels, fixed and adaptive evolution control. In a fixed evolution control, metamodels are used for a prefixed number of generations and updated afterward with a predefined criterion. On the other hand, in an adaptive evolution control, the frequency of using the metamodel is adjusted according to the accuracy of the metamodel. Therefore, the model management is also concerned with when to update the metamodel.
In fitness inheritance (Sastry et al. 2001;Zheng et al. 1997), the fitness values of offspring are evaluated using fitness values of the parents. In fitness imitation (Kim and Cho 2001), individuals are grouped into several clusters, e.g., using K-means clustering (Jain et al. 1999) in the decision space and only those individuals are evaluated which represent the clusters (e.g., individuals closest to the centroid of each cluster). The fitness values of other individuals are estimated using the fitness values of the representative individuals.
Next, we discuss the general function approximation framework and describe various algorithms based on the steps of the framework.

Approximation-based algorithms
In approximation-based algorithms, the computationally expensive element of the problem is replaced with an approximation which consumes less computation time. As classified in Jin (2005), an approximation can be applied in three ways in multiobjective optimization problems: problem, function and fitness approximation. In problem approximation, the original problem is replaced with a simplified problem which is faster to solve. In function approximation, an approximation of a computationally expensive function is formed which is faster to evaluate. On the other hand, in fitness approximation, the fitness value referring typically to the function value of an individual is derived from the fitness values of the existing evaluated individuals in its vicinity. The function approximation is more widely used than other approximations and algorithms applying it are discussed next in Sect. 3.1. Approaches using problem and fitness approximation are discussed in Sect. 3.2.

Function approximation
As said, function approximation is the most commonly used approach among approximation-based algorithms and there, an explicit or an implicit approximation of a computationally expensive, function is formed, which is faster to evaluate. In what follows, we refer to the functions of the original, computationally expensive problem as original functions. In the literature, metamodel is often used for all objective functions; therefore, all objective functions are assumed as computationally expensive. Neural networks (Kourakos and Mantoglou 2013;Mitra and Majumder 2011;Azzouz et al. 2014;Lim and Jin 2010), Kriging (Knowles 2006;Jang et al. 2009;Ponweiser et al. 2008a;Li et al. 2008) and polynomial regression (Husain and Kim 2010;Singh et al. 2010) are examples of algorithms used for function approximation.
The general steps for using function approximation in an EA can be divided into two stages consisting of ten steps, and most of the algorithms discussed in this section follow these steps. In what follows, we present a function approximation framework that captures the core of the EA proposed in the literature utilizing function approximation. Additionally, it is assumed that one type of metamodel is used for all objective and constraint functions involved although it is possible to use different metamodels for different objective functions. Stage 1 1. Initialize the population either randomly or using a sampling method. 2. Evaluate the individuals of the population using the original functions and add them to an archive. 3. If a prefixed number of generations is completed, go to stage 2. 4. Use EA operators (selection, crossover and mutation) to create a new population and go to step 2.
Stage 2 5. Build or update a metamodel for each computationally expensive objective and constraint function using the individuals from the archive. 6. Use EA operators (selection, crossover and mutation) to create a new population. 7. For each objective and constraint function, evaluate the individuals of the new population either using the metamodel or the original function (fixed or adaptive evolution control strategy). 8. If a stopping criterion is met, select the nondominated individuals from step 7 as the final population and stop. Otherwise, continue. 9. Select individuals from step 7 for re-evaluation using the original functions if needed. 10. Add the individuals from step 9 to the archive and go to step 5.
In stage 1, a general evolutionary algorithm works, and in stage 2, an EA with a metamodel is used. In step 1, a population is initialized either randomly or using a sampling method such as Latin hypercube sampling (Mckay et al. 2000). Individuals of the population are evaluated using the original functions and the evaluated individuals are added to an archive in step 2. If a prefixed number of generations is completed in step 3, stage 1 is terminated and the archive is carried over to stage 2. Otherwise, a new population (offspring population) is generated using EA operators such as selection, crossover and mutation in step 4. A prefixed number of generations is required to obtain an archive of a fixed size. In most of the papers cited in this paper, there is no explicit criterion mentioned for choosing the prefixed number of generations. However, it is mentioned in Arias-Montano et al. (2012) that the prefixed number of generations depends on the dimensions of the problem (in both objective and decision spaces) and the budget of evaluations with the original functions. For instance, in three popular algorithms known as ParEGO (Knowles 2006), SMS-EGO (Ponweiser et al. 2008a) and MOEA/D-EGO (Zhang et al. 2010), a data set of size 11n − 1 was used for training the metamodels. We provide the size of the data set used in different algorithms in Table 1. After completion of stage 1, a metamodel is created for each computationally expensive objective and constraint function to work with the EA algorithm for evaluating individuals in stage 2. A metamodel is created using the individuals from the archive in step 5. In step 6, an offspring population is generated using EA operators. In step 7, either the metamodel or the original functions are used to evaluate individuals, i.e., a fixed or an adaptive evolution control strategy is used to manage the metamodel. In step 8, if a termination criterion such as maximum number of generations or function evaluations is met, nondominated individuals from the last population are selected as the final population. Otherwise, some individuals are selected from step 7 and re-evaluated using the original functions in step 9. There are different criteria mentioned in the literature to select individuals for re-evaluation, e.g., selection of non-dominated individuals, using expected improvement Jones et al. (1998), expected hypervolume improvement Ponweiser et al. (2008a). These individuals are then added to the archive in step 10 to update or retrain the metamodel. If the size of the archive is prefixed, some individuals (e.g., random or dominated individuals) are eliminated from the archive using a predefined criterion.

Challenges in using metamodels
Before going into the details of different algorithms, we mention here major challenges when applying function approximation in multiobjective optimization when compared to single objective optimization. It is not straightforward to use a metamodel with an EA because several challenges exist which affect the performance of the metamodel used. These challenges are also the main differences in several algorithms in using function approximation.
1. Using the metamodel In case of single objective optimization, often one metamodel is used for the approximation of the objective and constraint function. Using metamodels in multiobjective optimization is not that trivial. For example, one can use different metamodels for different objective functions, one metamodel for all objective functions, different metamodels for one objective function or a single metamodel for scalarized objective functions. The same options may arise for different constraints. In the literature, most of the algorithms use a single metamodel for all objective functions. In addition, there are no guidelines for using different metamodels for different objective functions. 2. How to update the metamodel Updating the metamodel is an important step in both single objective and multiobjective optimization. In case of single objective optimization, different strategies have been used for updating the metamodel, e.g., lower confidence bound (Dennis and Torczon 1995), probability of improvement (Ulmer et al. 2003), expected improvement (Sasena et al. 2002;Jones et al. 1998). In case of multiobjective optimization, two goals, convergence and diversity have to be achieved in contrast to single objective optimization. In case of multiobjective optimization, one of the simplest ways is to select the nondominated individuals for re-evaluation (Jin et al. 2000(Jin et al. , 2002 using the original functions. Alternatively, a clustering method such as K-means clustering (Jain et al. 1999) can be applied to cluster the individuals in the objective space and the representative individuals are re-evaluated using the original functions Gräning et al. (2007); Jin and Sendhoff (2004). The representative individuals can be the individuals closest to each cluster center or the best individuals (with highest fitness values) in each cluster. Individuals with low approximation accuracy can also be selected for re-evaluation (Branke and Schmidt 2005;Emmerich et al. 2002) and this selection may be effective in improving the approximation accuracy of the metamodel. In the literature, criteria for updating the metamodel in single objective optimization are also modified to be used in multiobjective optimization. For instance, expected improvement in Knowles (2006), Jeong and Obayashi (2005), Keane (2006), Couckuyt et al. (2014 and Wagner et al. (2010) and probability of improvement in Couckuyt et al. (2014) are used to update the metamodel in multiobjective optimization. In addition, expected hypervolume improvement is used in Emmerich et al. (2006Emmerich et al. ( , 2011 and Shimoyama et al. (2013) to select the individuals for re-evaluation so that the metamodel can be updated. It is to be noted that infilling criterion is sometimes used as a synonym for updating criterion in the literature. 3. Training time for the metamodel Training time is also an issue in case of single objective optimization but when moving to multiobjective optimization, training or updating time becomes more influential. For example, different metamodels for different objective functions may take a different amount of training data (archive size) which increases the training time when compared to single objective optimization. In the literature, unfortunately, most of the algorithms do not mention the training time for the metamodel. This is due to the reason that the training time of the metamodels is usually assumed to be negligible compared to objective function evaluations. However, it may happen that the training time of the metamodels is substantial and the whole aim of reducing computation time is jeopardized. In addition to the three challenges above, two challenges exist which are common to both single objective and multiobjective optimization. They can influence the performance of an EA while using metamodels. 1. Choice of the metamodel In the literature, there is very little guidance about the choice of the metamodel for approximation of computationally expensive functions. A metamodel is either selected randomly or due to its popularity in the area with which the problem is associated. For instance, in Kourakos and Mantoglou (2013), radial basis neural network was used as a metamodel for approximation of objective functions in coastal aquifer management problem and it was mentioned that neural networks were popular for groundwater applications. Kriging was used in Jang et al. (2009) due to its promising nature for building accurate global approximations. In Palar et al. (2016), different metamodeling techniques were used for different problems, i.e., radial basis function (RBF) with linear kernel function for ZDT problems and Kriging for UF ) problems. The reason mentioned was that RBF was used due to simplicity of the function landscape in ZDT problems and Kriging was used because functions in UF problems are highly nonlinear and test with RBF showed a very low accuracy when used for UF problems. However, the algorithm was developed for black box computationally expensive problems, where function landscape of the problems in question is not known a priori. Moreover, Kriging is most popular technique when compared to others because of its ability to provide uncertainty for the approximated values. To overcome the problem of choosing of a metamodel, an ensemble of metamodels described in Sect. 2 is also used in the literature (Singh et al. 2010;Lim and Jin 2010). However, there are still some open challenges related to the ensemble of metamodels such as what should be the criterion for choosing different metamodels or how different metamodels can be used simultaneously? 2. When to update the metamodel It is also an important issue to decide when to update the metamodel or in other words, is there any need to update the metamodel. For example, before selecting individuals in step 9, one can check whether the existing metamodel is accurate enough to predict objective function values in the next generation. If so, the metamodel is not updated and the existing metamodel is used for approximation. An offspring population has to be generated before updating the metamodel so that the metamodel accuracy can be measured for the individuals of the offspring population. There may be a possibility that the metamodel is never updated and the metamodel trained after stage 1 is used in all generations. In that situation, steps 9 and 10 are eliminated from stage 2. In the literature, there is very little focus on the question of when to update the metamodel. 3. Handling constraints While using metamodels in constrained problems, an appropriate set of solutions is needed to train them. Most of the metamodel-based algorithms in the literature are not developed to handle constraints. In addition, few algorithms which are developed to handle constraints assume that feasible solutions are available to train metamodels for constraint function. In many problem, obtaining a feasible solution is not trivial, e.g., C1-DTLZ1 problem (Jain and Deb 2014) with very small feasible region. In such instances, using metamodels for constraints can be very challenging. We provide more details of handling constraints while describing different algorithms.
Steps 5, 9 and 10 representing the first three main challenges are written in italics in the function approximation framework to indicate the steps, where algorithms using function approximation in multiobjective optimization mainly differ from each other. Other three challenges, selection of a metamodel, when to update the metamodel and size of the data to train the metamodel are not the major differences in the algorithms.
In what follows, 30 algorithms and from the literature are discussed utilizing the three main steps (5, 9 and 10) of the function approximation framework and attention is paid to their efficiency in reducing the computation time. These algorithms are classified according to the number of metamodels they have used, e.g., single metamodel, multiple metamodels independently or ensemble of metamodels. Among single metamodel-based algorithms, Kriging is used more often than other metamodels. As mentioned, the main advantage from Kriging is the uncertainty information besides the predicted objective and/or constraint function values. This uncertainty information can be further utilized in the algorithm. For instance, uncertainty information is used while updating the metamodel in Emmerich et al. (2011), Chugh et al. (2016a,  and Knowles (2006).
In contrast to single metamodel-based algorithms, some algorithms use metamodels independently. Moreover, some algorithms use ensemble of metamodels to solve the problem of choosing a metamodel. In what follows, three categories, algorithms based on single metamodel, algorithms based on multiple metamodels and algorithms based on ensemble of metamodels, are used to describe different algorithms. It was found that the number of papers in the literature belonging to the function approximation are largely skewed toward the first category, i.e., algorithms based on a single metamodel. Hence a sub-classification based on the usage of metamodel within different algorithms is devised to enhance clarity and readability. Thus, we further classify the algorithms using single metamodel into Kriging-and non-Kriging-based algorithm. Due to a limited number of articles in other categories, such a sub-classification is not used for others. It is also worth mentioning that most of the algorithms are based on dominance-based EAs and are generic for using any metamodel. At the end of this section, a comparison is presented based on which metamodel and type of EAs are used, what is the evolution control strategy and what are the characteristics of the optimization problems.
Parameter values used in these algorithms in stages 1 and 2 of the function approximation framework are presented in Table 1, where "NA" indicates that the parameter is not applicable to the algorithm and "not given" means that this information is not given in the reference. The table collects information of the population size in step 1, the size of archive in step 2, the prefixed number of generations in step 3 in stage 1 and the number of individuals for re-evaluation in step 9 in stage 2, as reported for example problems solved in the papers cited.
All parameters mentioned in Table 1 are important and can influence the performance of the algorithm. The population size |N P | is a critical parameter and several studies like  and Coello Coello and Lamont (2004) show the effect of the population size on the performance of evolutionary algorithms. The size of the archive |N A | mainly affects the training time of the metamodels, and as can be seen from the table, most of the algorithms do not have a fixed size archive. The size of data set |N I | used to train metamodels is different in different algorithms. One can adjust this parameter based on the resources or maximum number of function evaluations available. In addition, the size of the data set |N U | to update the surrogates is also important and should be decided based on the resources available. For instance, if it is possible to do parallel evaluations, one can select the size accordingly. An example is given in Knowles (2009), where the ParEGO algorithm could not be applied in doing experi-ments for drug design because the algorithm did not have an option to select multiple sample points at a time. The maximum number of function evaluations FE max is also important when comparing different algorithms. Many algorithms use different numbers of function evaluations, which makes it difficult to select an algorithm to solve a given problem. The number varies from 50 to 30,000 in the literature. Next, we classify algorithms which use only a single metamodel, which are further classified into Kriging-and non-Krigingbased algorithms.

Algorithms based on a single metamodel
In this subsection, we discuss algorithms using one metamodel for all objective and/or constraint function. As mentioned, to make the a clear structure of the paper and also due to wide applicability of Kriging models, we further classify this subsection into Kriging-based algorithms and non-Kriging-based algorithms. Further these algorithms are discussed year wise, i.e., starting from the year 2008.

Kriging-based algorithms
In this section, we discuss algorithms using Kriging.
A DACE model (Sacks et al. 1989) which is a Krigingbased approximation was used in Ponweiser et al. (2008a). This algorithm is known as SMS-EGO which is an extension of efficient global optimization (EGO) (Jones et al. 1998) for multiobjective optimization problems. In step 5, the metamodel is built for each objective function and to update it in step 9, one individual having maximum contribution to the hypervolume is selected for re-evaluation using the original functions. This individual is then added to the archive in step 10, and the metamodel is updated in the next generation. The size of the archive is not fixed in this algorithm.
SMS-EGO was tested on five benchmark problems (2-5 objectives and 3-6 decision variables) and compared with ParEGO (Knowles 2006) and the algorithm proposed in Jeong and Obayashi (2005). The unary hypervolume indicator (Zitzler and Thiele 1998), the R2 indicator (Hansen and Jaskiewicz 1998) and the unary epsilon indicator (Zitzler et al. 2003) were used as the comparison criteria. The proposed algorithm performed better than the other algorithms in all cases when hypervolume was used as the performance measure. In terms of other comparison criteria, it performed worse than the other algorithms for three problems.
In Li et al. (2008), Kriging was used as a metamodel with a modified version of NSGA-II. In this algorithm which is known as K-MOGA, the metamodel is used for approximating each objective function in step 5. To update the metamodel in step 9, domination status is measured for each individual evaluated using the metamodel. Individuals which change the domination status are re-evaluated using the original functions and added to the archive in step 10. To check the domination status, minimum of minimum distance (MMD) is measured. To calculate MMD, individuals are partitioned into two sets in the decision space, nondominated (x nd ) and dominated (x d ). MMD is calculated as, wheref is the predicted objective vector using the metamodel. MMD is then projected to each objective function axis to obtain MMDf i (for i = 1, . . . , k). Thereafter, a threshold s i (x) ≤ MMDf i is specified by using the standard deviation obtained (s i (x)) from the Kriging model for each objective function. The individuals which do not satisfy this threshold are re-evaluated using the original functions. The size of the archive is not fixed in this algorithm. K-MOGA was tested on five benchmark (two objectives and 3-6 decision variables) and two real-world problems (two objectives and 2-4 decision variables) and compared with MOGA. For comparison, nondominated individuals from both algorithms were visualized and similar solutions were obtained in fewer function evaluations with K-MOGA.
The same algorithm was extended to KD-MOGA in  with one additional element. In KD-MOGA, a fixed number of individuals is generated using constrained maximum entropy design after step 9, which is an extension of unconstrained maximum entropy design (Currin et al. 1998). These individuals are then added to the population for the next generation.
KD-MOGA was tested on the same set of problems as K-MOGA with one more real-world problem (two objectives, four decision variables and four constraints). It was compared with both EAs (modified version of NSGA-II) and K-MOGA using visualization of the nondominated individuals. A similar performance was obtained in fewer function evaluations with KD-MOGA.
Kriging was used as a metamodel in Jang et al. (2009). This paper is based on an algorithm called combined AASO-AAMO (adaptive approximation in single objective optimization-adaptive approximation in multiobjective optimization) (Yang et al. 2002) which uses a multiobjective genetic algorithm (Yang et al. 2002) as an EA. The metamodel is built for each objective and constraint function in step 5. To update the metamodel in step 9, a fixed number of individuals is selected for re-evaluation using the original functions and added to the archive in step 10. The individuals having the best value according to the maxmin distance design criterion (Johnson et al. 1990), i.e., the largest value of minimum distance from individuals in the archive in the decision space, are selected for re-evaluation. The size of the archive is not fixed in this algorithm.
The proposed algorithm was tested on a fatigue design MOP with two objectives, 17 decision variables and 11 constraints. It took 70 h to complete 25 generations but the efficiency of the algorithm was not compared to any EA without using the metamodel. The authors mentioned that it was practically impossible to do the optimization without using approximation.
In Zhang et al. (2010), Kriging was used as a metamodel with the decomposition-based EA MOEA/D (Zhang and Li 2007). The algorithm is known as MOEA/D-EGO. In MOEA/D-EGO, after evaluating a fixed number of individuals in steps 1-4, the metamodel is built for the scalarized objective function. Chebyshev scalarizing function (Miettinen 1999) is used to convert multiobjective optimization problem into single objective optimization problems. After using the metamodel in step 7, the expected improvement (EI) (Jones et al. 1998) is calculated for each subproblem. Expected improvement is then maximized (for each subproblem) using MOEA/D-DE ) for a fixed number of generations. In other words, the scalarized problem is changed into another problem to maximize EI. To update the metamodel in step 9, firstly, all individuals after the local search, which are different from the individuals in the archive (in the decision space) are selected. After this, Kmeans clustering is used to cluster the weight vectors (used in MOEA/D initially) associated with the individuals selected above. From each cluster, an individual with maximum EI is selected and re-evaluated using the original functions. These individuals are then added to the archive in step 10. Moreover, to reduce the training time, a fuzzy clustering (Bezdek 1981) is used for selecting a fixed number of individuals for training or updating the metamodel.
MOEA/D-EGO was tested on 12 benchmark problems (2-3 objectives and 2-8 decision variables) and compared with ParEGO (Knowles 2006) and SMS-EGO ). Hypervolume and inverted generational distance (IGD) were used as the comparison criteria. The proposed algorithm outperformed on seven problems when compared to ParEGO and performed similar to SMS-EGO. It was also compared with MOEA/D on two problems (two objectives and 2-8 decision variables) and outperformed with IGD as the performance criterion.
A multiobjective variable-fidelity optimization (VFO) algorithm was proposed in Zhu et al. (2013), where Kriging was used as a metamodel with NSGA-II. Initially, a simplified or approximated problem (having low-fidelity functions) is used to replace the original MOP in stage 1. This simplified problem is then solved for a prefixed number of generations using NSGA-II. A fixed number of nondominated individuals obtained after this step is re-evaluated using the original functions and stored in an archive. These individuals are selected based on a simple inter-individual distance metric in the objective space (steps 1-4). These individuals are then used to construct a Kriging model for each objective function in step 5 of stage 2. For the following fixed number of generations, individuals generated by crossover and mutation in step 6 are evaluated either using the approximated problem functions with Kriging model or just the approximated problem functions (step 7) depending on the error (Gano et al. 2006) in the Kriging model. Thus, an adaptive evolution control strategy is used to manage the metamodel. In other words, if the metamodel is not accurate for an objective function, the original function is used to evaluate individuals. To update the metamodel in step 9, a fixed number of nondominated individuals evaluated using the metamodel is selected for reevaluation using the original functions. The individuals are selected by computing the error from the Kriging model and added to the archive in step 10. The size of the archive is not fixed in this algorithm.
The VFO algorithm was tested on the ZDT1-3 benchmark problems (with two objectives and five decision variables) and a structural engineering problem (with two objectives and six decision variables). For the ZDT problems, the efficiency of the VFO algorithm was not mentioned in terms of computation time or number of function evaluations.
For the structural engineering problem, an exhaustive search was carried out with different values of decision variables. Around 65 million combinations of decision variable values were evaluated using both original and approximated problem functions. Nondominated individuals after this search were identified, and nondominated individuals obtained using original functions were used to compare results of different studies mentioned next. Fourteen studies were performed with changes in the values of three parameters: number of prefixed number of generations in step 3 of stage 1, number of generations before updating the metamodel of stage 2 and number of individuals selected for re-evaluation in step 9 of stage 2. Two out of fourteen studies were performed using NSGA-II without using a metamodel with high-and low-fidelity functions (to be called case 1 and case 2, respectively). Results from these studies were then compared with the results of exhaustive search with different criteria [inverted generational distance, error ratio, crowding distance (Deb 2001) and span (Lee et al. 2005)]. Case 1 gave the best results and in comparison with case 1, and one of the studies out of thirteen gave similar results (a graphical comparison was performed) in fewer function evaluations. Liu and Collette (2014) extended the VFO algorithm, where instead of one global metamodel, multiple local metamodels were used. The authors mentioned that local metamodels are used for high-dimensional problems in the objective space. K-means clustering is used in the decision space to partition the data and to build multiple local metamodels. Other details are the same as in VFO.
This algorithm was tested on six benchmark problems (two objectives and 2-30 decision variables) and one realworld problem (three objectives and six decision variables). For benchmark problem, nondominated individuals obtained from this algorithm and with VFO and NSGA-II were visu-alized. The authors mentioned that the proposed algorithm outperformed on these problems. In case of the real-world problem, the same performance criterion was used as in the VFO algorithm. In this case too, the proposed algorithm performed better than NSGA-II and VFO in the same number of function evaluations. Singh et al. (2014) used Kriging and the main focus was to apply different strategies for objective and constraint functions while updating the metamodel. In step 5, a metamodel is built for each objective and constraint function. To update the metamodel in step 9, for objectives, selection is performed using hypervolume-based probability of improvement (POI hv ) (Couckuyt et al. 2014) and for constraints, probability of feasibility (POF) (Forrester and Keane 2009) is used. After this, γ = POI hv × POF is obtained and a fixed number of individuals with highest γ values is selected for re-evaluation in step 9. These individuals are then added to the archive in step 10. The size of the archive is not fixed in this algorithm.
The proposed algorithm was tested on two real-world problems with two objectives, 2-3 decision variables and 5-7 constraints. This algorithm used the surrogate modeling MATLAB toolbox (SUMO) from Gorissen et al. 2010 and was not compared with any other algorithm. Mlakar et al. (2015) used Kriging models with a differential evolution based EA to approximate the objective functions. After building the metamodels in step 5, individuals were generated using differential evolution operator in step 6 and metamodels were then used for approximation in step 7. For each offspring in step 7, uncertainties of the approximated values were compared to its corresponding parent. In other words, if the uncertainty value of offspring was less than that of parent, it would selected and kept in the population. On the other hand, if the uncertainties of both offspring and parent are comparable, both are kept in the population. All the offspring thus selected were re-evaluated with the original function to update the metamodels in step 9 and added the training archive in step 10.
The algorithm was tested on 12 benchmark problem out of which three had constraints. However, it was not mentioned in the article how the constraints were handled. All these 12 problems had two objectives and the number of decision variables for all problems was not mentioned. The algorithm was also tested on a continuous steel casting and an electrocardiography problem with four variables and 2-3 objectives. The algorithm was compared with differential-based EA Robic and Filipic (2005) and an algorithm based on NSGA-II-ANN Deb and Nain (2007). The algorithm obtained better performance when measured with hypervolume in 10000 function evaluations.
An algorithm called K-RVEA was proposed in Chugh et al. (2016a) to solve computationally expensive problems with more than three objectives. In this algorithm, the meta-models were updated based on the need of convergence or diversity. Angle penalized distance (Cheng et al. 2016) and uncertainty information from the Kriging models with the help of reference vectors are used to select individuals in step 9 for updating the metamodels. In addition, extra individuals are removed from the archive in step 10 to further reduce the computation time.
The proposed algorithm was tested on DTLZ and WFG benchmark problems with 3-10 objectives and 10 decision variables. It was also compared with ParEGO, MOEA/D-EGO and SMS-EGO using IGD and hypervolume. In addition to benchmark problems, the algorithm was also tested on a free-radical polymerization problem (Mogilicharla et al. 2014) and also compared with the state-of-the-art algorithms. In the given number of function evaluations, the proposed algorithm performed better than other algorithms.
The same algorithm was also extended to handle constraints in Chugh et al. (2016b). Three different approaches are used to while training metamodels based on the feasibility of solutions. The proposed algorithm was tested on constrained version of DTLZ problems (Jain and Deb 2014) with 3-10 objectives and 10 decision variables. The authors found out that infeasible solutions are having a vital role on the performance of surrogates. Roy and Deb (2016) used a high-dimensional model representation (HDMR) by Sobol (2001) as a metamodel for each objective function. Kriging was further used within HDMR to approximate its component functions. The algorithm was proposed to handle problems with large number of decision variables. In each HDMR model, n component functions were approximated using Kriging; therefore, n × k (n and k represent the number of decision variables and objectives respectively) Kriging models were built and k HDMR models were built. After building the metamodels in step 5, NSGA-II was used to from steps 6-7. To update HDMR models, a prefixed number of individuals were re-evaluated using the original functions in step 9. These individuals were selected by doing clustering in the decision space and then added to the training archive in step 10. In addition, bounds of the decision variables were updated to limit the search space after metamodels were updated.
The algorithm was tested on 17 benchmark problems from ZDT, DTLZ and CEC09 suite ) biobjective problems with 15-30 decision variables. It was compared with NSGA-II and Kriging-based NSGA-II (Luo et al. 2014) using IGD and performed better than others in 500 function evaluations.

Non-Kriging-based algorithms
In this subsection, algorithms using metamodels other than Kriging such as neural network, support vector regression and polynomial regression are discussed. Syberfeldt et al. (2008) proposed a multiobjective parallel surrogate-assisted evolutionary algorithm (MOPSA-EA), where a feedforward neural network was used as a metamodel with a steady state EA. The metamodel is built for each objective function in step 5. A fixed number of offspring individuals is generated in step 6 using crossover and mutation and evaluated using the metamodel in step 7. The fitness values of offspring individuals are altered as per the fitness values of parents. The authors mentioned that this approach was motivated by fitness inheritance, where the fitness values of offspring depend on the fitness values of parents.
To get the altered fitness values for offspring individuals, the parents which are used for creating the offspring individual are evaluated with the metamodel and the error is calculated between true fitness values and the approximated fitness values (using the metamodel) for parents. For example, if the errors (for a biobjective optimization problem) between the true fitness values and the approximated fitness values for two parents are (a e , b e ) and (c e , d e ) and an offspring individual is generated using these parents having fitness values (e, g), then the new altered fitness values for offspring are The weight coefficients w 1 and w 2 are selected based on the influence of parent individuals on an offspring individual during crossover. An individual is selected in step 9 after getting new fitness values for offspring individuals and added to the archive in step 10 to update the metamodel. To select this individual, a nondominated sorting for offspring individuals and individuals in the archive is performed and nondominated individuals in both sets are identified. Let O R1 and P R1 denote nondominated individuals in offspring and in the archive, respectively. These individuals are combined and individuals in O R1 dominating P R1 are identified. Among these individuals, individual having the largest Euclidean distance to its closest individual in P R1 is selected and added to the archive. A nondominated sorting is performed again for individuals of the archive and the worst individual (having the worst fitness value and the smallest crowding distance) is removed from the archive. The size of the archive is fixed in this algorithm.
The MOPSA-EA was tested on the ZDT1-4 and ZDT6 benchmark problems (2 objectives and not explicit information of the number of decision variables) and on a manufacturing MOP (with two objectives and 11 decision variables). Generational distance (GD), inverted generational distance (IGD), spread and hypervolume metrics were used to compare the proposed algorithm with SMS-EMOA (Emmerich et al. 2005), MAES (Emmerich et al. 2006) and NSGA-II-ANN (Nain and Deb 2005). For the same number of function evaluations, MOPSA-EA performed better than other algorithms in Y , Ω and S metrics and in Δ, SMS-EMOA and MOPSA-EA performed equivalent. As the Pareto front was not known for the manufacturing problem, only the S metric was used as the performance criterion and MOPSA-EA gave better results in S than the other algorithms in the same number of function evaluations.
A quadratic polynomial approximation was used as a metamodel with the EA μ-MOGA (Coello and Pulido 2001) in Liu et al. (2008). In this algorithm, the bounds of decision variables are updated after every generation using a trust region algorithm. In step 5, the metamodel is built for each objective and constraint function. To update the metamodel in step 9, a fixed number of uniformly distributed individuals (P a ) from nondominated individuals is re-evaluated using the original functions. The nondominated individuals in the decision space after re-evaluation are stored in a set P e . The set P o = P a ∩ P e is determined which is then used to calculate a reliability index N (P o )/N (P a ), where N (P o ) and N (P a ) are the numbers of individuals in P o and the number of nondominated uniformly distributed individuals evaluated using the metamodel, respectively. This reliability index is used to update the bounds of the decision variables in the next generation with a trust region algorithm. The trust region radius is updated according to Alexandrov et al. (1998), and the algorithm terminates if the trust region radius is smaller than a predefined limit or after a fixed number of generations (step 8). A Latin hypercube design is used for sampling the decision variables with updated bounds. These individuals with the updated bounds are evaluated with the original functions which are used to update the metamodel.
Step 10 is not applicable in this algorithm as individuals after step 9 are not added to the archive.
The proposed algorithm was compared with μ-MOGA using two benchmark problems (two objectives and 2-3 decision variables) and one structural engineering problem (two objectives, three decision variables and one constraint). Using generational distance as a performance metric for the benchmark problems and visual comparison for the structural engineering problem, the quality of the obtained set of nondominated solutions within a fixed budget of function evaluations using μ-MOGA with metamodel was better.
Mitra and Majumder (2011) used a feedforward neural network as a metamodel with NSGA-II. The metamodel is built for each objective function in step 5 if a threshold based on a predicted tolerance is satisfied. This predicted tolerance, which is an indication of the accuracy of the metamodel, is calculated (however, calculation of the predicted tolerance is not detailed in the paper) after every generation. If this predicted tolerance is less than a user-specified tolerance, then the metamodel is used in the next generation. Otherwise, the original functions are used to evaluate individuals. The predicted tolerance is updated after every generation and again a decision is to be made either to use the original functions or the metamodel, and thus, an adaptive evolution control strategy is used. Individuals obtained after every generation are added to the archive (steps 9 and 10), and as a result,  the size of the archive grows with generations. In this algorithm, the metamodel is not updated after every generation but after every generation, it is checked whether the existing metamodel is sufficient enough to predict function values to the extent of accuracy required.
This algorithm was tested on an iron induration MOP with two objectives, 22 decision variables and three constraints.
To check the efficiency of the proposed algorithm, a graphical comparison was presented to the solutions from the algorithm and NSGA-II without using any metamodel. Similar nondominated individuals were obtained in 50% fewer function evaluations, and for the same number of function evaluations, better nondominated individuals were obtained.  extended the algorithm proposed in Liu et al. (2008). The main differences include the type of metamodel, stopping criterion in step 8 and selection of individuals for re-evaluation for updating the metamodel in step 9. The algorithm proposed utilizes a RBF as a metamodel. In addition to the stopping criteria posed in Liu et al. (2008) mentioned earlier, the algorithm also terminates if the bounds of the decision variables are equal to predefined limit and the reliability index is equal to one. In step 9, the individuals are selected using an inherited Latin hypercube design (ILHD) (Wang 2003) and a local-densifying strategy (to reduce the possibility of an ill-conditioned RBF matrix) for re-evaluation and subsequently updating the RBF.
The proposed algorithm was tested on eight different benchmark problems (two objectives and 1-5 decision variables) and a structural engineering problem (two objectives and five decision variables). The results were compared with μ-MOGA without any metamodel and the algorithm proposed in Liu et al. (2008). In case of benchmark problems, the proposed algorithm obtained better values for spread and convergence metrics (Deb et al. 2002) in fewer function evaluations. For the structural engineering problem, the proposed algorithm was not tested using μ-MOGA without any metamodel. Kourakos and Mantoglou (2013) proposed a multiobjective surrogate-assisted (MOSA) algorithm, where a modular neural network (Kourakos and Mantoglou 2009) was used as a metamodel with NSGA-II. A metamodel is built for each objective function in step 5. A fixed number (equal to the population size of NSGA-II) of better performing individuals (the authors did not mention any criterion for defining better performing individuals) is used as the population of NSGA-II and the nondominated individuals (A) after this step are determined. The nondominated individuals after evaluating the offspring individuals (generated in step 6) are compared with A to select one individual for re-evaluation (step 9). To do this, the individuals evaluated using the metamodel are clustered into two sets A 1 and A 2 as shown in Fig. 1. The first set (A 1 ) consists of the individuals which dominate at least one individual of A, while the second set (A 2 ) consists of the individuals that do not dominate nor are dominated. For each offspring in A 1 , the Euclidean distances between the offspring and the individuals in A are calculated. As shown in left part of Fig. 1, individual (o 1 ) with the largest distance is selected for re-evaluation using the original functions. In case A 1 is empty, an individual is selected from the set A 2 for re-evaluation. To select the individual, for each offspring individual in A 2 , the normalized perimeter for rectangles created between consecutive individuals in A (as shown in the right part of Fig. 1) is calculated. An individual located in the rectangle with the largest perimeter (o 2 ) is selected for re-evaluation using the original functions. In case both A 1 and A 2 are empty, offspring are generated again. However, the authors did not mention about how this algorithm can be used for more than two objectives. The selected individual is added to the archive to update the metamodel in step 10. The size of the archive is not fixed in this algorithm.
The algorithm was used to solve a coastal aquifer management optimization problem with two objectives and eight decision variables. The time required to evaluate the original functions was 26 h while the time to train the metamodels was 63 min. A graphical comparison was presented between MOSA and NSGA-II and the proposed algorithm gave similar results in fewer function evaluations.
In Herrera et al. (2014), support vector regression was used as a metamodel for approximation of each objective function with NSGA-II. The main focus in this algorithm is to use different basis functions for different kinds of variables such as discrete, continuous and categorical variables. It is to be noted that this algorithm does not use an ensemble of metamodels as only one prediction is obtained for each objective function by using different basis functions for the different kinds of variables. To convert categorical values into real number, dummy coding is used. In step 5, the metamodel is built for each objective function. Steps 9 and 10 are not applicable to this algorithm as the metamodel and the archive are not updated.
The proposed algorithm was tested on one real-world problem (two objectives and 10 decision variables) and compared with NSGA-II. Out of 10 variables, 5 were continuous and 5 were categorical variables. A visualization of nondominated individuals in the objective space was performed to compare the two algorithms. The authors mentioned that the proposed algorithm performed similar to NSGA-II in fewer function evaluations.
A support vector regression was used a metamodel in Pilát and Neruda (2014). This algorithm is known as HO-MOMA, where a metamodel is built for each objective function in step 5. After evaluating new individuals in step 7, a local search is used to optimize the fitness function obtained from the metamodel evaluations. This is done as follows. Firstly, several nondominated fronts are obtained with nondominated sorting. Then, in each front, sorting is performed according to the first objective function values in the ascending order. A reference point is then calculated for each individual in each front using these sorted values. After calculating the reference point, a fitness value for each individual is calculated. Letf = f 1 ,f 2 be the predicted objective vector for one individual and r = {r 1 , r 2 } the reference point for that individual. Fitness is then calculated as where d(r ,f ) is the Euclidean distance between r andf . This fitness is then optimized using CMA-ES (Hansen and Ostermeier 2001) as the local search algorithm. To update the metamodel in step 9, nondominated individuals obtained after the local search are added to the archive in step 10. The size of the archive is fixed in this algorithm, and extra indi-viduals are removed from it randomly. The authors mention the study for more than two objectives is a future research. HO-MOMA was tested on 14 benchmark problems with two objectives and 10-30 decision variables. The algorithm was compared with NSGA-II and ASM-MOMA (Pilát and Neruda 2011a) with hypervolume as the performance criterion. The proposed algorithm performed better than the others in nine out of 12 problems.
An extreme learning-based algorithm called ELMOEA/D-DE was proposed in Pavelski et al. (2014Pavelski et al. ( , 2016, which was inspired by MOEA/D-RBF. The main focus in this algorithm is to use the metamodel for higher-dimensional problems in the decision space. Extreme learning is a single-layer feedforward neural network proposed in Huang et al. (2004). In step 5, the metamodel is built for each objective function. To update the metamodel in step 9, the same procedure is used as in MOEA/D-RBF. In addition, a minimum distance (in the decision space) is maintained between the individuals to be selected for re-evaluation. After adding these individuals in the archive, inferior solutions (in terms of scalarized single objective problem) are removed. Otherwise, the closest individual in the objective space is replaced by the individual obtained in step 9. This is done to ensure that every new solution is added to the archive for updating the metamodel. Pavelski et al. (2016), the algorithm was tested on ZDT, DTLZ and WFG problems with 2-5 objectives and 5-60 decision variables and on an airfoil shape optimization problem with two objectives and 12 decision variables. The algorithm was also compared using addictive epsilon and R 2 (Zitzler and Thiele 1998) as the performance indicators. For the given number of function evaluations, the proposed algorithm obtained better performance.

ELMOEA/D-DE was tested in Pavelski et al. (2014) on ZDT problems with 10-60 decision variables and compared with MOEA/D-DE and MOEA/D-RBF. Later in
In Akhtar and Shoemaker (2015), an algorithm called gap-optimized multiobjective optimization using response surfaces (GOMORS) was proposed, where radial basis functions were used to approximate the objective functions. After building the metamodels in step 5 for each objective function, an EA proposed in Vrgut and Robinson (2007) was used for optimization from steps 6-7. After step 7, an another optimization problem was solved using the same EA by reducing the bounds of the decision variables. This problem was referred as gap optimization problem in the article. In step 9 for updating the metamodels, four criteria were used and one individual corresponding to these criteria was selected and added to the training archive in step 10. Four criteria were based on hypervolume, distance to the individuals in the decision space, distance to the individuals in the objective space and hypervolume in the gap optimization problem. A maximum number of function evaluations was used as the termination criterion.
The algorithm was tested on 11 benchmark problems with 8-24 decision variables and two objectives. In addition, a groundwater remediation problem with 6-24 variables and two objectives was also solved. The algorithm was compared with NSGA-II and ParEGO using hypervolume in 200-400 function evaluations. The proposed algorithms obtained better performance than others in the given number of function evaluations.
In Palar et al. (2015), an algorithm called surrogateassisted local search memetic algorithm (SS-MOMA) was proposed, where RBF as a metamodel was build on the single objective problem after converting multiobjective optimization problem using a scalarization function. Two common ways of scalarization, i.e., Tchebycheff and weighted sum, were tested, where weights were generated randomly and the reference point in Tchebycheff scalarization was the current individual objective function values. After generating the offspring population in step 6, metamodels were built locally for each individual. For instance, if the offspring population of size 100 was generated, one metamodel for each individual (i.e., 100 in total) was built. As multiobjective optimization problem was converted into single objective using the scalarizing function, SQP algorithm was used to obtain the solutions. All solutions thus obtained were re-evaluated with the original functions and added to the training archive. Afterward, these solutions were combined with the parent individuals and nondominated sorting (Deb et al. 2002) was performed to obtain the population for the next generation. In this way, after every generation metamodels were updated with a number equal to the population size used.
The algorithm was tested on three benchmark problems with 15 variables and two objectives. The algorithm was not compared with any other algorithm. However, in 1606 function evaluations, Tchebycheff scalarization performed better than weighted sum in terms of generational distance (Deb 2001) and one diversity metric used. Datta and Regis (2016) used RBFs as metamodels for each objective and constraint function in step 5. A (μ + λ) evolution strategy mutation operator was used to generate new individuals in step 6 which were then evaluated using the metamodels in step 7. After using metamodels for each objective and constraint functions, feasible solutions were found and a nondominated sorting was performed on these feasible solutions. The best individuals, i.e., individuals in the first front were re-evaluated with the original functions to update the metamodels in step 9 and added to the training archive in step 10. Therefore, the maximum number of individuals to be updated was λ. In addition, initial training of metamodels in step 5 was performed without considering any feasibility or infeasibility of solutions. However, the authors clearly mentioned that the algorithm is not expected to work well on problems when the feasible region is empty.
The algorithm was tested on 15 benchmark problems with 2-15 decision variables, 2-5 objective functions and 2-13 constraint functions. In addition it was also tested on a manufacturing and robotics problem with 3-7 decision variables, 2-5 objectives and 2-8 constraints. Hypervolume was used to compare the performance of the algorithm against constrained version of (μ + λ) evolution strategy Bäck (1996) and NSGA-II. The algorithm performed better in 500-5000 function evaluations.
A RBF was used as a metamodel for each objective and constraint function in Regis (2016). It was assumed that at least one feasible solution is available to train the metamodels in step 5. To create new individuals in step 6, two different approaches were tested. In the first one, uniform random individuals were generated over the search space and in the second one, individuals were generated by adding Gaussian perturbation centered at the nondominated individual that has the most isolated objective function values. The most isolated individuals are identified by measuring the distance among nondominated individuals. Metamodels were then used to approximate the objective and constraint function values in step 7. Afterward, nondominated solutions with minimum constraint violations were found. Among these individuals, one individual was selected for re-evaluations in step 9 to update the metamodels. One individual was selected by the weighted sum (with equal wights) of two criteria. One criterion was the distance of solutions (in the decision space) obtained in step 7 from the individuals in the training archive and the second criterion was the distance from the nondominated individuals from the last generation in the objective space. Selected individual was then re-evaluated and added to the archive in step 10.
The algorithm was tested on 28 benchmark problems with 2-5 objectives, 2-15 decision variables and 1-11 constraint functions. the algorithm was compared with its two different versions (based on the generation of individuals in step 7), NSGA-II and DirectMultiSearch(DMS) (Custodio et al. 2011) using hypervolume as the performance indicator. Overall, the version where individuals were generated using Gaussian distribution performed better in the given number of function evaluations.

Algorithms based on multiple metamodels
In this subsection, algorithms using multiple metamodels are discussed. These metamodels are used independently to predict objective and/or constraint functions.
In Husain and Kim (2010), three independent case studies were performed, where three different metamodels (polynomial approximation, Kriging and radial basis function) were used with NSGA-II. The metamodel is built for each objective function in step 5. The nondominated individuals obtained after this step are improved using a local search algorithm with sequential quadratic programming. To perform local search, a variant of -constraint algorithm (Miettinen 1999) is used, i.e., one of the objectives is optimized and other objectives are converted to equality constraints. The improved individuals are combined with individuals from step 7 and dominated and duplicated individuals are eliminated. To update the metamodel in step 9, a fixed number of individuals is selected from the remaining individuals using K-means clustering in the objective space for re-evaluation using the original functions. These individuals are then added to the archive in step 10. The size of the archive is not fixed in this algorithm.
This algorithm was used to solve a heat sink MOP with two objectives and three decision variables. Three different studies were performed to compare the results while using different metamodels. In the first one, nondominated individuals were identified while using each metamodel involving steps 1-8 (i.e., without updating the metamodel) of the function approximation framework. Five representative individuals among nondominated individuals were selected using K-means clustering and re-evaluated with the original functions. The proposed algorithm with Kriging gave the least error in objective function values for these five individuals. In the second study, five representative individuals obtained using each metamodel were re-evaluated with other metamodels. For example, individuals obtained while using polynomial approximation were re-evaluated with Kriging and radial basis function. Individuals obtained while using Kriging when re-evaluated with other metamodels gave the least error in objective function values. In the third one, nondominated individuals were obtained utilizing steps 1-10 (i.e., by updating the metamodel in steps 9 and 10) and results were compared graphically in the objective space with the first study. While using Kriging and radial basis function, better results were obtained when compared to results of the first study and while using polynomial approximation, similar results were obtained. However, the proposed algorithm was not tested using NSGA-II without any metamodel and the efficiency of the algorithm was not mentioned in terms of computation time or the number of function evaluations.
Arias-Montano et al. (2012) used five RBFs with different basis functions as metamodels along with the EA MODE-LD+SS (Arias-Montano et al. 2010). However, these metamodels were used independently for each objective function but individuals from each metamodel evaluation were used to update all metamodels as discussed next. A metamodel is built for each objective function in step 5. In step 9, to select the individuals for updating the metamodel, one individual from each metamodel evaluation is selected for re-evaluation using the original functions and added to the archive in step 10. To select one individual, a set of uniformly distributed weight vectors (λ 1 , λ 2 , . . . , λ N ) (where N is the population size of EA) is defined. Next, from each metamodel evaluation, the individual is selected that minimizes the Chebyshev scalarizing function (Miettinen 1999) given . The values of objective functions obtained after step 7 and the minimum value of the objective function in the population of EA at the current generation are represented byf and f * , respectively. The authors mentioned that this updating criterion can balance the accuracy of the metamodel and the diversity among individuals. The size of the archive is fixed in this algorithm and extra individuals are removed from it based on their rank.
This algorithm was tested on five different aerodynamic shape optimization problems with 2-3 objectives and 12 decision variables. Hypervolume was used as the comparison criterion for the results of MODE-LD+SS with and without metamodels. The proposed algorithm gave a better hypervolume in fewer function evaluations for all five optimization problems.
In Palar et al. (2016), SS-MOMA described in Sect. 3.3.2 was tested with different metamodels for different problems, i.e., RBF with linear kernel function for ZDT problems and Kriging for UF problems because of different function landscape of the problems. Using a particular technique for a problem with an assumption that function landscape is known a priori raises the question of applicability of the algorithm to black box problems. In addition, one modification of achievement scalarizing function used, where reference point was replaced by the upper bounds of the objective function values at the current generation was also tested. An another local search algorithm called random mutation hill climber (Aarts and Lenstra 2003) was also tested after building the metamodels for each objective function. To handle constraints, metamodels were first built for each objective function and after scalarization, a constrained SQP algorithm was used to solve the single objective optimization problem.
The algorithm was tested on seven biobjective benchmark problems with 8-15 decision variables and a biobjective airfoil optimization problem with 16 variables and one constraint. Different versions of the algorithm with different scalarizing functions were compared with each other also with NSGA-II using IGD. Overall a low number of function evaluations were used to obtain a given IGD value when Tchebycheff function with reference point as the individual objective function values was used.

Algorithms based on ensemble of metamodels
In this subsection, we discuss algorithms using ensemble of metamodels. As defined in Sect. 2, either the weights are given to predicted output of the metamodel or a metamodel is selected based on its accuracy. In Lim and Jin (2010), an ensemble of metamodels (Kriging, polynomial regression and radial basis function) was used with the proposed generalized surrogate multiobjective memetic algorithm (GS-MOMA). In this algorithm, the offspring population is generated in step 6 before building the metamodel. A fixed number of individuals, equal to n + (n + 1)(n + 2)/2, (where n is the number of decision variables) is selected from the archive to build an ensemble of metamodels and a separate polynomial regression metamodel in step 5. The individuals are selected in the decision space using the Euclidean distance between the individuals in the offspring population and the individuals in the archive. The individuals which are closer to offspring individuals are used to build the metamodels.
The ensemble fitness values of offspring individuals are calculated as F ens (x) = m j=1 ω j F j (x), where ω j is the weight coefficient for fitness values F j (x) of the jth metamodel (step 7). The weight coefficient is assigned based on the accuracy of the metamodels which is calculated using root mean square error. A local search algorithm (sequential quadratic programming) is also used for single objective optimization of F j (x) to improve the individuals evaluated with the metamodel. The best found individuals after this step and the local search algorithm are combined with the parent population, and a selection mechanism is used to select individuals for the next generation. To update the metamodels, an offspring population is generated and individuals are selected from the initial archive based on the mechanism used to build the metamodel. Steps 9 and 10 are not applicable in this algorithm and the size of the archive is fixed.
This algorithm was tested on six benchmark problems (labeled as MF problems in the paper) with 2-3 objectives and with 10-50 decision variables. Three performance criteria, namely generalized distance (Veldhuizen and Lamont 1998), maximum spread (Zitzler 1999) and hypervolume ratio , were used for comparison between GS-MOMA and NSGA-II without a metamodel. GS-MOMA performed better in all performance criteria for all problems in the same number of function evaluations.
A surrogate-assisted simulated annealing (SASA) algorithm was proposed was proposed in Singh et al. (2010), where an ensemble of metamodels (quadratic polynomial and radial basis function) was used with constrained Pareto simulated annealing (C-PSA) as an EA. In this algorithm, one of two different metamodels is selected to evaluate individuals based on their accuracy which is calculated using the root mean square error. A fixed number of recently evaluated individuals in stage 1 is used to create the metamodels in step 5. In what follows, for each objective function, either one of the metamodels or the original functions are used to evaluate the offspring individuals in step 7, and thus, an adaptive evolution control strategy is used. In other words, if none of the metamodels are accurate enough (accuracy is compared with a predefined parameter), the original functions are used. To select the individuals for updating the metamodel in step 9, nondominated individuals after step 7 (in case a metamodel is used) are re-evaluated with the original functions and added to the archive in step 10. Dominated individuals are removed from the archive and the remaining individuals are used to update the metamodels. The size of the archive is fixed here and clustering is used via a linkage algorithm (Jain and Dubes 1998) to remove extra individuals.
The SASA algorithm was tested on eight benchmark problems (Deb 2001) with two objectives, 10 decision variables and 1-2 constraints. The hypervolume and displacement metric (Bandyopadhyay et al. 2008) were used to compare SASA with NSGA-II and the infeasibility driven evolutionary algorithm (IDEA) (Ray et al. 2009). For the same number of function evaluations, SASA performed better than the other algorithms for seven problems and for one problem, IDEA performed the best in both performance criteria.
In Martinez and Coello (2013), a similar algorithm to MOEA/D-EGO was proposed, where radial basis function was used as the metamodel. This algorithm is known as MOEA/D-RBF, where the metamodel is built for each objective function instead of scalarized objective function. An ensemble of metamodels (RBF with three different basis functions) is used for each objective function in step 5. A weighted sum of predictions from each metamodel is used in the ensemble of metamodels and the weights are decided based on the prediction error of the metamodel. To update the metamodel in step 9, the following procedure is adopted.
In MOEA/D let, W = w 1 , . . . , w N are a uniformly distributed set of weight vectors and R is the number of points to be re-evaluated in step 9. An another set of uniformly distributed weight vectors W R = w 1R , . . . , w N R is defined such that size of W R < size of W as shown in Fig. 2. For each w i R ∈ W R , a neighborhood is defined B R (w i R ) = w 1 , . . . , w Na , such that w 1 , . . . , w Na ∈ W are the N a closest weight vectors from W to w i R . After this, from each neighborhood, an individual which minimize the single objective optimization problem used in MOEA/D is selected for re-evaluation and added to the archive in step 10. In this algorithm, a penalty boundary intersection approach (Zhang and Li 2007) is used for decomposition of the MOP into single objective optimization problems. The size of the archive is fixed in this algorithm and the same procedure mentioned above is applied to eliminate extra individuals from the archive.
MOEA/D-RBF was tested on five benchmark (two objectives and 8-30 decision variables) and one real-world problem (two objectives and 11 decision variables). It was compared with MOEA/D (for both benchmark and realworld problems) and MOEA/D-EGO (only for benchmark problems) with hypervolume as the comparison criterion. The proposed algorithm obtained better performance in five out of six problems in the same number of function evaluations.

Comparison of function approximation-based algorithms
In what follows, a comparison of function approximationbased algorithms discussed so far is presented in Fig. 3. This figure represents the number of papers with respect to the metamodel (upper part of the figure), EA (lower part of the figure) and evolution control strategy used. The references for these three criteria are mentioned inside the bars of the figure. One should note that some algorithms mentioned in this section were not tested on computationally expensive MOPs. These algorithms are still cited here as they have been proposed for considering the computational burden in any MOP. The number of papers which used benchmark and realworld problems is also mentioned in Fig. 3. Seven algorithms (Lim and Jin 2010;Pavelski et al. 2014;Pilát and Neruda 2014;Singh et al. 2010;Syberfeldt et al. 2008;Zhang et al. 2010) were tested on benchmark problems, seven algorithms (Arias-Montano et al. 2012;Herrera et al. 2014;Husain and Kim 2010;Jang et al. 2009;Kourakos and Mantoglou 2009;Mitra and Majumder 2011;Singh et al. 2014) on real-world problems and seven algorithms Li et al. , 2008Liu et al. 2008;Liu and Collette 2014;Martinez and Coello 2013;Zhu et al. 2013;Chugh et al. 2016a) on both. The efficiency of the different algorithms in terms of computation time or number of function evaluations reduced is very important, especially in the case of real-world problems. In some cases, the authors mentioned that it was practically difficult to do optimization without approximations due to high computation time and the algorithm was not compared with an EA without any metamodel.
As far as selection of a metamodel is concerned, as mentioned in Sect. 3.1, there is no general rule or correlation between a metamodel for approximation and a particular problem to be solved. As shown in Fig. 3, a single metamodel is used more than multiple metamodels and ensemble of metamodels. In algorithms using a single metamodel, Kriging (Jang et al. 2009;Li et al. 2008Liu and Collette 2014;Singh et al. 2014;Zhang et al. 2010;Zhu et al. 2013) has been used more than other metamodels such as neural networks Kourakos and Mantoglou 2013;Mitra and Majumder 2011;Pavelski et al. 2014;Syberfeldt et al. 2008), support vector regression (Herrera et al. 2014;Pilát and Neruda 2014) and polynomial regression (Liu et al. 2008). Zhu et al. (2013) used Kriging as a metamodel and mentioned that within different metamodels studies to date, perhaps the most common type of metamodel used is the Kriging model. In , RBFs were used and the authors mentioned that RBF was a kind of approximation having a very good approximation accuracy. Arias-Montano et al. (2012) mentioned that RBFs were very powerful functions to represent complex fitness landscapes. In addition, the authors mentioned that Kriging had a strong mathematical basis, and is probably one of the most powerful interpolation algorithms currently available.
Managing the metamodels or evolution control is also very important as it affects the performance of the metamodel used. As shown in Fig. 3, the fixed evolution control strategy was used more than the adaptive evolution control strategy. As mentioned in Sect. 2, using an adaptive evolution control strategy depends on the accuracy of the metamodel used and using a fixed evolution control strategy implies that either the metamodel is accurate enough for approximation or the metamodel accuracy is not important or not checked. There are only five algorithms (Liu and Collette 2014;Mitra and Majumder 2011;Mlakar et al. 2015;Singh et al. 2010;Zhu et al. 2013), where adaptive evolution control was used. For instance, Liu and Collette (2014) compared the uncertainties of approximated values with a predefined value and if that uncertainty was acceptable, then only Kriging models were further used otherwise original function were used. A similar strategy was followed by Mitra and Majumder (2011), where accuracy of neural networks was measured after every generation and compared with a predefined value. In Mlakar et al. (2015), uncertainties of offspring were compared with the uncertainties of parents and based on the comparison, the decision was made to use Kriging models or original functions. In Singh et al. (2010) and Zhu et al. (2013) used the approximation accuracy of metamodels for using them in subsequent evaluations. Moreover, the training time for different metamodels was not considered in many papers though training time can be substantially high in some cases, particularly, when the metamodel approximation accuracy is important. Among EAs, dominance-based algorithms are more widely used than other algorithms, e.g., decompositionbased or indicator-based. Moreover, in dominance-based EAs, NSGA-II was used more often than other EAs.

Problem and fitness approximation
In addition to function approximation, problem and fitness approximations are also used in the literature to reduce the computation time in multiobjective optimization problems. In problem approximation, the original problem is replaced by a simplified problem which is faster to solve. The main goal in problem approximation is to reduce the computational complexity of the problem. For example, in computational fluid dynamics or structural analysis, the governing equations can be modified to reduce the computation time. As an example, replacing three-dimensional Navier-Stokes equations by two-dimensional Euler equations Lattarulo et al. (2013) reduces the computational complexity of the problem. Euler equations were used instead of Navier-Stokes equations in Lattarulo et al. (2013) and Lee et al. (2008) to solve aerodynamic shape optimization problems. A similar approach was followed in Mengistu and Ghaly (2008); Oyama et al. (2009), where two-dimensional Navier-Stokes equations were used for solving aerodynamic shape optimization problems.
Recently, data-driven optimization algorithms Wang et al. (2016) are also proposed to solve problems, where an analytical forms or simulations models for the objective functions are not available and some data are available obtained through some physical experiments. Therefore, to obtain Pareto optimal solutions, one has to rely upon the data available. In Wang et al. Wang et al. (2016), a MOP was formulated using the data available from a trauma system. In addition, the authors found out that accuracy and size of the data used in optimization are conflicting.
In fitness approximation, we have algorithms which use metamodels for approximating some element instead of objective functions. For instance, approximating a rank, classifications of individuals into nondominated and dominated sets and distance to the nondominated individuals. Moreover, we consider the papers in the categories of fitness inheritance and fitness imitation (defined in Sect. 2) in this section.
Fitness inheritance was originally proposed in Smith et al. (1995) to improve the performance of genetic algorithms. Two types of fitness inheritance were proposed in that paper. In the first one, known as averaged inheritance, the fitness values of the offspring were calculated by taking the average of the fitness values of the parents. In the second one, known as proportional inheritance, weighted average of the fitness values of the parents were assigned to offspring and weights were assigned according to common elements between offspring and parents. This algorithm was tested on OneMax problem Ackley (1987), and the algorithm with fitness inheritance gave better performance (a graphical comparison was performed) than without inheritance in fewer function evaluations.
In what follows, some papers are cited here which use the concept of fitness inheritance. One should note that all these algorithms were not tested on computationally expensive MOPs. Chen et al. (2002) proposed analytical algorithms for computing convergence time and population size for using the fitness inheritance in MOPs. Ducheyne et al. (2003) used fitness inheritance with a binary genetic algorithm (GA) for the ZDT benchmark problems. The authors found a similar performance of the binary GA with and without using the fitness inheritance for problems having convex and continuous Pareto fronts. In case of nonconvex and discontinuous Pareto fronts, the binary GA without inheritance performed better and the authors mentioned that fitness inheritance can only be applied to problems having convex and continuous Pareto fronts. In contrast, Reyes-Sierra and Coello (2005) found out that the fitness inheritance can also be applied to problems having nonconvex and discontinuous Pareto fronts. Loshchilov et al. (2009) used a metamodel to predict the class for the individuals. This class was the indication that individuals were either nondominated or dominated. This algorithm is known as Pareto SVM, where one class SVM is used for the classification of individuals in the decision space into nondominated and dominated. In addition, support vector regression is used to predict those individuals (after classification in the decision space) to a target value in the objective space. Aggregate surrogate metamodel terminology is used in the paper as two metamodels are used, one in the decision space and another in the objective space. Instead of binary classification in the decision space, all nondominated individuals mapped to the single value ρ with tolerance and all dominated individuals are mapped to ] − ∞, ρ − [. In this way, the individuals belonging to [ρ + , +∞[ are in the unexplored region. In step 9, an updating criterion for the metamodel is inspired from EA PESA-II (Corne et al. 2001). Firstly, individuals obtained after step 7 are added to the archive and duplicate individuals are removed. Then, the objective space is partitioned into a fixed number of equal sized hyperboxes and one individual or nondominated individual (if the box contains nondominated) is selected from each box randomly. These individuals are re-evaluated using the original functions. The size of the archive is fixed in this algorithm. In addition, to generate offspring in step 6, firstly, a fixed number of individuals (more than the population size) is generated using the variation operators. Then a fixed number of individuals is selected based on their distance to the current nondominated individuals in the decision space.
The proposed algorithm was tested on eight benchmark problems with two objectives and 10-30 decision variables. Hypervolume was used as the performance criterion for comparison of Pareto SVM with (μ + λ) − S−NSGA-II (Emmerich et al. 2005) and μ × (1 + λ)−MO-CMA-ES (Igel et al. 2007). Pareto SVM obtained similar performance in fewer function evaluations.
In Loshchilov et al. (2010), Pareto SVM was extended to a rank-based aggregate surrogate model. In this algorithm, instead of two classes, various classes were obtained in the decision space based on the rank of individuals. In this way, nondominated individuals at the current generation are not constrained to the bounds. A pairwise dominance relation was used to classify individuals. For example, two possibilities exist, if an individual dominates another individual or is dominated by it. The case for nondominated individuals is mentioned as the future work. Other details are the same as in Pareto SVM.
The proposed algorithm was tested on the same benchmark problems. It was compared with (μ+λ)−S−NSGA-II, μ × (1 + λ)−MO-CMA-ES and Pareto SVM with hypervolume as the performance criterion. In comparison with Pareto SVM, it performed similar and in comparison with other two algorithms, it performed better in fewer function evaluations.
Instead of predicting objective functions, Pilát and Neruda (2011a) used a metamodel to predict the distance to the current nondominated individuals in the decision space. Three metamodels, linear regression, support vector regression and multilayer perceptron, were used independently. This algorithm is known as ASM-MOMA, where the metamodel is built for the distance to the nondominated individuals in step 5. In addition, a local search algorithm is used to improve the individuals obtained after step 7. To update the metamodel in step 9, nondominated individuals at the current generation are added to the archive in step 10. The size of the archive is fixed in this algorithm, and extra individuals are removed from it randomly.
ASM-MOMA was tested on four benchmark problems (two objectives and 15 decision variables) and compared with NSGA-II and IBEA (Zitzler and Kunzli 2004) with hypervolume as the performance criterion. The proposed algorithm obtained similar performance in fewer function evaluations than other algorithms while using linear regression as the metamodel.
In Pilát and Neruda (2011b), ASM-MOMA was extended to high-dimensional problems (in the objective space) with two different additional elements. Instead of using one global metamodel, multiple local metamodels were used. This algorithm is known as LAMM-MMA, where the training data in the decision space are partitioned into different sets based on the distance of individuals to the current nondominated individuals. Different from ASM-MOMA, the algorithm allocates weights to samples based on the distance to nondominated samples in the decision space. Other details are the same as in ASM-MOMA. This algorithm was tested on four benchmark problems with 5-15 objectives and 20 variables and compared with IBEA and ASM-MOMA with hypervolume as the performance criterion. LAMM-MMA performed better than IBEA and similar to ASM-MOMA in the same number of function evaluations.
A metamodel was used in Azzouz et al. (2014) to predict the contribution to the hypervolume instead of the objective functions. This algorithm is known as NN-SS-IBEA, where a metamodel (neural network) is built for the contribution to the hypervolume of individuals in step 5. After evaluating the new individuals in step 7, one individual having the maximum contribution to the hypervolume is re-evaluated in step 9 and added to the archive in step 10. The size of the archive is fixed in this algorithm, and extra individuals are eliminated from the archive based on their hypervolume contribution.
The proposed algorithm was tested on 12 benchmark (2-3 objectives and 8-30 decision variables) and one realworld problem (two objectives and 11 decision variables). The algorithm was compared with MOEA/D-RBF (Martinez and Coello 2013) and IBEA (Zitzler and Kunzli 2004). The hypervolume was used as the comparison criterion. In case of benchmark problems, NN-SS-IBEA performed better in eight out of 12 problems with the same number of function evaluations. In case of the real-world problem, the proposed algorithm performed better than others for a given number of function evaluations.
A similar algorithm to Loshchilov et al. (2010) was proposed in Bandaru et al. (2014, where nondominated individuals were also used during classification. In this algorithm, a metamodel is built for three classes while doing pairwise dominance comparison. Here, three possibilities exist, one solution dominates another, one solution is dominated by another, or both solutions are nondominated. Ten different classification algorithms were used independently. The metamodel is updated in step 9 after a fixed number of generations θ 1 θ 2 rank 1 individuals rank 2 individuals rank 3 individuals Fig. 4 Assignment of ranks to individuals with the previously evaluated individuals in step 7. The information on the size of the archive is not mentioned. The proposed algorithm was tested on three benchmark problems (two objectives and 5-20 decision variables). However, this algorithm was not compared with any other algorithm, but the results of ten classification algorithms were compared using the training time and the accuracy as the comparison criteria. The authors mentioned that SVM, classification trees, k-neural network and quadratic discriminant analysis were the most preferred among other metamodels as these metamodels found the knee region in the objective space. Seah et al. (2012) used a metamodel to predict the rank of individuals in the objective space. In this algorithm, the metamodel is used to create boundaries between fronts in the objective space. As shown in Fig. 4, θ 1 and θ 2 represent the two boundaries created using individuals of three different fronts. The ranks of offspring individuals are then predicted with an indicator function mentioned in the paper with boundaries created (steps 5-7). For example, individuals that lie below θ 1 are of rank 1, individuals between θ 1 and θ 2 are of rank 2 and individuals above θ 2 are of rank 3. The first-rank individuals are re-evaluated using the original functions in step 9 and added to the archive in step 10. The size of the archive is not fixed in this algorithm.
This algorithm was tested on 19 benchmark problems (the authors did not mention explicitly the number of objectives and decision variables) and compared with NSGA-II, SPEA2 (Zitzler et al. 2002) and MOEA/D (Zhang and Li 2007). Three performance criteria [generational distance, inverted generational distance and hypervolume (Durillo et al. 2010)] were used to compare the proposed algorithm and EAs for the same number of function evaluations. Out of 19 problems, the proposed algorithm performed better than the other algorithms in 16, 14 and 12 problems in generational distance, inverted generational distance and hypervolume, respectively. In the rest of the problems, MOEA/D and NSGA-II performed better than the proposed algorithm. Toscano and Deb (2016) compared four metamodeling techniques and three different scalarizing functions using two different approaches of approximation. Therefore, for each approximation, 12 different experiments were performed on different problems. A data set of prefixed size was used for training and validation. The error after the validation was used to compare different studies. As no algorithm was proposed in the current study, other steps of the function approximation were not applied. Four metamodels were based on Kernel ridge regression (KRR), Kriging, generalized linear models and tree-based regression. Three different scalarizing functions tested were based on Tchebycheff (TCH), boundary intersection (PBI) and weighted sum (WS). As mentioned, two different approaches of approximation were used; in the first one, the metamodel was used for the scalarizing function, and in the second one, the metamodel was used to approximate the rank of individuals.
In the study, experiments were conducted on DTLZ suite with 4-10 objectives and 8-19 decision variables. In both types of approximation used, Kriging performed the best. In comparing different scalarizing functions, TCH performed the best in the first type of approximation and WS in the second type of approximation. Developing a algorithm was considered as a future work.

Discussion
In this section, we first summarize promising elements that we found in the algorithms proposed in the literature. Furthermore, we also discuss ways which can be used to design enhanced EAs to handle computationally expensive MOPs. Finally, we discuss the main issues we have observed in the algorithms in the literature with respect to using an approximation in an EA and the numerical settings used to test their efficacy.

Promising elements in the literature
We discuss here promising elements with respect to the choice of the metamodel to be used, building/updating the metamodels and using different types of approximations together to reduce the computation time. Furthermore, we discuss the potential of using hybrid algorithms with function approximation-based algorithms to enhance the rate of convergence of EAs near Pareto optimal solutions. Also, one can use these elements while designing an algorithm for using an approximation with EA to reduce the computation time.
Often when metamodels are used, the type of metamodel to be used is a random choice. It is often difficult to know a priori the best metamodel to be used in the solution process. An ensemble of metamodels (see, e.g., Lim and Jin 2010;Singh et al. 2010) could be an easy fix, where the algo-rithm is not limited to one single metamodel. Lim and Jin (2010) assigned weights to the fitness values predicted by each metamodel and the weighted sum is then used to get the final fitness value. In Singh et al. (2010), one metamodel having the highest accuracy is used for evaluating individuals of the EA population. Moreover, Arias-Montano et al. (2012) present an improved way of using multiple metamodels that focus on different parts of the Pareto front. These metamodels are updated using some solutions that were more and less accurate. Selection of these solutions enables the algorithm proposed by Arias-Montano et al. (2012) to converge faster and explore the search space for promising candidate Pareto optimal solutions.
One more promising element observed in some of the algorithms in the literature is to use a combination of different approximations. In Syberfeldt et al. (2008), fitness inheritance (a type of fitness approximation) is used with function approximation to reduce the computation time. Problem approximation is used with function approximation in Zhu et al. (2013) to decrease the computational complexity of the problem in question and the numbers of function evaluations. Moreover, the use of multilevel approaches (Giannakoglou and Kampolis 2010;Kampolis and Giannakoglou 2008) which depend on the problem to be solved needs more attention from both researchers and practitioners. In these approaches, either different solvers were used to solve governing equations (e.g., Navier-Stokes equations) or different problems were solved at different levels. For more details about multilevel approaches, see Giannakoglou and Kampolis (2010).
In addition, metamodels have been used in the literature to predict the quality of solutions instead of the objective functions. For instance, Azzouz et al. (2014) used a metamodel to predict the hypervolume contribution of the individuals. In Pilát and Neruda (2011a, b) used a metamodel to predict the distance from the current nondominated individuals. Moreover, metamodels were used by Bandaru et al. (2014); Loshchilov et al. (2009Loshchilov et al. ( , 2010 to classify the individuals into different classes based on their dominance relations. In Seah et al. (2012), a metamodel was used to predict rank of individuals by creating boundaries between individuals in the objective space. This way of using a metamodel other than function approximation could be affective as instead of actual objective functions, quality of individuals (e.g., using hypervolume or rank) is predicted.
Hybrid EAs with a combination of global and local searches are used to enhance the rate of convergence toward the Pareto front (Ishibuchi et al. 2008(Ishibuchi et al. , 2010Sindhya et al. 2011Sindhya et al. , 2013. However, such algorithms are not commonly used together with function approximation to handle computationally expensive MOPs. In the literature, e.g., in Husain and Kim (2010), Lim and Jin (2010) and Turco (2011), a local search is used with function approximation. This way of using local search algorithm with function approximation can inherit advantages of both hybrid algorithms and function approximation-based algorithms, i.e., fast convergence and reducing the numbers of computationally expensive function evaluations.

Performance metrics and benchmark problems
It is also important to point out that different algorithms use different performance metrics for measuring their efficiency compared to other algorithms. These metrics, e.g., generational distance (GD), inverted generational distance (IGD) and hypervolume (or S metric) are widely used in the literature. For definitions, see Coello et al. (2007) and Deb (2001). In Table 2, we mention the performance metrics used in different articles. As can be seen in the table, some articles do not use any performance metric which is notated by no preference measure and some articles compare solutions obtained with different algorithms visually, e.g., by plotting a scatter plot of nondominated solutions.
Moreover, different algorithms have been tested on different benchmark problems as also shown in Table 2. We mention the name of the benchmark problem suite, e.g., ZDT, DTLZ and WFG in the table. For more details about these problems, see Coello et al. (2007). From the names of the suites, one can observe the trend of using different benchmark problems in the field of surrogate-assisted multiobjective optimization.

Guidelines for selecting an algorithm
For practitioners and engineers, it is very important to select an appropriate algorithm for solving a real-world problem. In the literature, different algorithms use different EAs, surrogate techniques and other relevant parameters. Therefore, it is challenging to generalize the efficiency of an algorithm or choose an algorithm. However, one can consider the following points to select an algorithm: 1. Dimensions in both objective and decision spaces and 2. Budget, e.g., maximum number of expensive function evaluations available.
The first point in selecting an algorithm is the dimensions of the problem to be solved in both the objective and decision spaces. One can observe from Table 2 that some algorithms, e.g., (Arias-Montano et al. 2012;Azzouz et al. 2014;Liu and Collette 2014;Veldhuizen et al. 1999) were tested with up to three objectives, few (i.e., Chugh et al. 2016a;Ponweiser et al. 2008a;Pilát and Neruda 2011b) were tested with more than three and most of the algorithms were tested only on biobjective optimization problems. A similar pattern can be observed for decision variables: only two algorithms, i.e., Lim and Jin 2010;Pavelski et al. 2014, were tested on problems with more than 30 variables. Therefore, references provided in the table give options for selecting an algorithm based on the dimensions of the problem to be solved.
The second key point in selecting an algorithm should be the capability of the algorithm to obtain solutions in a given limited budget. In many real-world problems, the number of expensive function evaluations is limited, e.g., because of a time limit and it may be infeasible to do initial experimental runs to select an algorithm. As shown in Table 1, some algorithms in the literature use up to 50,000 function evaluations to solve benchmark problems and such a high number may not be realistic in solving real problems. Although the efficiency of different algorithms with the number of function evaluations is not visible in most of the articles, one can look at the maximum number of function evaluations used to choose an algorithm.

Issues with the present algorithms
Let us discuss next the main issues we have observed in the algorithms in the literature related to using an approximation in an EA and the numerical settings used to test their efficacy. These main issues include: 1. Type of test problems used (benchmark or real world), 2. Dimensions of the test problems in both objective and decision spaces considered, 3. Scarce use of constrained problems for numerical testing, 4. Neglecting the training time for fitting the metamodel used when reporting results, 5. Less focus on the accuracy of the metamodel, 6. Not well-explored updating criterion of the metamodel, 7. Not well-structured ensemble of metamodels, and 8. Less emphasis to algorithms that consider problem information when available.
To augment the discussion on the issues enumerated above, we present in Table 2 an overview of the test problems considered in the literature with their dimensions in both objective and decision spaces and the type of metamodel used. Issues 1-3 and 4-8 address the shortcomings in numerical testing and solution algorithms used, respectively. In Table 2, it can be clearly seen that both benchmark and real-world problems have been widely used by researchers to test their proposed EAs. Benchmark problems available in the literature are designed to be simple to implement, fast to evaluate, have known global optimal solutions, pose varied challenges to EAs such as several locally optimal solutions, etc. Benchmark problems solved are not computationally expensive and used just to test the efficiency of a particular algorithm. For more details about the characteristics of benchmark problems, see (Deb 2001; Coello et al. 2007). On the other hand, real-world problems either are actual problems considered in the industry or an emulation of them. Often, real-world problems are computationally expensive to evaluate, have uncertainties in their input and output variables, which in turn does not allow a thorough testing of an EA in practice using several real-world problems. Thus, it could be a good practice, if a numerical test setting involving several benchmark and a few real-world problems could be used to thoroughly test EAs. However, in the literature very few researchers have adopted this practice. One reason for this could be attributed to the lack of open availability of real-world problems to all researchers.
Real-world problems in industries often involve a large number of objectives, a large decision space and a large number of constraints. It is shown in Table 2 that most of the problems considered in the literature had low dimensions in both decision and objective spaces. Additionally, only six algorithms in Table 2 considered constraints besides the bounds for the decision variables. Thus, significant attention is needed toward issues 2 and 3 and more constrained benchmark and real-world problems need to be included in the numerical test settings. The number of objectives considered did not usually exceed three, and thus, the applicability of existing surrogate-based algorithms is limited. In the literature, the term many-objective optimization is also used as a synonym for optimizing (usually) more than three objectives. Moreover, the maximum number of decision variables considered is under 60. In fact, only three algorithms were used on more than three objective functions, four algorithms were used for three objective functions and the rest of the algorithms were used for solving biobjective optimization problems.
Furthermore, it is worth noting the connection between the type of metamodels used and the number of decision variables associated with the problem. In one instance, neural network was used when the number of decision variables was high (60). Next, in the decreasing order of popularity are algorithms using ensemble of metamodels, Kriging, support vector regression and polynomial approximation (including linear regression) for handling high-dimensional (decision space) problems. What is missing in the literature is a study with each of the algorithms about the efficacy of using different metamodels when the dimensions of the decision space increases, especially in algorithms where the type of the metamodel to be used within the algorithm is not binding. Hence, further research toward dimensions and types of metamodels is needed.
The issues 4-7 concern the use of metamodels. The major concern in the literature is how to alleviate the computational cost of the MOP. However, the training time needed to build/update a metamodel is completely neglected when reporting the results. This is mainly based on the assumption by the researchers that the computation time needed to compute objective and constraint function values is significantly higher as compared to the building/updating time of the metamodel used. This assumption may not hold in all cases and, specifically, when a large amount of data are used to build/update the metamodel.
Most of the algorithms discussed in this survey used a fixed evolution control strategy, i.e., the accuracy of the metamodel is not considered as an important issue. In this survey, we defined the adaptive evolution control, where if the metamodel is not accurate, the original functions are used for evaluating individuals. A desired accuracy for a metamodel can be predefined (e.g., in terms of a parameter) or changed by the user during the solution process. The accuracy can be calculated by different statistical measurements such as root mean square error, coefficient of determination, analysis of variance (ANOVA) (Tenne and Armfield 2008). As the accuracy of the metamodel is not considered in many algorithms in the literature, this is a valid future research direction.
The updating criterion of the metamodel can influence the quality of nondominated solutions. Considering the function approximation framework mentioned in Sect. 3, selected individuals from step 7 are re-evaluated with the original functions. The number and selection criterion of individuals to be used for re-evaluation from step 7 are not well explored in the literature; however, few algorithms considered both exploration and exploitation while updating the metamodel. We consider this as an important issue as selection of these individuals affects the updating time and accuracy of the metamodel, which may finally affect the quality of the set of nondominated solutions generated by the algorithm.
Although we consider an ensemble of metamodels to be a very promising element in the literature, we consider it as an issue as well. In most of the algorithms proposed in the literature, only the accuracy of the metamodel is used as the criterion to select a particular metamodel. However, there can be other criteria that can be considered in addition to the accuracy of the metamodel, e.g., diversity among individuals, training time of the metamodel. Considering these criteria for structuring an ensemble of metamodel is a future research direction.
Finally, the information about the MOP such as a different computation time required for evaluating different objective and/or constraint functions must be considered by the EA. Here we suggest two types of problem information that can be considered: firstly, if any objective function is not computationally expensive, a metamodel need not be created for that function. However, many real-world problems involve black box functions, and therefore, an explicit information about them can be difficult to obtain. Secondly, sometimes performing an experiment to calculate the value of an objective function may take less time or may give better results when compared to using a numerical model to evaluate the same. For example, Knowles (2009) performed a closedloop optimization, where values of objective functions are obtained by performing real experiments and selection while crossover and mutation operations are performed by using an EA. For more details about closed-loop evolutionary multiobjective optimization, see Knowles (2009). Thus, both types of problem information can lead to significant savings in computation time.

Conclusions
In this paper, we have presented a survey of different algorithms to reduce the computation time while solving computationally expensive multiobjective optimization problems. Altogether, 45 algorithms were found in the literature published in years 2008-2016, which were classified based on the type of approximation used, i.e., problem, function or fitness approximation. We found that function approximation or using surrogates is the most common approach for reducing the computation time. Different algorithms were summarized based on the steps of a unified function approximation framework involving an evolutionary algorithm. Additionally, six important challenges were identified, involving using the metamodel, updating the metamodel, training time, type of the metamodel to be used, when to update the metamodel and constraint handling for implementing the proposed approximation framework. Subsequently, the first three important challenges were used in the survey to highlight the differences among different algorithms.
A proper management of the metamodel makes an algorithm efficient to solve a computationally expensive problem. As different algorithms used different strategies to manage the metamodels and solved problems with different characteristics, it is difficult to generalize the efficiency of a method. However, the efficiency of an algorithm can be attributed by 1. Number of function evaluations used, 2. Dimensions (in both objective and decision spaces) of problems solved and 3. Characteristics of the problems solved, e.g., multimodality or disconnected Pareto front.
As metamodel-based algorithms are generally developed for black box problems, where characteristics of the problems to be solved are not known a priori, one can measure the efficiency of an algorithm by its ability to provide meaningful solutions in a least number of function evaluations. We also provided the dimensions and the number of function evaluations used to solve a particular problem.
We also compared these algorithms with respect to different criteria such as type of metamodel and evolutionary algorithm used, type of problem (benchmark or real-world) solved and evolution control. Some of the major findings were: 1. Kriging and neural networks were the most commonly used metamodels, 2. Most of the algorithms were based on dominance-based evolutionary algorithms, 3. Most of the algorithms solved the problems with no more than three objectives, 4. Most of the algorithms were not developed to handle constraints, 5. The number of decision variables was also limited especially when using Kriging, 6. Only few algorithms used an ensemble of metamodels, 7. Many algorithms were tested only on benchmark problems which were not at all computationally expensive.
We discussed problem and fitness approximation-based algorithms with respect to their potential toward reducing the computational complexity and the number of function evaluations. We also identified some promising elements and issues among algorithms in the literature. Promising elements involve using an ensemble of metamodels, hybrid algorithms with function approximation-based algorithms and different types of approximations together. We plan to address several issues observed related to shortcomings in numerical testings and algorithms in our future research.