The Multiple-Update-Infill Sampling Method Using Minimum Energy Design for Sequential Surrogate Modeling

Abstract

Computer experiments are widely used to evaluate the performance and reliability of engineering systems with the lowest possible time and cost. Sometimes, a high-fidelity model is required to ensure predictive accuracy; this becomes computationally intensive when many computational analyses are required (for example, inverse analysis or uncertainty analysis). In this context, a surrogate model can play a valuable role in addressing computational issues. Surrogate models are fast approximations of high-fidelity models. One efficient way for surrogate modeling is the sequential sampling (SS) method. The SS method sequentially adds samples to refine the surrogate model. This paper proposes a multiple-update-infill sampling method using a minimum energy design to improve the global quality of the surrogate model. The minimum energy design was recently developed for global optimization to find multiple optima. The proposed method was evaluated with other multiple-update-infill sampling methods in terms of convergence, accuracy, sampling efficiency, and computational cost.

1. Introduction

The traditional development of an engineering system requires repeated experiments involving evaluation of variability in material properties and external conditions (for example, loading and excitation). Consequently, such development is a time-consuming and cost-expensive development process undertaken with limited resources, and some experiments are infeasible due to technical limitations (for example, extreme loadings such as earthquakes and collisions). In this regard, computer experiments, using the lowest amount of time and the lowest cost possible, are widely used to evaluate the performance reliability of engineering systems [1]. The underlying idea of a computer experiment is to simulate the performance and reliability of the engineering systems using computational models instead of time-consuming experiments. In order to ensure the predictive accuracy of computational models, a high-fidelity model is commonly required. High-fidelity models are typically time-consuming and computationally intensive. In addition, many computational analyses are required for some applications (that is, uncertainty or inverse analysis). As a result, applications using high-fidelity models become computationally expensive and infeasible under limited resources.

A surrogate model can play a valuable role in addressing these computational issues [2,3,4,5]. The surrogate model is constructed by fitting the input–output relations of a high-fidelity model, thus constructing a fast approximation of the high-fidelity model. Once the surrogate model is accurately constructed, the high-fidelity model is replaced by this constructed surrogate model to produce outputs at certain inputs. Constructing an accurate surrogate model needs adequate input–output samples (training samples) that can properly capture the characteristics of the response surface. Conventionally, training samples are generated by space-filling designs including Latin hypercube designs (LHDs) [6], orthogonal LHDs [7], uniform designs [8], and generalized LHDs [9]. A surrogate model is constructed at once using training samples (one-stage method).

Another way of constructing a surrogate model is the sequential sampling (SS) method. The SS method is an efficient way to construct a surrogate model. If the number of model analyses is not limited and a minimal number of training samples (unknown) is desired, the SS method is preferable for space-filling designs (that is, the one-stage method). The SS method typically involves two stages. Firstly, the surrogate model is constructed using initial training samples from space-filling designs. The second stage sequentially selects an additional sample (infill sample) and places it among the current training samples to update the surrogate model [10,11,12]. For the SS method, infill criteria have been developed to update the surrogate model by learning information about the response surface. The SS method can involve single-update-infill sampling [5,12,13,14] and multiple-update-infill sampling [10,15]. Single-update-infill sampling consists in the selection of only a single infill sample at each infill stage, whereas the multiple-update-infill sampling method consists in the selection of multiple infill samples to accelerate the improvement of global surrogate modeling.

This paper proposes a multiple-update-infill sampling method using a minimum energy design (MED) to improve the global quality of the surrogate model. The MED was recently developed for global optimization to find multiple optima [12]. The feasibility of the multiple-update-infill sampling method using an MED was evaluated using a mathematical test function (the Friedman model), a toy test function (the borehole model), and an engineering application (the finite element model of the 2D beam frame structure). The proposed method was compared to the multiple-update-infill sampling method using mean squared errors (MSEs) [14] and expected improvement (EI) [16]. Based on the comparison, we drew conclusions in terms of convergence, accuracy, sampling efficiency, and computational cost for global surrogate modeling. The multiple-update-infill sampling method using an MED seems to be promising and efficient with a minimum number of final training samples for general cases, except for response surfaces with noisy components or strong non-stationarity.

The remainder of this paper is organized as follows. Section 2 provides a brief introduction of the kriging model and space-filling design. Section 3 presents the concept of the SS method. Thereafter, three infill criteria and their multiple-update-infill sampling methods are described. In Section 4, three case studies (the Friedman model, the borehole model, and the finite element model of the frame structure) are considered to evaluate the performance of the multiple-update-infill sampling method. Finally, Section 5 provides conclusions from the case studies. Hereafter, the boldface letters indicate vectors or matrices. In the remaining of this paper, the bold format denotes a vector or matrix.

2. The Surrogate Model for Computer Experiments

The surrogate model is a fast approximation of the computational model. It is also known as the response surface, the meta-model, and the emulator. The procedure of surrogate modeling is shown in Figure 1. Firstly, input variables and target outputs are defined. Then, training samples (input variables) are generated by the design of experiments (for example, space-filling design). Based on the training samples, computational model analysis is performed to obtain target outputs. A mathematical model is used to construct the surrogate model by fitting the input–output relationships (response surface). Surrogate models can be categorized by regression-based models (for example, polynomial functions) and interpolation-based models (for example, radial basis functions and kriging). In this study, the response surfaces of interest are confined to those from the deterministic computation model. Stated differently, the response surfaces of interest do not contain noisy components (that is, small oscillations with a high-frequency component), discontinuity, or strong non-stationarity. In view of the response surfaces of interest and deterministic computation model, the kriging model was chosen to represent the interpolation-based model.

2.1. The Kriging Model

The kriging model was originally introduced in Geostatistics to estimate the properties of the most likely distribution of gold based on a given set of sampled sites [17]. The kriging model is also known as the Gaussian process model. The kriging model is an interpolation method for which the unknown predictions are modeled by the Gaussian process, which is expressed by a prior covariance. Stated differently, the unknown predictions ($\widehat{y}\left({\mathit{x}}_{\mathit{new}}\right)$) are modeled as a realization of the Gaussian process with the covariance conditional on the training samples. The kriging model can be expressed as

$$\widehat{y}\left(\mathit{x}\right)=\mathit{f}\left(\mathit{x}\right)+\mathit{Z}\left(\mathit{x}\right)\mathbf{with}\mathit{Z}\left(\mathit{x}\right)~GP\left(\mathbf{0},\mathit{C}\left(\mathit{x},{\mathit{x}}^{\prime }\right)\right)$$

(1)

