An interactive surrogate-based method for computationally expensive multiobjective optimisation

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY-NC-ND 4.0 https://creativecommons.org/licenses/by-nc-nd/4.0/ An interactive surrogate-based method for computationally expensive multiobjective optimisation © 2018 Operational Research Society Published version Tabatabaei, Mohammad; Hartikainen, Markus; Sindhya, Karthik; Hakanen, Jussi; Miettinen, Kaisa Tabatabaei, M., Hartikainen, M., Sindhya, K., Hakanen, J., & Miettinen, K. (2019). An interactive surrogate-based method for computationally expensive multiobjective optimisation. Journal of the Operational Research Society, 70(6), 898-914. https://doi.org/10.1080/01605682.2018.1468860 2019 Full Terms & Conditions of access and use can be found at https://www.tandfonline.com/action/journalInformation?journalCode=tjor20 Journal of the Operational Research Society ISSN: 0160-5682 (Print) 1476-9360 (Online) Journal homepage: https://www.tandfonline.com/loi/tjor20 An interactive surrogate-based method for computationally expensive multiobjective optimisation Mohammad Tabatabaei, Markus Hartikainen, Karthik Sindhya, Jussi Hakanen & Kaisa Miettinen To cite this article: Mohammad Tabatabaei, Markus Hartikainen, Karthik Sindhya, Jussi Hakanen & Kaisa Miettinen (2019) An interactive surrogate-based method for computationally expensive multiobjective optimisation, Journal of the Operational Research Society, 70:6, 898-914, DOI: 10.1080/01605682.2018.1468860 To link to this article: https://doi.org/10.1080/01605682.2018.1468860 © 2018 Operational Research Society Published online: 20 May 2018. Submit your article to this journal Article views: 421 View Crossmark data https://doi.org/10.1080/01605682.2018.1468860 OPEN ACCESS An interactive surrogate-based method for computationally expensive multiobjective optimisation Mohammad Tabatabaei, Markus Hartikainen, Karthik Sindhya, Jussi Hakanen and Kaisa Miettinen University of Jyvaskyla, Faculty of Information Technology, Finland ABSTRACT Many disciplines involve computationally expensive multiobjective optimisation problems. Surrogate-based methods are commonly used in the literature to alleviate the computational cost. In this paper, we develop an interactive surrogate-based method called SURROGATE-ASF to solve computationally expensive multiobjective optimisation problems. This method employs preference information of a decision-maker. Numerical results demonstrate that SURROGATEASF efficiently provides preferred solutions for a decision-maker. It can handle different types of problems involving for example multimodal objective functions and nonconvex and/or disconnected Pareto frontiers. ARTICLE HISTORY Received 4 October 2016 Accepted 20 April 2018 KEYWORDS Multiple criteria decision-making (MCDM); interactive methods; computational cost; black-box functions; metamodeling techniques; achievement scalarising function 1. Introduction Many multiobjective optimisation problems (MOPs), arising e.g., in engineering applications, often involve multipleincommensurable,highlynonlinear,black-box and/or multimodal objective functions. Function evaluations in such problems may be conducted through time-consuming experiments and/or simulators. In the literature, these problems are called computationally expensivemultiobjectiveoptimisationproblems.MOPs typically have several (in many cases infinitely many) optimal solutions known as Pareto optimal solutions. The set of all Pareto optimal solutions in the objective space (called Pareto frontier) can be nonconvex and/or disconnected. From a mathematical point of view and without any preference consideration, Pareto optimal solutions are equally acceptable for an MOP. Therefore, when solving an MOP, a decision-maker (DM) is required to provide preference information and to choose his/her preferred solutions. Multiobjective optimisation methods are often categorised according to a decision-maker’s role in the solution process, i.e., non-interactive and interactive methods (Miettinen,1999). According to Miettinen (2008), in non-interactive methods, the DM either is not involved or provides preference information before or after the actual solution process. In interactive methods, the DM plays an essential role and the intention is to support him/her in the search for the most preferred solution. In such methods, steps of an iterative solution algorithm are repeated and the DM progressively specifies preference information so that the most preferred CONTACT Mohammad Tabatabaei [email protected] solution can be found. Examples of types of specifying preference information are reference points and classification of objective functions. What is noticeable is thattheDMcandetermineand alter his/her preferences between each iteration and in the meantime learn about the interdependencies in the problem as well as about one’s own preferences. This is a significant advantage of interactive methods in light of the fact that becoming more acquainted with the problem, its possibilities and limitations is often very valuable for the DM. See Miettinen (2008), Miettinen, Ruiz, and Wierzbicki (2008) for details of non-interactive and interactive methods. Surrogate-based methods are commonly used in the literature to alleviate computational cost (Jones,2001; Regis & Shoemaker,2012;Audet,2014;Tan,2015;Dengiz et al.,2009;Chih,2013;ReisdosSantos&Reis dos Santos,2011;Hurrion,2000;Kleijnen & van Beers, 2013;Van Beers & Kleijnen,2003;Yang & Tseng,2002; Mehdad & Kleijnen,2015;Park, Yum, Hung, Jeong, & Jeong,2016). The basic idea in such methods is to introduce a computationally less expensive problem called a surrogate problem and to replace the original problem with the surrogate one. In the literature, methods have been developed in which surrogate problems are built independently of the type of the optimisation algorithmsemployed in them. In Messac andMullur(2008), a non-interactive surrogate-based method is developed. This surrogate-based method does not consider the role of a DM and cannot handle multimodal functions or disconnected Pareto frontiers. The surrogate-based methods proposed in Yun, Yoon and Nakayama (2009), © 2018 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http://creativecommons.org/licenses/ by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way. JOURNAL OF THE OPERATIONAL RESEARCH SOCIETY 2019, VOL. 70, NO. 6, 898–914 ORIGINAL ARTICLE Revised 23 January 2018 Kitayama, Srirat, Arakawa, and Yamazaki (2013) consider the role of a DM in the solution process. However, they are not interactive methods. This means that if a DMwishes to provide newpreferences, the entire methods should be run again. Therefore, the DM should wait for a long time to get new solutions corresponding to his/her preferences. In (Hartikainen, Miettinen, &Wiecek,2012;Hartikainen & Lovison,2015;Ruiz, Sindhya,Miettinen, Ruiz, & Luque,2015)methods have been developed by which the DM can find preferred solutions quickly. These methods require a set of approximated solutions a priori generated by some other surrogate-based methods. InTabatabaei,Hakanen,Hartikainen,Miettinen,and Sindhya (2015), we surveyed surrogate-based methods, and observed some shortcomings including inability to (1) involve a DM in the solution process, (2) deal with multimodal functions and (3) capture a nonconvex and disconnected Pareto frontier. To overcome these shortcomings, in this paper, we develop an interactive surrogate-based method. In our method, a DM provides his/her preferences in the form of a reference point containing aspiration levels representing desirable values for objective functions. According to Larichev (1992), the type of preference information in the form of a reference point has been regarded to be understandable for a DM. The method proposed is called SURROGATE-ASF and involves two phases, i.e., initialisation and decisionmaking phases. In the initialisation phase, a set of nondominatedsolutionsinthedecisionandobjectivespaces are generated. Using these solutions, hyper-boxes are formed in the decision space. Corresponding to each hyper-box, a single objective surrogate function is built by approximating an achievement scalarising function (ASF)(Wierzbicki,1986). In the literature, this function is used to compute the (weakly) Pareto optimal solution for a given reference point. In the decision-making phase, the single objective surrogate functions built in the initialisation phase are employed in generating solutions reflecting the preferences of the DM. One should note that the DM is only involved in the latter phase.BasedonthecomparisoninMuller and Shoemaker (2014), a cubic radial basis function (RBF) with a linear tail is employed as the metamodeling technique in this paper, but SURROGATE-ASF is not limited to RBF. SURROGATE-ASF has been developed for computationally expensive MOPs with box constraints. The novelty of SURROGATE-ASF can be summarised as follows: (1) the DM does not need to wait for a long time to obtain his/her solutions corresponding to his/her preferences, (2) SURROGATE-ASF utilises reference points given by a DM as input, (3) the DM can explore different regions, a particular interesting region or the entire Pareto frontier, (4) it provides approximated solutions in both the decision and objective spaces simultaneously. The rest of this paper is organised as follows. In Section 2, the basic concepts used in this paper are addressed. The SURROGATE-ASF method is presented in Section 3. Numerical results of evaluating the performance of SURROGATE-ASF on a practical shape optimisation problem of designing an airfoil as well as some benchmark problems are presented in Section 4. Finally, in Section 5, we draw our conclusions and discuss paths for future research. 2. Basic concepts In this section, we introduce the concepts and notations used in this paper. We consider multiobjective optimisation problems of the form: minimise x∈S{f1(x),...,fk(x)},(1) where fi:S→Rare k(≥2)conflicting, computationally expensive objective functions, S={x∈Rn: xl c≤xc≤xu c,c=1, ...,n}is a nonempty feasible decision set which is a subset of the decision space Rn. Asolutionx=(x1,...,xn)T∈Sis called a feasible decision (variable) vector,wherexc,c=1, ...,n, are decision variables and, xl cand xu care the lower and upper bounds of xc, respectively. The image of xin the objective space Rkis called a feasible objective vector denoted by f(x).TheimageofSin the objective space is called the feasible objective set denoted by Z(=f(S)). A feasible solution x∗∈Sand the corresponding f(x∗)∈Zare termed weakly Pareto optimal for problem (1), if there does not exist another feasible solution x∈Ssuch that fi(x)<f i(x∗)for all i=1, ...,k. Correspondingly, they are Pareto optimal for problem (1), if there does not exist another feasible solution x∈Ssuch that fi(x)≤fi(x∗)for all i=1, ...,k,and fj(x)<f j(x∗)for at least one index j∈{1, ...,k}.The set of all Pareto optimal solutions in the objective space is called a Pareto frontier. Let the set Xd={x1,...,xd} be an arbitrary subset of feasible solutions in S,and Fd={f(x1),...,f(xd)}, the corresponding objective vectors in Z. A solution xi(or f(xi)), i=1, ...,d, that satisfies the definition of Pareto optimality with respect to all solutions in Xd(or Fd), is called a nondominated solution in Xd(or Fd). A solution x(or f(x))iscalledalocally non-dominated solution if there exist a non-empty set S⊆Ssuch that, x(or f(x)) satisfies the definition of Pareto optimality with respect to all points in S(or f(S)). A Pareto optimal solution is a non-dominated solution, but a non-dominated one is not necessarily Pareto optimal. We also define a feasible solution xi e∈argminx∈S{fi(x)}for i=1, ...,k.Theith JOURNAL OF THE OPERATIONAL RESEARCH SOCIETY 899 extreme point (solution) for i=1, ...,k, is defined as zi e=f(xi e). These extreme points in the objective space are also extremes of the Pareto frontier. Based on the extreme points, a vector of lower bounds of the objective function values in the Pareto frontier is defined as the ideal (objective) vector and denoted by zideal =(zideal 1,...,zideal k)Twhere zideal i=fi(xi e)for i=1, ...,k.Theutopian (objective) vector zutp is a vector in which its components are calculated by subtracting some small positive scalar (e.g., 10−6)from the components of zideal. A vector of upper bounds of the objective function values in the Pareto frontier is defined as the nadir (objective) vector and denoted by znadir =(znadir 1,...,znadir k)T. The components of the nadir vector can be approximated by e.g., a pay-off table usingtheextremepoints. More information of the ideal, utopian and nadir vectors is given e.g., in Miettinen (1999). We define the difference operator a−b= (a1−b1,...,ak−bk),wherea=(a1,...,ak),b= (b1,...,bk)∈Rk. In SURROGATE-ASF, the DM provides his/her preferences in the form of a reference point z∗=(¯z∗ 1,...,¯z∗ k)T,where¯z∗ iis an aspiration level representing a desirable value for the objective function fi.Onecanfindapreferredsolutionforareferencepoint given by a DM by applying an appropriate scalarisation. Scalarising problem (1) means formulating a single objective optimisation problem such that its (globally) optimalsolution is aParetooptimalsolutionfor (1).Inthis paper, we consider the following widely used achievement scalarising function (ASF) (Wierzbicki,1986)as an important element of SURROGATE-ASF: ASF: S×Rk→R (x,z∗)→ max i=1,...,k(wi(fi(x)−¯z∗ i)),(2) where wi≥0, for i=1, ...,k, are non-negative fixed weights which actually set a direction where z∗is projected onto the Pareto frontier. In this paper we set wi=1 znadir i−zutp i ,fori=1, ...,k, which are widely used (Wierzbicki,1986). A new (but still computationally expensive) single objective optimisation problem is formulated as minimise x∈SASF(x,z∗). (3) The reference point in problem (3) can be feasible or infeasible, i.e., inside or outside of the feasible objective set. As proved in Wierzbicki (1986), theoretically by solving problem (3), a (weakly) Pareto optimal solution corresponding to the reference point is obtained, regardless of the feasibility or infeasibility of the reference point. This does not hold for all other scalarising approaches (Miettinen,1999 ). Different (weakly) Pareto optimal solutions can be obtained by changing Figure 1. Projection of a given reference point onto the convex hull. the reference point in problem (3). One can add an augmentationtermto(2) to avoid weakly Pareto optimal solutions as discussed and proved in Miettinen (1999), Wierzbicki (1986). Numerically, in cases with a limited number of function evaluations, the optimal solution of problem (3) is a locally non-dominated solution. In what follows, an optimal solution of problem (3) corresponding to a reference point given by a DM and a predetermined reference point is called a preferred solution and a reference solution, respectively. Moreover, problem (1)is referred to as the original problem and its computationally expensive functions as original functions. 3. The SURROGATE-ASF method In the SURROGATE-ASF method, a DM can provide preferred ranges of the objective functions. Then, they are considered as the utopian and the nadir vectors corresponding to the region which is interesting for the DM in the objective space. Alternatively, estimates of the utopian and nadir vectors can be incorporated into the method. In fact, if the utopian and nadir vectors are given, then the DM can explore the entire Pareto frontier by providing different reference points through multiple interactions in the decision-making phasetobediscussed in Section 3.3. SURROGATE-ASF involves two phases, i.e., initialisation and decisionmaking phases. In the initialisation phase, a set of nondominated solutions within the preferred ranges given by the DM is generated. Then, by calling Algorithm 1 to be discussed in the following subsection, a finite number of hyper-boxes is formed using the non-dominated solutions in the decision space. For each individual hyper-box, a single objective surrogate function is built by approximating the ASF (2) where the decision variables of problem (1) and the aspiration levels appearing in the ASF are treated as variables of the surrogate function. In the decision-making phase to be discussed in Section 3.3, these computationally inexpensive surrogate functions are utilised to interact with the DM. 900 M. TABATABAEI ET AL. (a) (b) Figure 2. (a): A convex hull with 5 points for a bi-objective problem with q=4 divisions. (b): A convex hull with 15 points for a three-objective problem with q=4 divisions. 3.1. Surrogate building module If the DM provides the preferred ranges, the extreme points of the corresponding region in the objective space must be calculated. To be more specific, the corresponding solutions in the decision space are found by solving problem (3) using the extreme points as reference points. One should note that if the extreme points are formed based on the estimates of the utopian and nadir vectors, then their corresponding solutions in the decision space are already available. The approximation of the ASF requires sample points for the aspiration levels and the decision variables. However, before discussing the procedure of selecting sample points, we first look at an interesting property regarding problem (3). This property assists to identify regions in the decision and objective spaces where to select sample points. Suppose that His the convex hull of all individual extreme points. This is constructed by a convex combination of all extreme points in the objective space. Moreover, let z∗be a reference point given by a DM, and x∗be the corresponding optimal solution of problem (3). Thus, x∗∈argmin x∈S [max i=1,...,k(wi(fi(x)−z∗ i))]. On the other hand, x∗∈argmin x∈S [max i=1,...,k(wi(fi(x)− z∗ i)) +t],forallt∈R. Therefore, we have: x∗∈argmin x∈Smax i=1,...,k(wi(fi(x)−z∗ i)) +t =argmin x∈Smax i=1,...,k(wi(fi(x)−z∗ i)+t) =argmin x∈Smax i=1,...,k(wi(fi(x)−z∗ i+t wi )) =argmin x∈Smax i=1,...,k(wi(fi(x)−(z∗ i−t wi ))) =argmin x∈Smax i=1,...,k(wi(fi(x)−z∗∗ i)), where z∗∗ is treated as an arbitrary reference point on the ray passing through z∗in parallel with w= (w1,...,wk)T. It means that the optimal solution of problem (3) for a reference point z∗given by a DM is also the optimal solution of problem (3)forany reference point on the ray passing though z∗in parallel with w(see Figure 1). The following problem gives the projected reference point z∗∗ in H: z∗∗ ∈argmin h∈H    h−z∗ w   (4) where .is the Euclidean norm. Problem (4) applies if and only if there exists a real number t∈Rsuch that z∗=h+twfor some h∈H. Since any reference point within the preferred ranges of objective functions given by the DM can be projected onto the convex hull H(i.e., finding the closest point on Halong the direction w), sample points for the aspirationlevelsasreferencesamplepointsareselectedonthis convex hull. To do this, first, a set of evenly distributed predetermined reference points including the extreme points is generated on the convex hull Hand denoted by ZP. The reference solutions (non-dominated solutions) in the decision space corresponding to these reference points are obtained by solving problem (3)foreach individual predetermined reference point. This is conducted by employing a surrogate-based single objective optimisation method. The choice of this method has an impact on the performance of SURROGATE-ASF. One can use an appropriate state-of-the-art method to solve problem (3) efficiently. These reference solutions are evaluated with the original functions. In order to generate the predetermined reference points (without involving any DM), we use the method presented in Das and Dennis (1998) that places points on a convex hull of all individual extreme points (a (k−1)- dimensional simplex). If hzpredetermined reference points (including the extreme points in the objective space) are considered, the qdivisions along JOURNAL OF THE OPERATIONAL RESEARCH SOCIETY 901 each objective coordinate axis in the objective space of ak-objective problem can be calculated using: hz=k+q−1 q.(5) Figures 2(a) and 2(b) depict two convex hulls for biand three-objective problems with q=4andhz= 2+4−1 4=5andhz=3+4−1 4=15 predetermined reference points (black circles), respectively. The sets of reference solutions in the decision and the objective spaces are denoted by XPand FP, respectively. By making use of the set of predetermined reference points ZPand reference solutions XP, the convex hull Hand the decision space are decomposed into a finite number of sub-regions and hyper-boxes, respectively. In what follows, a=1, ...,r, represents the index corresponding to the ath hyper-box or sub-region where r is the number of hyper-boxes (sub-regions). For k=2 and 3, we have r=qand q2, respectively. Algorithm 1 starts by constructing sub-regions on H using the predetermined reference points in ZP. Each sub-region is formed by selecting the knearest neighbor points on Hstarting from one of the extreme points. Figures 2(a) and 2(b) represent one possible way of forming the sub-regions and numbering them where k=2andk=3, respectively. Then, reference sample points (the stars in Figures 2(a) and 2(b)) are selected within these sub-regions. To select reference sample points,the predetermined reference pointscorresponding to the ath sub-region (i.e., ZP,a={zP,a,1,...,zP,a,k}) are considered as the extreme points of this sub-region. Then,by using the method presented in Das and Dennis (1998), a set of evenly distributed reference sample points within this sub-region on the convex hull is generated. The larger the number of generated reference sample points, the higher is the accuracy of the surrogate function  ASFacorresponding to the ath subregion. These reference sample points along with the reference points in ZP,aare considered as sample points corresponding to this sub-region and denoted by Za. To decompose the decision space, hyper-boxes are built corresponding to sub-regions on the convex hull H. Considering the ath sub-region on Hand the set of predetermined reference points forming this subregion, i.e., ZP,a={zP,a,1,...,zP,a,k}, the set of nondominated solutions corresponding to these predetermined reference points in the decision and the objective spaces are denoted by XP,a={xP,a,1,...,xP,a,k}and FP,a={f(xP,a,1),...,f(xP,a,k)}, respectively. Then, the lower and the upper bounds of the decision variable xc, for c=1, ...,n, within the ath hyper-box denoted by Sacorresponding to this sub-region are calculated as follows: Figure 3. Hyper-boxes in the decision space corresponding to the sub-regions on the convex hull in Figure 2(a). For example, reference solutions xP,1,1 and xP,1,2 form S1. min{xP,a,1 c,...,xP,a,k c}≤xc≤max{xP,a,1 c,...,xP,a,k c}, for all c=1, ...,n.(6) Figure 3shows a simple example of hyper-boxes for a problem with two decision variables. Once Sais formed, a small number of initial sample points within this hyper-box is selected using some sampling technique such as Latin hypercube sampling (LHS) (Helton, Johnson, Sallaberry, & Storlie,2006) and evaluated with the original functions. These evaluated points along with the non-dominated solutions in XP,aare considered as the initial sample points in the decision space corresponding to Saand denoted by Xa. To build the surrogate function  ASFacorresponding to the hyper-box Sa, a cubic RBF with a linear tail is employed to approximate ASF (its description is given 902 M. TABATABAEI ET AL. e.g., in Muller, Shoemaker, and Piche (2013)). To do this, the following Cartesian product of sets Xaand Za is formed as input data denoted by Tafor RBF: Ta=Xa×Za={(x,z)|x∈Xaand z∈Za}.(7) The ASF values of all elements in Taare calculated. One should note that the objective function values of all sample points in Xaare available. Then, Tais divided into validation and training sets using the leaveone-out cross-validation (LOOCV) or the k-fold crossvalidation method discussed in Muller and Shoemaker (2014) (note that it differs from the index kfor the kth objective function). Let |Ta|be the cardinality of Ta. According to Muller and Shoemaker (2014), if |Ta|≤ 50, we use LOOCV. If |Ta|>50, then we use the kfold cross-validation method where k=10, for 50 < |Ta|≤100, k=20, for 100 <|Ta|≤150 and k=30, for 150 <|Ta|≤200. Having the elements in the training set and their corresponding ASF values as input data, RBF is employed to build  ASFa. When using either LOOVC or the k-fold cross-validation method, the prediction accuracy of  ASFais evaluated using standard error measures like Root Mean Squared Error (RMSE) (Giunta, Watson, & Koehler,1998) and/or R2 (Jin, Chen, & Simpson,2001). The steps of the decomposition procedure of the decision space into a finite number of hyper-boxes and building corresponding surrogate functions are shown in Algorithm 1. Note that the prediction accuracy of  ASFamay not be satisfactory after steps 1-10 of Algorithm 1 and, in that case, it needs to be improved. In the following subsection, we discuss an update strategy called Algorithm 2 by which new sample points within the hyper-box Sacan be selected. The process of updating  ASFais repeated until a desired accuracy is achieved. The non-dominated solutions used to form hyperboxescorrespondtothepredeterminedreferencepoints on the convex hull. It may happen that at least two predetermined reference points have the same corresponding non-dominated solution. In such cases, to form sub-regions and hyper-boxes, out of all predetermined reference points with the same non-dominated solution (i.e., the same objective values but different decision variable values), only one reference point (selected arbitrarily) and its corresponding non-dominated solution are considered along with other different predetermined reference points and non-dominated solutions. Then, sub-regions and hyper-boxes are formed. It may also occur that hyper-boxes overlap each other. In these cases, some sample points selected fall into more than one hyper-box. However, sub-regions and hyper-boxes are still formed as mentioned in Algorithm 1. In practice, obtaining a desired accuracy may not be possible. In such a case, depending on the computational budget, one can consider a maximum number of function evaluations for each individual hyper-box. Then, either achieving a desired accuracy or reaching the maximum number of function evaluations can be considered as a stopping criterion. 3.2. Accuracy improvement module To improve the accuracy prediction of  ASFacorresponding to Sa, we discuss Algorithm 2. When updating  ASFa, two matters are considered, i.e, sampling non-dominated solutions within Saand having a good diversity among the sample points in Xa. In order to sample non-dominated solutions, we employ the sampling function discussed in Kitayama et al. (2013) to generate new sample points within Sa. The authors of Kitayama et al. (2013) modify the Pareto fitness function (Schaumann, Balling, & Day,1998)and introduce a sampling function by using RBF such that by minimising this function, the optimal solution can be a non-dominated solution. See Kitayama et al. (2013) for details of building this sampling function. For the sake of simplicity, we refer to this sampling function as the modified Pareto fitness function (MPF). By minimising the MPF iteratively, non-dominated solutions within Sacan be generated. Nevertheless, the JOURNAL OF THE OPERATIONAL RESEARCH SOCIETY 903 Figure 4. Selecting sample points within Sato update  ASFa. diversity of the points in Xato cover the hyper-box should be considered. Thus, in some iteration of updating  ASFa, some solution obtained by minimising the MPF is not evaluated with the original functions. This requires a criterion to assess whether the optimal solution of the MPF should be evaluated with the original functions and added into Xaor not. We discuss this criterion here. Since some points are not evaluated with the original functions, we introduce XPaand FPaas sets of archived points generated during the updating process in the decision and the objective spaces, respectively. These sets include all points in Xaand Fa(the points evaluated with the original functions) and sample points which are not evaluated with the original functions. In what follows, a procedure is discussed to assign objective function values to the points not evaluated with the original functions. Building the MPF is based on the points in XPaand FPa. To do this, all points in XPaand their corresponding MPF values are considered as input data to train a cubic RBF with a linear tail. Once the MPF is built, it is minimised and the optimal solution xcand within Sa is obtained. Then, as depicted in Figure 4, the closest point in Xain terms of the Euclidean distance to xcand is found. This point is denoted by xcand close , the corresponding objective vector by fcand close and the Euclidean distance between xcand and xcand close by dcand close . The closest point in Xa\{xcand close }to xcand close is also found. We denote this point by xclose close and the Euclidean distance between xclose close and xcand close by dclose close. To have diversity among sample points within Xa, the criterion dcand close >dclose close 2is checked. This means that if xcand is outside the circle centered on xcand close with radius dclose close 2(see Figure 4), then xcand is evaluated with the original functions and the objective vector fcand is calculated. Before adding xcand and fcand into Xaand Fa, respectively, the prediction accuracy of  ASFais assessed. To accomplish this, a test set as defined in Ta(7) is formed where xcand and all reference sample points in Zaare utilised to form the Cartesian product set. All elements in the test set are evaluated with  ASFaand ASF, and the prediction accuracy of  ASFais checked. If the accuracy is not acceptable, xcand is added into Xa and XPaand fcand into Faand FPa, respectively. Then, the points in Xaare utilised to update  ASFaand the points in XPato update the MPF. The process of minimising the MPF and checking the diversity criterion is repeated. If the prediction accuracy is acceptable, then the process of updating  ASFais terminated. If dcand close ≤dclose close 2,xcand is inside the circle centered on xcand close with radius dclose close 2,thenxcand is not evaluated with the original functions and is not added into Xa.However, to generate new sample points by the MPF, the objective vector corresponding to xcand is required. In this case, to avoid generating xcand again, the objective vector corresponding to this point is artificially treated as a dominated point in the objective space by the objective vector fcand close . To do this, an artificial objective vector denoted by fcand artif is considered for this point by adding a random vector with positive components to the objective vector fcand close . This random vector is a fraction of znadir −fcand close , e.g., znadir−fcand close 10 .Bythis choice, fcand close is non-dominated with respect to fcand artif , and is inside the feasible objective space. Then, xcand and fcand artif are added into XPaand FPa. Note that these points are not added into Xaand Fa.Then,theMPF values of all points in XPaare calculated and the MPF is updated. The updating process is repeated until the stopping criterion regarding the accuracy of  ASFais met. 3.3. SURROGATE-ASF and decision-making In this subsection, we present both the decision-making phase as well as the SURROGATE-ASF algorithm. As discussed in Section 3.1, in the initialisation phase, surrogate functions are built by calling Algorithm 1. These surrogate functions are employed to interact with the DM. In the decision-making phase, if the DM did not provide any preferred ranges, the utopian and the nadir vectors are shown to him/her. The DM is supposed to specify reference points within the ranges. In any case, the DM provides a reference point ¯z∗. The interaction with the DM in SURROGATE-ASF corresponds to that of the reference point method given in Wierzbicki (1982). The reference point ¯z∗is projected onto the convex hull Hby solving problem (4) as discussed in Section 3.1. In practice, to find the projected reference point denoted by z∗∗, we do not solve problem (4). A simple approach is to generate a large number of uniformly distributed points on H. Then, the point 904 M. TABATABAEI ET AL. MATSuMoTo performed better than SURROGATEASF with DIRECT. Overall, considering μcand σcgiveninTables2a, 3a,3b,4a and 4b, SURROGATE-ASF (based on both MATSuMoTo and DIRECT) performed very well in the numerical tests when compared to the two other methods. However, only in Table 2b, the Yun method outperformed SURROGATE-ASF based on DIRECT. In addition, SURROGATE-ASF based on MATSuMoTo performed better than SURROGATE-ASF based on DIRECT on the benchmark problems considered. In terms of computational burden we only report the maximum CPU time among all CPU time taken by the methods to solve the benchmark problems. SURROGATE-ASF with MATSuMoTo and DIRECT required at most 22 and 18 seconds, respectively to build surrogate functions in the initialisation phase. In the decision-making phase, the DM found the preferred solution corresponding to each reference point in less than 1 second. In PAINT with MATSuMoTo and DIRECT, the surrogate problem was built in 183 and 169 seconds, respectively. Then, the DM was able to find the preferred solutions in less than 1 seconds for all reference point except the last once. In the last iteration, to see the corresponding preferred solution, the DM waited for 4 seconds. In the Yun method, the DM managed to see the preferred solution for each reference point in 255 seconds. As can be understood, SURROGATE-ASG with MATSuMoTo and DIRECT outperformed other methods in terms of CPU time to provide preferred solutions for the DM. 5. Conclusions and future research directions In this paper, we developed an interactive surrogatebased method called SURROGATE-ASF to solve computationally expensive MOPs. This method consists of two phases: initialisation and decision-making phases. In the initialisation phase, the decision space is decomposed into a finite number of hyper-boxes. For each hyper-box, a single objective surrogate of the achievement scalarising function is built by using e.g., a cubic RBF with a linear tail. In the decision-making phase, iterations with a decision-maker are conducted in an interactive fashion. In this phase, at each iteration a reference point is specified by the decision-maker. Then, a surrogate single objective optimisation problem is formulated and solved by any appropriate single objective optimisation method. The optimal solution is an approximation of the preferred solution in the decision space. This solution is evaluated with the original functions and shown to the decision-maker. The interaction with the decision-maker is repeated until the most preferred solution for him/her is found. In SURROGATE-ASF, the approach of building a computationally inexpensive surrogate of the achievement scalarising function is applied instead of a typical approach of building a surrogate of each individual computationally expensive objective function and then forming the achievement scalarising function as discussed in Yun et al. (2009). Therefore, when interacting with the DM, only a single objective optimisation problem is solved. The surrogate assists the DM to find his/her preferred solutions in both the decision and the objective spaces quickly. Numerical results confirmed that SURROGATE-ASF performed very well, in particular, on problems with multimodal objective functions, non-convex and/or disconnected Pareto frontier and it fills a gap in the selection of methods available for computationally expensive problems. Solving a real-world computationally expensive airfoil optimisation problem demonstrated that SURROGATE-ASF reduced the computational burden significantly in comparison with the typical approach by which the DM had to wait for a long time to find the most preferred solution. The speed can be further improved by applying parallel computing when building single objective surrogate functions. By increasing the number of objective functions, the number of sub-problems is also increased. Therefore, with a fixed budget of function evaluations, the number of function evaluations for each sub-problem is decreased. This may lead to a sacrificed accuracy of the sub-problems. A possible idea is to merge hyper-boxes and sub-regions to reduce the number of sub-problems. This is a future research direction when SURROGATEASF is to be applied for problems with a large number of objective functions. SURROGATE-ASF with the current setting along with an appropriate single objective surrogate-based method can be employed for problem with a high-dimensional decision space. However, the accuracy of the surrogate functions for each hyper-box may not be satisfactory. A possible approach is to apply a dimensionality reduction method within each hyper-box to obtain a low-dimensional hyper-box. Then, SURROGATE-ASF with the current format can be applied for each new hyper-box. This is another future research direction. Abbreviations and Notations MOP Multiobjective Optimization Problem DM Decision-Maker ASF Achievement Scalarising Function RBF Radial Basis Function SVR Support Vector Regression LOOCV The leave-one-out cross-validation MPF Modified Pareto Fitness kNumber of (computationally expensive) objective functions nNumber of decision variables fiith computationallyexpensiveobjectivefunction SFeasible decision set JOURNAL OF THE OPERATIONAL RESEARCH SOCIETY 911 xDecision (variable) vector xccth decision variable xl c,xu cLower and upper bounds of xc RkObjective space f(x)Objective vector Z=f(S)Feasible objective set zi eith extreme point zideal Ideal (objective) vector zutp Utopian (objective) vector znadir Nadir (objective) vector z∗Referencepointgivenbythedecision-maker HConvex hull of all individual extreme points z∗∗ Projected reference point on Hcorresponding to z∗ ZPSet of predetermined reference points generated on H XPSet of reference solutions in the decision space FPSet of reference solutions in the objective space qNumber of divisions along each objective coordinate axis rNumber of sub-regions (or hyper-boxes) Saath hyper-box zP,a,iith predetermined reference point forming the ath sub-region on H ZaSet of reference sample points corresponding to the ath sub-region XP,aSet of non-dominated solutions in the decision space corresponding to the predetermined reference points in ZP FP,aSet of non-dominated solutions in the objective space corresponding to the predetermined reference points in ZP xP,a,iith of the non-dominated solutions defining Sa XaSet of sample points in the decision space corresponding to Sa  ASFaSurrogate function corresponding to Sa TaCartesian product of sets Xaand Za XPaSet of archived points in the decision space generated during the updating process FPaSet of archived points in the objective space generated during the updating process xcand Minimiser of MPF xcand close Closest point in Xato xcand xclose close Closest point in Xa\{xcand close }to xcand close dcand close Euclidean distance between xcand and xcand close dclose close Euclidean distance between xcand close and xclose close fcand close Objective vector corresponding to xcand close fcand artif Artificial objective vector Acknowledgements This work was partly funded by the COMAS Doctoral Program at the University of Jyvaskyla, the Academy of Finland [project No. 287496], Early Career Scheme (ECS) sponsored by the Research Grants Council of Hong Kong [project No. 21201414 (Dr. Matthias Hwai Yong Tan)] and the KAUTE Foundation.MohammadTabatabaei thanksDr. AlfredoArias Montano for providing the airfoil optimisation problem code and Prof. Raino Mäkinen for fruitful discussions on this problem. References Arias-Montano, A., Coello Coello, C.A. & Mezura-Montes, E. (2012). Multi-objective airfoil shape optimization using a multiple-surrogate approach. Proceedings of the 2012 IEEE Congress on Evolutionary Computation (pp. 1–8). IEEE. Audet, C. (2014). A survey on direct search methods for blackbox optimization and their applications. In P. M. Pardalos & T. M. Rassias (Eds.), Mathematics without boundaries (pp. 31–56). New York, NY: Springer. Chih, M. (2013). A more accurate second-order polynomial metamodel using a pseudorandom number assignment strategy. Journal of the Operational Research Society, 64(2), 198–207. Das, I., & Dennis, J. (1998). Normal-boundary intersection: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems. SIAM Journal on Optimization, 8(3), 631–657. Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002). Scalable multi-objective optimization test problems. In Proceedings of the 2002 IEEE Congress on Evolutionary Computation (Vol. 1, pp. 825–830). IEEE. Dengiz, B., Alabas-Uslu, C., & Dengiz, O. (2009). Optimization of manufacturing systems using a neural network metamodel with a new training approach. Journal of the Operational Research Society, 60(9), 1191–1197. Drela, M. (1989). Xfoil: An analysis and design system for low Reynolds number aerodynamics. Giunta, A., Watson, L. T., & Koehler, J. (1998). A comparison of approximation modeling techniques: Polynomial versus interpolating models. In Proceedings of 7th AIAA/USAF/NASA/ISSMO symposium on Multidisciplinary Analysis and Optimization (Vol. 1, pp. 392–404). St. Louis, MO. Goldberg, D. E. (1989). Genetic algorithms in search, optimization and machine learning. Boston: AddisonWesley Longman Publishing Co., Inc. Hartikainen, M., & Lovison, A. (2015). Paint-SiCon: constructing consistent parametric representations of Pareto sets in nonconvex multiobjective optimization. Journal of Global Optimization, 62(2), 243–261. Hartikainen, M., Miettinen, K., & Wiecek, M. M. (2012). PAINT: Pareto front interpolation for nonlinear multiobjective optimization. Computational Optimization and Applications, 52(3), 845–867. Helton,J.C.,Johnson,J.D.,Sallaberry,C.J.,&Storlie,C.B. (2006). Survey of sampling based methods for uncertainty and sensitivity analysis. Reliability Engineering & System Safety, 91(10–11), 1175–1209. Hurrion,D.R.(2000). A sequential method for the development of visual interactive meta-simulation models using neural networks. Journal of the Operational Research Society, 51(6), 712–719. Jin, R., Chen, W., & Simpson, T. W. (2001). Comparative studies of metamodelling techniques under multiple modelling criteria. Structural and Multidisciplinary Optimization, 23(1), 1–13. 912 M. TABATABAEI ET AL. Jones, D. R. (2001). A taxonomy of global optimization methods based on response surfaces. Journal of Global Optimization, 21(4), 345–383. Jones, D. R., Perttunen, C. D., & Stuckman, B. E. (1993). Lipschitzian optimization without the Lipschitz constant. Journal of Optimization Theory and Applications, 79(1), 157–181. Kitayama, S., Srirat, J., Arakawa, M., & Yamazaki, K. (2013). Sequential approximate multi-objective optimization using radial basis function network. Structural and Multidisciplinary Optimization, 48(3), 501–515. Kleijnen, C. J. P., & van Beers, M. W. C. (2013).Monotonicitypreserving bootstrapped kriging metamodels for expensive simulations. Journal of the Operational Research Society, 64(5), 708–717. Larichev, O. I. (1992). Cognitive validity in design of decisionaiding techniques. Journal of Multi-Criteria Decision Analysis, 1(3), 127–138. Mehdad, E., & Kleijnen, C. J. P. (2015). Classic Kriging versus Kriging with bootstrapping or conditional simulation: classic Kriging’s robust condence intervals and optimization. Journal of the Operational Research Society, 66(11), 1804–1814. Messac, A., & Mullur, A. A. (2008). A computationally efficient metamodeling approach for expensive multiobjective optimization. Optimization and Engineering, 9(1), 37–67. Miettinen, K. (1999). Nonlinear multiobjective optimization. Boston: Kluwer Academic Publishers. Miettinen, K. (2008). Introduction to multiobjective optimization: Noninteractive approaches. In J. Branke , K. Deb, K. Miettinen, & R. Slowinski (Eds.), Multiobjective Optimization: Interactive and evolutionary approaches (pp. 1–26). Berlin Heidelberg: Springer-Verlag. Miettinen, K., Ruiz, F., & Wierzbicki, A. P. (2008). Introduction to multiobjective optimization: Interactive approaches. In J. Branke, K. Deb, K. Miettinen, & R. Slowinski (Eds.), Multiobjective optimization: Interactive and evolutionary approaches (pp. 27–57). Berlin Heidelberg: Springer-Verlag. Muller, J., & Shoemaker, C. A. (2014). Inuence of ensemble surrogate models and sampling strategy on the solution quality of algorithms for computationally expensive blackbox global optimization problems. Journal of Global Optimization, 60(2), 123–144. Muller, J., Shoemaker, C. A., & Piche, R. (2013). SO-MI: A surrogate model algorithm for computationally expensive nonlinear mixed-integer black-box global optimization problems. Computers & Operations Research, 40(5), 1383– 1400. Park, T., Yum, B., Hung, Y., Jeong, Y.-S., & Jeong, K. M. (2016). Robust kriging models in computer experiments. Journal of the Operational Research Society, 67(4), 644– 653. Regis, R. G., & Shoemaker, C. A. (2012). Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization. Engineering Optimization, 45(5), 529–555. Reis dos Santos, I. M., & Reis dos Santos, M. P. (2011). Construction and validation of distributionbased regression simulation metamodels. Journal of the Operational Research Society, 62(7), 1376–1384. Ruiz, A. B., Sindhya, K., Miettinen, K., Ruiz, F., & Luque, M. (2015). E-NAUTILUS: A decision support system for complex multiobjective optimization problems based on the nautilus method. European Journal of Operational Research, 246(1), 218–231. Schaumann, E., Balling, R., & Day, K. (1998). Genetic algorithms with multiple objectives. Proceedings of 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization (Vol. 3, pp. 2114–2123). St. Louis, MO. Tabatabaei, M., Hakanen, J., Hartikainen, M., Miettinen, K., & Sindhya, K. (2015). A survey on handling computationally expensive multiobjective optimization problems using surrogates: non-nature inspired methods. Structural and Multidisciplinary Optimization, 52(1), 1–25. Tan, M. (2015). Sequential Bayesian polynomial chaos model selection for estimation of sensitivity indices. SIAM/ASA Journal on Uncertainty Quantification, 3, 146–168. Van Beers, M. W. C., & Kleijnen, C. J. P. (2003). Kriging for interpolation in random simulation. Journal of the Operational Research Society, 54(3), 255–262. Wierzbicki, A. P. (1982). A mathematical basis for satisficing decision making. Mathematical Modelling, 3(5), 391–405. Wierzbicki, A. P. (1986). On the completeness and constructiveness of parametric characterizations to vector optimization problems. OR Spectrum, 8(2), 73–87. Yang, T., & Tseng, L. (2002). Solving a multi-objective simulation model using a hybrid response surface method and lexicographical goal programming approach-a case study on integrated circuit ink-marking machines. Journal of the Operational Research Society, 53(2), 211–221. Yun, Y., Yoon, M., & Nakayama, H. (2009). Multi-objective optimization based on metamodeling by using support vector regression. Optimization and Engineering, 10(2), 167–181. Zitzler, E., Deb, K., & Thiele, L. (2000). Comparison of multiobjective evolutionary algorithms: Empirical results. Evolutionary Computation, 8(2), 173–195. JOURNAL OF THE OPERATIONAL RESEARCH SOCIETY 913 Table A1. Parameter settings for SURROGATE-ASF. Initial number of sample points within each hyperbox n+1 Number of reference sample points within each sub-region k=2 30 (q = 29) k=3 66 (q = 10) k=4 84 (q = 6) k=5 125 (q = 5) Table A2. Parameter settings for MATSuMoTo. Metamodeling technique Cubic RBF Sampling technique Latin hypercube sampling New candidate points to be evaluated with the computationally expensive functions adding random perturbations to the best point found so far Initial number of sample points 2(n+1) Table A3. Parameter settings for DIRECT. Maximum number of function evaluations To solve problem (3) differs for each reference point and benchmark problem To solve problem (8) 500 Maximum number of iterations 500 500 Maximum number of rectangle divisions 1000 1000 Appendix 1. 914 M. TABATABAEI ET AL.