where $\mathit{x}\in R^k$ is a sample with $k$ input variables; $\mathit{f}\left(\mathit{x}\right)$ can be constant or polynomial function to approximate a global trend over the entire parameter space; $\mathit{Z}\left(\mathit{x}\right)$ is a stochastic component to represent a local spatial deviation and assumed to be a stationary Gaussian process with a zero-mean and covariance ($\mathit{C}\left(\mathbf{·},\mathbf{·}\right)$) as

$$\mathit{C}\left(\mathit{Z}\left({\mathit{x}}_{\mathit{i}}\right),\mathit{Z}\left({\mathit{x}}_{\mathit{j}}\right)\right)={\sigma }^2\psi \left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)$$

(2)

where ${\sigma }^2$ is the process variance; the subscript $i$ and $j$ indicate the $i$th and $j$th training sample, respectively; $\psi \left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)$ is a spatial correlation function to control the smoothness of the kriging model. A typical choice of $\psi \left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)$ is the $k$-dimensional Gaussian correlation function given by the following:

$$\psi \left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)=exp\left(-{\displaystyle \sum }_{p=1}^k{\theta }_p{\Vert x_i^{\left(p\right)}-x_j^{\left(p\right)}\Vert }^2\right)$$

(3)

where $\Vert x_i^{\left(p\right)}-x_j^{\left(p\right)}\Vert {}^2$ is the Euclidean distance measure between two samples. $\mathit{\theta }=\left[{\theta }_1,\cdots ,{\theta }_p\right]$ is an unknown correlation parameter vector to scale the correlation length in each input. The correlation length also reflects the relative significance of each input [2].

Different types of the kriging model can be modeled as follows: (1) the ordinary kriging model assumes a constant mean only ($\mathit{f}\left(\mathit{x}\right)=\mu \left(\mathit{x}\right)$) over the neighborhood of $x$, and (2) the universal kriging model assumes a general polynomial trend model ($\mathit{f}\left(\mathit{x}\right)={{\displaystyle \sum }}_{i=0}^k{\beta }_i{x}_i$). Since the ordinary kriging model is used in this study, the ordinary kriging model is briefly described. Hereafter, the kriging model denotes the ordinary kriging model for simplicity.

Suppose that we have $n$ observed inputs ${\mathbf{X}}_{\mathit{n}}={\left[{\mathit{x}}_{\mathbf{1}},\cdots ,{\mathit{x}}_{\mathit{n}}\right]}^T$ and the corresponding outputs ${\mathbf{Y}}_{\mathit{n}}={\left[y\left({\mathit{x}}_{\mathbf{1}}\right),\cdots ,y\left({\mathit{x}}_{\mathit{n}}\right)\right]}^T$. To construct an ordinary kriging model using ${\mathbf{X}}_{\mathit{n}}$ and ${\mathbf{Y}}_{\mathit{n}}$, the unknown parameters of the kriging model ($\mu$, ${\sigma }^2$ and $\mathit{\theta }$) are estimated by maximum likelihood estimation. The log-likelihood function (Equation (4)) is expressed with these unknown parameters:

$$\ln \left(L\left(\mu ,{\sigma }^2,\mathit{\theta }|{\mathbf{X}}_{\mathit{n}},{\mathbf{Y}}_{\mathit{n}}\right)\right)\approx -\frac{n}{2}\ln \left({\sigma }^2\right)-\frac{1}{2}\ln \left|\mathsf{\Psi }\right|-\frac{{\left(Y_t-\mathbf{1}\mu \right)}^T{\mathsf{\Psi }}^{-1}\left({\mathbf{Y}}_t-\mathbf{1}\mu \right)}{2{\sigma }^2}$$

(4)

where $\mathsf{\Psi }$ is an $n\times n$ correlation matrix with the element $\psi \left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)$ and 1 is an n-by-1 unit vector. By taking the derivatives of Equation (4) with respect to $\mu$ and ${\sigma }^2$ and setting the derivatives to zero, their maximum likelihood estimators are derived as

$$\widehat{\mu }={\left({\mathbf{1}}^T{\mathsf{\Psi }}^{-\mathbf{1}}\mathbf{1}\right)}^{-1}{\mathbf{1}}^T{\mathsf{\Psi }}^{-1}{\mathbf{Y}}_t$$

(5)

$${\widehat{\sigma }}^2-\frac{{\left(Y_n-\mathbf{1}\widehat{\mu }\right)}^T{\mathsf{\Psi }}^{-1}\left({\mathbf{Y}}_{\mathit{n}}-\mathbf{1}\widehat{\mu }\right)}{n}.$$

(6)

Equations (5) and (6) show that the only parameter is the correlation parameters ($\mathit{\theta }$) in $\mathsf{\Psi }$. By substituting Equations (5) and (6) into Equation (4), the concentrated log-likelihood function is obtained as

$$\ln \left(L\left(\mathit{\theta }|{\mathbf{X}}_{\mathit{n}},{\mathbf{Y}}_{\mathit{n}}\right)\right)\approx -\frac{n}{2}\ln \left({\widehat{\sigma }}^2\right)-\frac{1}{2}\ln \left|\mathsf{\Psi }\right|.$$

(7)

The optimal correlation parameters ($\widehat{\mathit{\theta }}$) are estimated by maximizing Equation (7) as

$$\widehat{\mathit{\theta }}=\arg \underset{\mathit{\theta }}{\max }\ln \left(L\left(\mathit{\theta }|{\mathbf{X}}_{\mathit{n}},{\mathbf{Y}}_{\mathit{n}}\right)\right).$$

(8)

Once the optimal correlation parameters ($\widehat{\mathit{\theta }}$) are estimated, the unknown predictions ($\widehat{y}\left({\mathit{x}}_{\mathit{new}}\right)$) can be estimated by a weighted linear combination of the training samples (${\mathbf{X}}_{\mathit{n}}$ and ${\mathbf{Y}}_{\mathit{n}}$) as

$$\widehat{y}\left({\mathit{x}}_{\mathit{new}}\right)=\widehat{\mu }+{\widehat{\mathit{\psi }}}^T{\mathsf{\Psi }}^{-1}\left(Y_{\mathit{n}}-\mathbf{1}\widehat{\mu }\right)$$

(9)

where $\widehat{\mathit{\psi }}={\left[\psi \left({\mathit{x}}_{\mathit{new}},{\mathit{x}}_{\mathbf{1}}\right),\cdots ,\psi \left({\mathit{x}}_{\mathit{new}},{\mathit{x}}_n\right)\right]}^T$ represents the correlation vector between the unknown sample (${\mathit{x}}_{\mathit{new}}$) and ${\mathbf{X}}_{\mathit{n}}$. Equation (9) is the mean of the posterior Gaussian distribution. The prediction variance will be presented in Section 3.1. For more information about the kriging model, please refer to Jones [11] and Forrester et al. [2].

2.2. Space-Filling Design

Conventional surrogate modeling uses the training samples at once to fit a response surface. The training samples are generated by design of experiments (DOE). A space-filling design is the most widely used in computer experiments and it uniformly fills the samples over the input space by maximizing the distance between the samples as far apart as possible (that is, uniformity). The space-filling designs include Latin hypercube designs (LHDs) [6], maximum LHDs [18], orthogonal LHDs [7], uniform designs [8], generalized LHDs [9], and maximum projection designs [19], according to their own criterion for uniformity. Among the various space-filling designs, the maximum LHD (Mm LHD) seems to be the most commonly used space-filling design in practice thanks to its simplicity and availability in software packages [19]. Mm LHD searches an LHD by maximizing the minimum distance among the training samples, and this can be obtained by minimizing the following criterion (Morris–Mitchel criterion).

$$\underset{D}{\min }{\left\{{\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^n\frac{1}{d^k\left({\mathit{x}}_i,{\mathit{x}}_j\right)}\right\}}^{1/k}$$

(10)

where $d^k\left({\mathit{x}}_i,{\mathit{x}}_j\right)={\left({\displaystyle \sum }\text{}\left({\left({\mathit{x}}_i-{\mathit{x}}_j\right)}^k\right)\right)}^{1/k}$.

It is worth noting that space-filling designs were not originally designed for sequential surrogate modeling. Space-filling designs have been widely employed under a specific number of the training samples (due to a limited number of the model evaluation). In the sequential sampling method, the space-filling designs are typically used to generate initial training samples. Within the initial training samples, an initial surrogate model is constructed and used to select an infill sample for the sequential sampling method.

3. The Sequential Sampling Method

The concept of the sequential sampling (SS) method is illustrated in Figure 2, and this shows that the surrogate model converges to a true function by way of adding the infill samples to the current training samples. A core technique in the SS method is for the infill criteria to search the infill samples and refine the surrogate models for its own purpose (for example, the optimization or global quality of the surrogate model). There are numerous studies to evaluate and compare the infill criteria for sequential surrogate sampling [11,15,20,21,22]. These infill criteria include the mean squared error (MSE) [14], expected improvement (EI) [16], the probability of improvement (PoI) [23], statistical lower bound (SLB) [24], and so on. The PoI and SLB criteria are not suitable for the global surrogate modeling since they were developed to improve the prediction of the best value (pure exploitation). The MSE criterion can be used to add infill samples to the region that has high uncertainty in the predictions of the surrogate model (pure exploration), whereas the EI criterion can add infill samples for both exploitation and exploration of the response surface (balanced exploitation and exploration). Both criteria can be used as multiple-update-infill sampling methods to accelerate the surrogate modeling [10,15]. In this study, a minimum energy design (MED) is investigated to evaluate its feasibility for the SS method. The MED was recently developed to find multiple alternative optima rather than to accelerate the search for global optimum [12]. This section describes three infill criteria for the multiple-update-infill sampling method to improve the global quality of the surrogate model. These infill criteria are provided with their concepts in the following sequential order: (1) the mean squared error (MSE), (2) expected improvement (EI), and (3) the minimum energy design (MED).

3.1. Infill Criterion I: Mean Squared Error

The mean squared error (MSE) criterion is the prediction variance (${\widehat{s}}^2\left(\mathit{x}\right)$) of the kriging model. The MSE is derived from the curvature of the augmented log-likelihood function as [2,11,13,14]

$${\widehat{s}}^2\left(\mathit{x}\right)={\widehat{\sigma }}^2\left(\mathbf{1}-{\widehat{\mathit{\psi }}}^T{\mathsf{\Psi }}^{-1}\widehat{\mathit{\psi }}\right).$$

(11)

As shown in Figure 3, the predictions of the unknown samples can be represented as the Gaussian process with a mean ($\widehat{y}\left({\mathit{x}}_{\mathit{new}}\right)$) in Equation (9) and prediction variance (${\widehat{s}}^2\left(\mathit{x}\right)$) in Equation (11). The MSE can be used as an infill criterion to search the infill sample that has a large prediction uncertainty (larger ${\widehat{s}}^2\left(\mathit{x}\right)$). In Figure 3, there are two unknown samples ($x=3$ and $x=4$) with their prediction distributions (represented by purple lines). Based on Equation (11), ${\widehat{s}}^2\left(3\right)$ is much smaller than ${\widehat{s}}^2\left(4\right)$. This indicates that the bigger MSE is obtained far apart from the training samples. Stated differently, The MSE can be used as a measure of the sparseness of the input space. The infill sample of the MSE (${\mathit{x}}_{MSE}^{n+1}$) is chosen by maximizing the MSE as

$${\mathit{x}}_{MSE}^{n+1}=\arg \underset{\mathit{x}}{\max }{\widehat{s}}^2\left(\mathit{x}\right)$$

(12)

3.2. Infill Criterion II: Expected Improvement

The expected improvement (EI) is an infill criterion to compute how much improvement of the current kriging model can be achieved if an infill sample is added. Let $Y\left(\mathit{x}\right)$ be a random variable representing the uncertainty in the prediction of sample $\mathit{x}$. The uncertainty in the prediction is represented by the Gaussian distribution with a mean ($\widehat{y}\left(\mathit{x}\right)$) and prediction variance (${\widehat{s}}^2\left(\mathit{x}\right)$). If the current best value is the minimum value ($y_{min}$) thus far, then the improvement on the minimum ($\mathit{I}\left(\mathit{x}\right)$) can be defined as

$$\mathit{I}\left(\mathit{x}\right)=y_{min}-Y\left(\mathit{x}\right).$$

(13)

Since $Y\left(\mathit{x}\right)$ is the Gaussian distribution, Equation (13) can be integrated to compute the expectation of $\mathit{I}\left(\mathit{x}\right)$. Then, EI on minimum ($E{I}_{min}\left(\mathit{x}\right)$) can be defined as

$$E{I}_{min}\left(\mathit{x}\right)=\left(y_{min}-\widehat{y}\left(\mathit{x}\right)\right)\left[\frac{1}{2}+\frac{1}{2}\mathrm{erf}\left(\frac{y_{min}-\widehat{y}\left(\mathit{x}\right)}{\widehat{s}\left(\mathit{x}\right)\sqrt{2}}\right)\right]+\widehat{s}\left(\mathit{x}\right)\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{{\left(y_{min}-\widehat{y}\left(\mathit{x}\right)\right)}^2}{2{\widehat{s}}^2\left(\mathit{x}\right)}\right]$$

(14)

where $\mathrm{erf}(·)$ denotes the error function for the cumulative Gaussian distribution of the Kriging model at sample $\mathit{x}$. The first term in Equation (14) tends to be bigger when the mean ($\widehat{y}\left(\mathit{x}\right)$) is likely to be better than the current minimum (that is, exploitation). On the other hand, the second term in Equation (14) plays a role in exploration by getting bigger with higher prediction variance values (${\widehat{s}}^2\left(\mathit{x}\right)$). In this context, $E{I}_{min}\left(\mathit{x}\right)$ searches an infill sample which improves both exploitation and exploration in the response surface. The infill sample of the $E{I}_{min}\left(\mathit{x}\right)$ is chosen as

$${\mathit{x}}_{E{I}_{min}}^{n+1}=\arg \underset{\mathit{x}}{\max }E{I}_{min}\left(\mathit{x}\right).$$

(15)

EI on maximum ($E{I}_{max}\left(\mathit{x}\right)$) can be defined by replacing $y_{min}-Y\left(\mathit{x}\right)$ with $Y\left(\mathit{x}\right)-y_{max}$. Then, $E{I}_{max}\left(\mathit{x}\right)$ can be defined as the following equation:

$$E{I}_{max}\left(\mathit{x}\right)=\left(\widehat{y}\left(\mathit{x}\right)-y_{max}\right)\left[\frac{1}{2}+\frac{1}{2}\mathrm{erf}\left(\frac{\widehat{y}\left(\mathit{x}\right)-y_{max}}{\widehat{s}\left(\mathit{x}\right)\sqrt{2}}\right)\right]+\widehat{s}\left(\mathit{x}\right)\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{{\left(\widehat{y}\left(\mathit{x}\right)-y_{max}\right)}^2}{2{\widehat{s}}^2\left(\mathit{x}\right)}\right].$$

(16)

The infill sample of the $E{I}_{max}\left(\mathit{x}\right)$ is chosen as

$${\mathit{x}}_{E{I}_{min}}^{n+1}=\arg \underset{\mathit{x}}{\max }E{I}_{max}\left(\mathit{x}\right).$$

(17)

Figure 4 shows the $E{I}_{min}\left(\mathit{x}\right)$ and $E{I}_{max}\left(\mathit{x}\right)$ from the kriging model constructed by the three training samples. The values of $E{I}_{min}\left(\mathit{x}\right)$ are becoming larger near $x=5$ since this region has the greatest prediction variance (${\widehat{s}}^2\left(\mathit{x}\right)$) and is near the best observed value ($\mathrm{y}\left(6\right)$) thus far. In other words, the region near $x=5$ maximizes both terms in Equation (14) (Figure 4a). A similar observation can be found in $E{I}_{max}\left(x\right)$ (Figure 4b).

3.3. Infill Criterion III: Minimum Energy Design

The minimum energy design (MED) was recently developed as a new space-filling design to explore unknown regions of particular interest [12]. The key idea of the MED originates from the laws of electrostatics. In the MED, samples in the input space are considered as charged particles inside a box. Let these charged particles be positively charged. By setting the same charge for all particles, they will repel each other and occupy the positions inside the box to minimize the total potential energy (Figure 5). Let $q\left({\mathit{x}}_i\right)$ be the charge of the $i$th particle. The potential energy between the $i$th particle and $j$th particle (${\mathit{E}}_{ij}$) is proportional to the distance and charge as

$${\mathit{E}}_{\mathit{i}\mathit{j}}=\frac{q\left({\mathit{x}}_{\mathit{i}}\right)q\left({\mathit{x}}_{\mathit{j}}\right)}{d^2\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)}$$

(18)

where $d^2\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)$ is the Euclidean distance between the $i$th particle and the $j$th particle. The charge function $q\left(\mathit{x}\right)$ will be described later since this function is differently defined depending on the desired value of the output (that is, the maximum or minimum). For $n$-charged particles, the total potential energy is defined by

$$\mathit{E}={\displaystyle \sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}-\mathbf{1}}{\displaystyle \sum }_{\mathit{j}=\mathit{i}+\mathbf{1}}^{\mathit{n}}\frac{q\left({\mathit{x}}_{\mathit{i}}\right)q\left({\mathit{x}}_{\mathit{j}}\right)}{d^2\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)}.$$

(19)

In the MED, the box represents the input space, the positions taken by the charged particles are the training samples, and the charge of each particle reflects the value of the target output. Unlike the conventional space-filling designs, the MED tries to place the training samples as far apart from each other in the region of interest by reflecting their values (that is, the numerator of Equation (19)).

Joseph et al. [12] proposed a sequential strategy for implementing an MED. Suppose that we have $n$ observed inputs ${\mathbf{X}}^{\left(n\right)}={\left[{\mathit{x}}_{\mathbf{1}},\cdots ,{\mathit{x}}_{\mathit{n}}\right]}^T$ in a $k$-dimensional input space and the corresponding outputs ${\mathbf{Y}}^{\left(n\right)}={\left[y\left({\mathit{x}}_{\mathbf{1}}\right),\cdots ,y\left({\mathit{x}}_{\mathit{n}}\right)\right]}^T$. The kriging model of the $n$th stage ($KR{G}^{\left(n\right)}\left(\mathit{x}\right)$) can be constructed using the training samples (${\mathbf{X}}^{\left(n\right)}$ and ${\mathbf{Y}}^{\left(n\right)}$). Since the value of the output can be negative, an auxiliary function ($\widehat{g}\left(\mathit{x}\right)$) should be used instead of directly using a value of the output ($\widehat{y}\left(\mathit{x}\right)$).

$$\widehat{g}\left(\mathit{x}\right)=\widehat{y}\left(\mathit{x}\right)-{\widehat{y}}_{min}^{\left(n\right)}$$

(20)

where $\widehat{y}\left(\mathit{x}\right)$ is the mean value from $KR{G}^{\left(n\right)}\left(\mathit{x}\right)$ for the sample $\mathit{x}$, and ${\widehat{y}}_{min}^{\left(n\right)}$ is the estimated minimum value from $KR{G}^{\left(n\right)}\left(\mathit{x}\right)$. Then, the $\left(n+1\right)$th infill sample (${\mathit{x}}_{MED}^{n+1}$) is obtained by

$${\mathit{x}}_{MED}^{n+1}=\arg \underset{\mathit{x}}{\min }{\displaystyle \sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}-\mathbf{1}}{\displaystyle \sum }_{\mathit{j}=\mathit{i}+\mathbf{1}}^{\mathit{n}}\frac{q\left({\mathit{x}}_{\mathit{i}}\right)q\left({\mathit{x}}_{\mathit{j}}\right)}{d^2\left({\mathit{x}}_{\mathit{i}},{\mathit{x}}_{\mathit{j}}\right)}+{\displaystyle \sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}\left(\frac{q\left(\mathit{x}\right)q\left({\mathit{x}}_{\mathit{j}}\right)}{d^2\left(\mathit{x},{\mathit{x}}_{\mathit{j}}\right)}\right).$$

(21)

Since the first term of Equation (21) is constant regardless of the sample $\mathit{x}$, this term can be ignored. Thus, Equation (21) is reformulated as

$${\mathit{x}}_{MED}^{n+1}=\arg \underset{\mathit{x}}{\min }{\displaystyle \sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}\left(\frac{q\left(\mathit{x}\right)q\left({\mathit{x}}_{\mathit{j}}\right)}{d^2\left(\mathit{x},{\mathit{x}}_{\mathit{j}}\right)}\right).$$

(22)

The charge function ($q\left(\mathit{x}\right)$) can be divided into two types according to the desired value of the output. When the smaller value of the output (the minimum) is the region of interest, the charge ($q\left({\mathit{x}}_{\mathit{i}}\right)$) is proportional to $\widehat{g}\left(\mathit{x}\right)$ as

$$q_{max}\left(\mathit{x}\right)=\widehat{g}{\left(\mathit{x}\right)}^{2k}.$$

(23)

For the large value of the output (the maximum) is described, the charge ($q\left({\mathit{x}}_{\mathit{i}}\right)$) is inversely proportional to $\widehat{g}\left(\mathit{x}\right)$ as

$$q_{max}\left(\mathit{x}\right)=\frac{1}{\widehat{g}{\left(\mathit{x}\right)}^{2k}}.$$

(24)

The infill samples of the $ME{D}_{min}\left(x\right)$ and $ME{D}_{max}\left(x\right)$ are obtained by Equations (25) and (26), respectively.

$${\mathit{x}}_{ME{D}_{min}}^{n+1}=\arg \underset{\mathit{x}}{\min }{\displaystyle \sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}\left(\frac{q_{min}\left(\mathit{x}\right)q_{min}\left({\mathit{x}}_{\mathit{j}}\right)}{d^2\left(\mathit{x},{\mathit{x}}_{\mathit{j}}\right)}\right).$$

(25)

$${\mathit{x}}_{ME{D}_{max}}^{n+1}=\arg \underset{\mathit{x}}{\min }{\displaystyle \sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}\left(\frac{q_{max}\left(\mathit{x}\right)q_{max}\left({\mathit{x}}_{\mathit{j}}\right)}{{\mathit{d}}^{\mathbf{2}}\left(\mathit{x},{\mathit{x}}_{\mathit{j}}\right)}\right).$$

(26)

Figure 6 shows the $ME{D}_{min}\left(\mathit{x}\right)$ and $ME{D}_{max}\left(\mathit{x}\right)$ from the kriging model constructed by the three training samples. Both infill criteria search the regions that are highly populated near the global and local optima (exploitation), and they also try to be as far away from each other as possible (exploration).

3.4. The Multiple-Update-Infill Sampling Method

The multiple-update-infill sampling method can be defined as a combination of the infill criteria. For efficient global surrogate modeling, the SS method should contain two conflicting parts [10,25,26]: (1) global exploration and (2) local exploitation. The global exploration discovers interesting regions that have not been visited before. By doing so, it can reduce the bias in sampling due to the incomplete exploration of the entire input space. On the contrary, the local exploitation plays a key role in finding the regions of particular interest (maximum and minimum values or large prediction errors). Our previous study proposed the multiple-update-infill sampling method with three infill criteria ($MSE$ + $E{I}_{min}$ + $E{I}_{max}$) to contain both the global exploration and local exploitation for global surrogate modeling [10]. The rationale of the previous multiple-update-infill sampling method is as follows: (1) $MSE$ can select infill samples from the regions where training samples are spares, and (2) $E{I}_{min}$ and $E{I}_{max}$ can explore and exploit the potential regions near maximum and minimum to update the surrogate model.

The proposed method uses the MED to improve the capability for the balanced exploration and exploitation (that is, a similar role of $E{I}_{min}$ and $E{I}_{max}$). The potential problem in using the EI criteria is that the infill samples are sometimes clustered, and this results in the waste of the SS sampling and potential numerical instability. This is the motivation of considering the MED instead of the EI criteria in the proposed multiple-update-infill sampling method. The MED can select infill samples from the regions of particular interest (that is, maximum and minimum) by maximizing the inter-distance among the training samples. This feature of the MED is good for global surrogate modeling in order to avoid ill-conditioning due to clustering of the infill samples to some regions in the input space.

In this study, five multiple-update-infill sampling methods are considered to add multiple infill samples to ensure the global exploration and local exploitation. As a representative of the single-update-infill sampling method, the MSE is used to compare the multiple-update-infill sampling methods. The multiple-update-infill sampling methods are considered as follows: (1) $MSE$, (2) $E{I}_{max}+E{I}_{min}$, (3) $ME{D}_{max}+ME{D}_{min}$, (4) $MSE+E{I}_{max}+E{I}_{min}$, (5) $MSE+ME{D}_{max}+ME{D}_{min}$, and (6) $MSE+E{I}_{max}+E{I}_{min}+ME{D}_{max}+ME{D}_{min}$ (All).

4. Case Study on the Multiple-Update-Infill Sampling Methods

This section evaluates the sampling efficiency and performance of the SS methods based on six infill criteria combinations. Since the purpose of this study is to evaluate the feasibility of the proposed method for a minimal number of the training samples, three case studies with stationary responses surfaces are considered as follows: (1) one mathematical test function with 5 inputs, (2) one toy test function with 8 inputs, and (3) an engineering application with 13 inputs. Firstly, this section presents the preliminary description of the case study and four performance measures. Then, the three case studies are described with the landscapes of their response surface (that is, contour plots). Lastly, the results of the three case studies are discussed, and we draw conclusions in terms of convergence rate, accuracy, sampling efficiency, and computational cost to provide a guideline for less experienced users.

4.1. Preliminary Description of the Case Studies

Three case studies (that is, the Friedman model, the borehole model, and the 2D beam frame structure) were performed to evaluate the performances of the proposed SS methods. In surrogate modeling, it is important to normalize the input space (for example, ${\left[-11\right]}^p$ or ${\left[01\right]}^p$). By doing so, domination by a certain input (that is, with larger magnitude) can be avoided. In other words, normalization of the input space is employed to have an equal effect on the output of all case studies.

All multiple-update-infill sampling methods were performed under identical conditions 10 times. At each repetition, the different initial samples were generated, and these initial samples were identically used for all SS methods. To evaluate the performance of the SS method, the same validation samples (the number and their location) were used for each case study. The performance measures were calculated at each infill stage for all SS methods until the maximum number of the training samples was reached. For searching the infill sample and finding the minimum of the current kriging model in the MED (${\widehat{y}}_{min}^{\left(n\right)}$ of Equation (20)), the firefly algorithm was used as an effective global optimizer [27]. The flowchart of the multiple-update-infill sampling methods is tabulated in

Table 1. The flowchart of the multiple-update-infill sampling methods in the case studies.

Input:	Number of initial training samples $n_0$; Validation samples ${\mathit{X}}_{\mathit{v}\mathit{a}\mathit{l}}$; Number of final training samples $n_{final}$
1. Generate validation samples over input space
(1) Generate ${\mathit{X}}_{val}$ using Mm LHD (Equation (10))
(2) Compute outputs ${\mathit{Y}}_{val}$ using the computational model of ${\mathit{X}}_{val}$
2. Perform sequential sampling method
For $i=1,\cdots ,10$
(1) Set current infill stage to $n_s=1$
(2) Generate ${\mathit{X}}_i^{\left(n_s\right)}$ using Mm LHD (Equation (10)) with $n_0$
(3) Compute outputs ${\mathit{Y}}_i^{\left(n_s\right)}$ using the computational model of ${\mathit{X}}_i^{\left(n_s\right)}$
(4) Construct kriging model ($KR{G}_i^{\left(n_s\right)}$) using ${\mathit{X}}_i^{\left(n_s\right)}$ and ${\mathit{Y}}_i^{\left(n_s\right)}$
(5) Set # of training sample to $n_{train}=n_0$
Until ($n_{train}=n_{final}$)
(6) Predict outputs of ${\mathit{X}}_{val}$ (${\widehat{\mathit{Y}}}_{val}$) using $KR{G}_i^{\left(n_s\right)}$
(7) Evaluate and save performance measures using ${\mathit{Y}}_{val}$ and ${\widehat{\mathit{Y}}}_{val}$
(8) Search the infill sample(s) (${\mathit{X}}_{infill}$) based on $KR{G}_i^{\left(n_s\right)}$
(9) Compute outputs ${\mathit{Y}}_{infill}$ using computational model of ${\mathit{X}}_{infill}$
(10) Add ${\mathit{X}}_{infill}$ and ${\mathit{Y}}_{infill}$ to ${\mathit{X}}_i^{\left(n_s\right)}$ and ${\mathit{Y}}_i^{\left(n_s\right)}$
(11) Update kriging model ($KR{G}_i^{\left(n_s\right)}$) using ${\mathit{X}}_i^{\left(n_s\right)}$ and ${\mathit{Y}}_i^{\left(n_s\right)}$
(12) Set $n_{train}=n_{train}+N\left({\mathit{X}}_{infill}\right)$ and $n_s=n_s+1$
End
(13) Save final results: $KR{G}_i^{\left(final\right)}$, ${\mathit{X}}_i^{\left(final\right)}$, and history of the performance measure
End

Note: $\mathrm{N}\left({\mathit{X}}_{infill}\right)$ denotes the number of the infill sample(s).

.

This study considers four performance measures as follows: (1) root mean squared error ($RMSE$), (2) minimum distance ($d_{min}$), (3) pairwise correlation (${\rho }^2$), and (4) computational efficiency ($T_{norm}$). The $RMSE$ aggregates the magnitude of the residuals (that is, squared difference) into a single measure of the prediction accuracy. A smaller $RMSE$ indicates a higher prediction accuracy, as shown by the following equation:

$$RMSE=\sqrt{\frac{{\left({\mathit{Y}}_{val}-{\widehat{\mathit{Y}}}_{val}\right)}^T\left({\mathit{Y}}_{val}-{\widehat{\mathit{Y}}}_{val}\right)}{n_{val}}}$$

(27)

where ${\overline{Y}}_{val}$ is the average value of ${\mathit{Y}}_{val}$, and $n_{val}$ is the number of the validation samples. $d_{min}$ is the minimum value of the pairwise Euclidean distance between pairs of the final training samples (${\mathit{X}}_i^{\left(final\right)}$) as given in Equation (28). Larger values of $d_{min}$ indicate that the training samples are placed as far apart from each other as possible. In this context, a larger value of $d_{min}$ indicates efficient sampling performance in terms of the uniformity.

$$d_{min}=\underset{x_i,x_j\in {\mathit{X}}_i^{\left(final\right)}}{\min }{d}^2\left({\mathit{x}}_i,{\mathit{x}}_j\right).$$

(28)

${\rho }^2$ evaluates the goodness of the space-filling with respect to pairwise correlations [28]. A higher value of ${\rho }^2$ indicates that the projection capability in all subspace of the input variables is becoming poor. Therefore, the lower value of ${\rho }^2$ represents superior uniformity in the input space. ${\rho }^2$ is calculated as

$${\rho }^2=\frac{{{\displaystyle \sum }}_{i=2}^k{{\displaystyle \sum }}_{j=1}^{i-1}{\rho }_{ij}^2}{k\left(k-1\right)/2}$$

(29)

where ${\rho }_{ij}^2$ is the linear correlation between columns $i$ and $j$. The last performance measure is $T_{norm}$. Since single-update-infill sampling (that is, $MSE$) requires more infill stages than others, this can be used as the reference of the computational time. This measure is the normalized time by the average of the computational time based on the MSE:

$$T_{norm,i}=\frac{T_i}{1/N{{\displaystyle \sum }}_{k=1}^{10}{T}_k^{MSE}}$$

(30)

where $T_i$ is the computational time of the SS method, and $T_k^{MSE}$ denotes the computation time for the $k$th SS method based on MSE.

4.2. Case Study I: The Friedman Function

The Friedman function has been widely used to evaluate the performance of surrogate modeling (Friedman 1991). This test function has five input variables as given in Equation (31):

$$f\left(x\right)=10\sin \left(\pi x_1{x}_2\right)+20{\left(x_3-0.5\right)}^2+10x_4+5x_5.$$

(31)

Table 2. The nomenclature of the Friedman function.

Input Variable	Baseline	Lower Bound	Upper Bound
$x_1$	0.5	0	1
$x_2$	0.5	0	1
$x_3$	0.5	0	1
$x_4$	0.5	0	1
$x_5$	0.5	0	1

shows a baseline value and the lower and upper bound (range) of the input variables in the Friedman function.

The baseline values and ranges were used to generate a filled contour plot of the Friedman function. Each subplot in Figure 7 shows a contour of the Friedman function versus two of the five variables by fixing the remaining three variables to the baseline value. Figure 7 shows that all input variables sufficiently affect the output. There is a strong interaction effect due to the component of $10\sin \left(\pi x_1{x}_2\right)$.

4.3. Case Study II: The Borehole Model

A borehole model simulates water flow through a borehole drilled from the ground surface through two aquifers. This model is commonly used for testing various methods in computer experiments [29,30,31]. The input variables and their ranges are tabulated in

Table 3. The nomenclature of the borehole model [30].

Input Values		Baseline	Lower Bound	Upper Bound
$x_1$	Radius of borehole	0.1	0.05	0.15
$x_2$	Radius of influence	25,050	100	50,000
$x_3$	Transmissivity of upper aquifer	89,335	63,070	115,600
$x_4$	Potentiometric head of upper aquifer	1050	990	1110
$x_5$	Transmissivity of lower aquifer	89.55	63.1	116
$x_6$	Potentiometric head of lower aquifer	760	700	820
$x_7$	Length of borehole	1400	1120	1680
$x_8$	Hydraulic conductivity of borehole	8250	9855	12,045

. The output of this model is defined as

$$f\left(x\right)=\frac{2\pi x_3\left(x_4-x_6\right)}{\ln \left(x_2/x_1\right)\left(1+\frac{2x_7{x}_3}{\ln \left(x_2/x_1\right)x_1^2{x}_8}+\frac{x_3}{x_5}\right)}$$

(32)

where the output is water flow rate.

Figure 8 shows the contour plots for the response surfaces of the borehole model. Some of the input variables ($x_2,x_3$, and $x_5$) seem to be non-influential variables for the output. As a result, some of the contour plots are flat (for example, $x_2$ versus $x_3$). The contour plot between $x_1$ and $x_8$ exhibits the non-linear response surface, while the remaining contour plots show the linear response surface. Therefore, the borehole model seems as though the complex response surface is made in the eight-dimensional space.

4.4. Case Study III: The FE Model Based on a 3-Bay-5-Story Frame Structure

As an example of an engineering application, a three-bay–five-story frame structure is investigated. The frame structure was modeled using the Euler–Bernoulli beam with a lumped mass system and the finite element model of the frame structure was constructed in MATLAB [32]. As a target output, the displacement of the top floor is considered. For input variables, 13 model parameters and their lower/upper bounds were selected as tabulated in

Table 4. The five-story frame structure variables.

Input Variable		Baseline	Lower Bound	Upper Bound
P1	Load (KN)	133.45	53.38	333.63
P2		88.97	35.59	222.43
P3		71.18	28.47	177.95
E4	Young’s Modulus (MPa)	21,738	10,869	32,606
E5		23,796	11,898	35,659
I6	Moment of Inertia ($m^4$)	0.0081	0.0041	0.0122
I7		0.0115	0.0058	0.0173
I8		0.0232	0.0116	0.0348
I9		0.0259	0.0130	0.0389
A10	Sectional Area ($m^2$)	0.0312	0.0156	0.0468
A11		0.3716	0.1858	0.5574
A12		0.3725	0.1863	0.5588
A13		0.4181	0.2091	0.6272

. Figure 9 shows the frame structure with the input variables and the target output.

Figure 10 shows the contour plots of the response surfaces. The input variable P1 is the most influential variable among the input variables since it has larger values than the others and it acts from the third to fifth story of the frame. The response surfaces are smooth without any discontinuity. Therefore, the finite element model seems to be such that the smooth response surfaces are generated in the high-dimensional space.

4.5. Results and Discussion

All methods were performed with 30 initial samples for the case studies. To generate different initial samples, Mm LHD was performed 10 times. For the validation samples, 20,000 samples were generated using a quasi-random low-discrepancy sequence (that is, a sobol sequence). Based on the initial samples of each repetition, all SS methods were conducted until the number of the final samples was reached. The number of the final samples is tabulated in

Table 5. The preliminary setups for all SS methods.

Case Study		# Initial Samples	# Final Samples	# Validation Samples
1	Friedman Function	30	210	20,000
2	Borehole Model		240
3	FE Model Based on 3-bay-5- story Frame Structure		300

.

Based on the results of the case studies, the performances of all methods were evaluated in terms of convergence, accuracy, sampling efficiency, and computational cost. Firstly, we evaluated the convergence rate by the decreasing rate of the $RMSE$ values. A rapid decreasing rate of the $RMSE$ values indicated that the method accelerates the accuracy of global surrogate modeling. Secondly, the accuracy of the surrogate model was evaluated by the $RMSE$ values at the final samples. The lower $RMSE$ value indicated the superior global accuracy of the surrogate model. Next, the pairwise correlation and minimum distance were used to evaluate the measure of clustering the infill samples for sampling efficiency. Lastly, the computational cost was compared in terms of normalized time.

Convergence: The results of the case studies are shown in Figure 11, Figure 12, and Figure 13, respectively. For all case studies, $ME{D}_{max}+ME{D}_{min}$ and $MSE+ME{D}_{max}+ME{D}_{min}$ showed the best convergence performance with small variability. In Case Study 2, $MSE$ and $E{I}_{max}+E{I}_{min}$ were the worst SS methods in convergence (that is, with slow convergence and larger variability). In Case Study 3, $E{I}_{max}+E{I}_{min}$ was the worst performance in terms of convergence (that is, with divergence and larger variability). $MSE+E{I}_{max}+E{I}_{min}$ and “All” also provided convergence performances that were superior to those of $MSE$ and $E{I}_{max}+E{I}_{min}$. However, their convergence performance was not consistent for all case studies.

Accuracy: Accuracy was evaluated based on $RMSE$ values. In Figure 11b, Figure 12b, and Figure 13b, the last value of the $RMSE$ indicates the measure of the accuracy. Stated differently, lower $RMSE$ values indicate superior accuracy. For simplicity, the $RMSE$ values of the final surrogate model were used. $ME{D}_{max}+ME{D}_{min}$ and $MSE+ME{D}_{max}+ME{D}_{min}$ provided minimal $RMSE$ values for all case studies. In addition, their $RMSE$ values were less dispersed than the other SS methods as shown in Figure 11, Figure 12, and Figure 13. “All” also provided minimal $RMSE$ values but had a larger dispersion of $RMSE$ values than $ME{D}_{max}+ME{D}_{min}$ and $MSE+ME{D}_{max}+ME{D}_{min}$. The other SS methods were not as good as the abovementioned SS methods. In Case Study 3, $MSE$ showed excellent performance in terms of convergence; however, it showed poor performance in accuracy.

Sampling efficiency: The sampling efficiency is measured by the uniformity of the samples in the final surrogate model, as shown in Figure 14, Figure 15, and Figure 16. As expected, $MSE$ had the smallest values of ${\rho }^2$ with large $d_{min}$ values. Similar performances can be observed in $ME{D}_{max}+ME{D}_{min}$ and $MSE+ME{D}_{max}+ME{D}_{min}$. This indicates that the final samples are located over the input space as far apart from each other. On the contrary, the other SS methods have larger values of ${\rho }^2$ and smaller values of $d_{min}$. This implies that some of the training samples are closely placed. In other words, some samples are clustered to certain regions in the input space. This may result in the numerical instability in the construction of the surrogate model.

Computational efficiency: The computational cost of each SS method is shown in Figure 14, Figure 15, and Figure 16. As expected, the values of $T_{norm}$ became smaller in accordance with the increase of the infill samples at each infill stage. The optimization scheme (firefly algorithm) was implemented to search the infill samples, and numerous evaluations of the current surrogate model were required. The computational cost for the optimization scheme would be negligible in applications with computationally intensive models.

General discussion: $MSE+ME{D}_{max}+ME{D}_{min}$ showed the best performance in terms of convergence, accuracy, sampling, and computational efficiency, while $E{I}_{max}+E{I}_{min}$ showed the worst performance. Since $E{I}_{max}$ and $E{I}_{min}$ were originally developed to find the global optima, they tend to locate more infill samples in the regions near the optima. As the surrogate model becomes more accurate, infill criteria make the infill samples near the global optima for exploitation. On the contrary, $ME{D}_{max}$ and $ME{D}_{min}$ can allocate the infill samples that are far apart from the existing samples. Since more multiple infill criteria can save computational costs as discussed above, $MSE+ME{D}_{max}+ME{D}_{min}$ can accelerate the learning capability of the unknown response surface and can avoid the numerical instability caused by samples being too close to each other. It is worth noting that MED does not require a specific surrogate model. For example, many infill criteria (for example, MSE and EI) can be applied to only the kriging model. On the contrary, the multiple-update-infill sampling method using the MED criteria (that is $ME{D}_{max}+ME{D}_{min}$) can be applied to any surrogate models (penalized spline regression model [33], radial-basis function [34], and multivariate adaptive regression spline (MARS) [22]).

Limitations: The proposed method has a limitation for some applications. Since the kriging model in this study is based on the stationary correlation specification, it cannot suitably represent strong non-stationary response surfaces or noisy response surfaces (that is, high-frequency components due to approximation and truncation errors). Under such response surfaces, the kriging model with stationary correlation specification may result in numerical instabilities and poor predictive performances by violating the stationary assumption. To deal with such response surfaces, treed GPs [35] or a non-stationary covariance-based kriging method [36] are needed to ensure the predictive performance of the kriging model.

5. Conclusions

Many engineering applications sometimes require a larger number of time-consuming and computationally intensive models (that is, a high-fidelity model). Surrogate models are solutions that address this computational burden. The SS methods add the infill samples to the current training samples using infill criteria and refine the current surrogate model in a sequential framework. One advantage of SS methods is the generative learning capability on the underlying and unknown response surface. Numerous infill criteria have been developed for SS methods for their exploitation, exploration, and balanced exploitation and exploration. Multiple-update-infill sampling simultaneously uses multiple infill criteria to accelerate global surrogate modeling.

This study conducted case studies on several multiple-update-infill sampling methods for the global quality of the surrogate model. Three infill criteria were selected as follows: MSEs (pure exploration), EI, and an MED (balanced exploitation and exploration). Based on these infill criteria, six combinations of these infill criteria were considered. The performances of six different infill sampling methods have been compared through three case studies by increasing the dimensionality of the inputs.

Based on the results of the case studies, some observations have been made: Firstly, reducing the sparseness of the input space (that is $MSE$) may not result in the improvement on the global quality of the final surrogate model. In Case Study 2, $MSE$ consistently provides higher $RMSE$ values than other SS methods, except for $E{I}_{max}+E{I}_{min}$. This indicates that it is more efficient to make infill samples from the particular regions of interest (maximum and minimum) rather than from pure exploration ($MSE$). Secondly, some combinations of the infill criteria may result in the numerical instability because of samples that are too closely spaced. Lower values of $d_{min}$ with higher values of ${\rho }^2$ were observed in the three combinations, including $E{I}_{max}$ and $E{I}_{min}$. These infill criteria can increase the risk of making a singularity in the matrix manipulation in the surrogate model construction since they occasionally locate the infill samples near the optima to exploit the response surface. In this context, $MSE$, $ME{D}_{max}+ME{D}_{min}$, and $MSE+ME{D}_{max}+ME{D}_{min}$ are preferred in terms of numerical stability and sampling uniformity. Thirdly, multiple infill criteria can save computational time. As expected, more infill samples at each stage can further reduce computational time. In addition, one further benefit of the multiple infill criteria is a utilization of parallel computing in the high-fidelity model. Considering that the analysis of the high-fidelity model is the most intensive and expensive part in the SS methods, the availability of the parallel implementation can significantly reduce the computational time for the SS methods.

Based on the three case studies, the multiple-update-infill sampling methods using MED ($ME{D}_{max}+ME{D}_{min}$ and $MSE+ME{D}_{max}+ME{D}_{min}$) show the best performances in terms of convergence, accuracy, sampling, and computational efficiency for the global quality of the surrogate model. Implementation of both $ME{D}_{max}$ and $ME{D}_{min}$ can simultaneously generate infill samples with larger and smaller output values with larger inter-sample distances so that they can generate a dense set in the particular region by avoiding the risk of numerical instability. $MSE+ME{D}_{max}+ME{D}_{min}$ can save more computational time than $ME{D}_{max}+ME{D}_{min}$ with the parallel implementation of the high-fidelity model. It is worth noting that the MSE and EI are only available for the kriging model, but the MED can be applied to any surrogate model. In this context, the MED can be good infill criteria for global surrogate modeling to be automatically customized for various response surfaces ranging from smooth to complex ones.