Ensemble of Regression-Type and Interpolation-Type Metamodels

Abstract

Metamodels have become increasingly popular in the field of energy sources because of their significant advantages in reducing the computational cost of time-consuming tasks. Lacking the prior knowledge of actual physical systems, it may be difficult to find an appropriate metamodel in advance for a new task. A favorite way of overcoming this difficulty is to construct an ensemble metamodel by assembling two or more individual metamodels. Motivated by the existing works, a novel metamodeling approach for building the ensemble metamodels is proposed in this paper. By thoroughly exploring the characteristics of regression-type and interpolation-type metamodels, some useful information is extracted from the feedback of the regression-type metamodels to further improve the functional fitting capability of the ensemble metamodels. Four types of ensemble metamodels were constructed by choosing four individual metamodels. Common benchmark problems are chosen to compare the performance of the individual and ensemble metamodels. The results show that the proposed metamodeling approach reduces the risk of selecting the worst individual metamodel and improves the accuracy of the used individual metamodels.

1. Introduction

Metamodels, which are also referred to as surrogate models, are essentially approximate mathematical models of real physical systems. In the past decade, metamodels have become increasingly popular in the field of energy sources because of their significant advantages in reducing the computational cost of time-consuming tasks [1,2]. Melo et al. [3] pointed out that researchers in many countries are developing metamodels to estimate the energy performance of the building stock. Bornatico et al. [4] used a kind of metamodel to optimize energy systems, and found that the metamodel converged to the same solution at 150 times the speed of the fine model. Westermann and Evins [5] summarized and discussed recent studies on the application of metamodels in sustainable building design. Ferrero Bermejo et al. [6] reviewed and compared two typical metamodels, namely the artificial neural networks and the support vector machine, for energy forecasting and condition-based maintenance in PV plants.

Actually, a good metamodel mainly depends on its accuracy and generality for different design tasks. To enhance the performance of metamodels, researchers have carried out a lot of studies over the past few decades [7,8,9,10,11]. As a result, a large number of metamodels have been proposed, of which several types have gained wide acceptance in various applications. They are polynomial response surface (PRS) [12,13,14], support vector regression (SVR) [15,16,17], radial basis functions (RBF) [18,19], extended radial basis functions (E-RBF) [20], moving least squares (MLS) [21], artificial neural networks (ANN) [22,23], multivariate adaptive regressive splines (MARS) [24] and Kriging (KRG) [25,26]. These different metamodels give us more options for different tasks. However, lacking the prior knowledge of the actual physical systems, it is challenging to find a suitable metamodel in advance for a new task. In particular, the worst metamodel may be chosen for the task.

A simple way to overcome the difficulty is to build a series of metamodels based on a given training dataset at first, and then select the best one on the basis of some statistical techniques like the cross-validation method. Another favorite way is to construct an ensemble metamodel, which assembles two or more individual metamodels by introducing weight factors. The basic idea of such an ensemble metamodel can be traced back to 1990s [27,28], and currently it has become a research hotspot [8,29]. According to the characteristics of the weight factors, the techniques for building the ensemble metamodels can be mainly categorized into methods based on local errors, methods based on global errors, and methods based on regression.

In the first category, the weight factors $\left({\omega }_i={\omega }_i\left(\mathbf{x}\right)\right)$ are functions of design space, which are determined by the local errors of individual metamodels at the point of interest. Zerpa et al. [30] introduced a local weighted average model for the optimization of alkaline-surfactant-polymer flooding processes by using the prediction variances of three individual metamodels (PRS, KRG, and RBF). Sanchez, Pintos, and Queipo [31] proposed a general approach toward the ensemble of kernel-based models based on the local prediction variances. Acar [32] investigated the efficiency of methods based on the local errors, and developed a new approach to determine the weight factors by using the pointwise cross-validation errors instead of the prediction variances. Zhang, Chowdhury, and Messac [33] proposed a new metamodeling technique called adaptively hybrid functions, whose weight factors are determined based on the local measure of accuracy in the pertinent trust region. Lee and Choi [34] presented a new pointwise ensemble of metamodels, of which the weight factors are calculated by using the v nearest points cross-validation errors.

In the second category, the weight factors $({\omega }_i=C_i,\forall \mathbf{x})$ are constant values in the entire design space, which are determined by the global errors of individual metamodels. Goel et al. [35] studied a global weight factor selection approach based on the generalized mean square cross-validation errors (GMSE). Acar and Rais-Rohani [36] developed an accurate ensemble of metamodels by solving an optimization problem that minimizes GMSE or root mean square errors (RMSE). Viana, Haftka, and Steffen [37] obtained the optimal weight factors of the optimization problem by using the Lagrange multipliers. This method was also employed by Toal and Keane [38] to construct an ensemble of ordinary, universal, non-stationary and limit KRG models. Additionally, Acar [39] performed the simultaneous optimization of the weight factors and the shape parameters in the ensemble of RBFs.

It should be noted that in the first two categories the weight factors of individual metamodels are restricted to a positive range $\left({\omega }_i>0\right)$ and the sum of these factors is equal to 1 $\left({\sum }_{i=1}^M{\omega }_i=1\right)$. Since they are different from the first two categories, the techniques in the third category mainly use the regression methods (like least squares) to determine the weight factors. Accordingly, there is no longer any restriction on the weight factors, which may even have negative values. Polynkin and Toropov [40] introduced a novel mid-range metamodel assembly for the large-scale optimization problems, which is constructed based on the linear regression method. Ferreira and Serpa [41] developed an augmented least-square approach for creating the ensemble of metamodels, which can be extended to the efficient global optimization. Zhou and Jiang [42] constructed an ensemble of four individual metamodels (PRS, KRG, SVR, and RBF) from the view of the polynomial regression, and proposed a metamodel selection method on the basis of the stepwise regression to eliminate the redundant ones from the set of the candidate metamodels.

Motivated by these existing works, this paper proposes a different method for constructing the ensemble metamodels, which combines the advantages of regression-type and interpolation-type metamodels. The regression-type metamodels have better global trend fitting capacity than the interpolation-type metamodels, while the interpolation-type metamodels perform better than the regression-type metamodels in the vicinity of the sampling locations. By thoroughly exploring the characteristics of regression-type and interpolation-type metamodels, the proposed method could extract some useful information from the feedback of the regression-type metamodels to further improve the functional fitting capability of the ensemble metamodels.

2. Proposed Ensemble of Metamodels

2.1. Motivation and Basic Characteristics

The existing individual metamodels can be classified into regression-type and interpolation-type metamodels. The regression-type metamodels aim to fit the global trend of the underlying functions of the real physical systems in the entire design space, while the interpolation-type metamodels aim to achieve the local accuracy in the vicinity of the sampling locations. Accordingly, the regression-type metamodels can build smooth surfaces that pass across all the training points, while the interpolation-type metamodels can construct models that go through each training point. That is to say, for the regression-type metamodels there may be obvious deviations between the actual responses and the approximate responses at the sampling locations, while for the interpolation-type metamodels there is no deviation. These different characteristics make the two types of metamodels possess different advantages and limitations. For example: (i) the regression-type metamodels have better global trend fitting capacity than the interpolation-type metamodels, while (ii) the interpolation-type metamodels perform better than the regression-type metamodels in the vicinity of the sampling locations.

It should be noted that obtaining the training dataset required for constructing the metamodels may be time-consuming. Therefore, as much information as possible should be extracted from these data. However, for the regression-type metamodels, there are apparent deviations between the actual responses and the approximate responses at the sampling locations, from where some useful information may be still extracted to further improve the performance of these metamodels. Exploring the underlying knowledge of the training dataset and combining the characteristics of regression-type and interpolation-type metamodels, this paper proposes a novel metamodeling approach for the ensemble metamodels. The flowchart of the proposed metamodeling technique is shown in Figure 1, which involves four main steps as follows.

Step 1:

An appropriate design of experiment (DOE) should be first chosen to generate n sampling locations $({\mathbf{x}}^1,{\mathbf{x}}^2,\cdots ,{\mathbf{x}}^n)$, at where the actual responses $(y^1,y^2,\cdots ,y^n)$ are obtained by conducting experiments or simulations. By using the initial training dataset $({\mathbf{x}}^i,y^i)$$(i=1,\cdots ,n)$, a regression-type metamodel ${\widehat{y}}_1\left(\mathbf{x}\right)$ in Equation (1) is subsequently constructed to approximate the actual model $y\left(\mathbf{x}\right)$.

$${\widehat{y}}_1\left(\mathbf{x}\right)\approx y\left(\mathbf{x}\right),\mathbf{x}={(x_1,x_2,\cdots ,x_k)}^{\mathrm{T}}$$

(1)

where $\mathbf{x}$ denotes any point of interest.

Step 2:

We suppose that there is a deviation function $y_d\left(\mathbf{x}\right)$. It is obtained by subtracting the approximate model ${\widehat{y}}_1\left(\mathbf{x}\right)$ from the actual model $y\left(\mathbf{x}\right)$.

$$y_d\left(\mathbf{x}\right)=y\left(\mathbf{x}\right)-{\widehat{y}}_1\left(\mathbf{x}\right)$$

(2)

Some useful information may be still extracted from the deviation function $y_d\left(\mathbf{x}\right)$. To approximate the deviation function, the training dataset should be updated. In detail, this paper first uses the established regression-type metamodel in Equation (1) to predict the approximate responses $({\widehat{y}}_1^1,{\widehat{y}}_1^2,\cdots ,{\widehat{y}}_1^n)$ at the initial sampling locations. Subsequently, the deviations $(y_d^1,y_d^2,\cdots ,y_d^n)$ between the actual responses and approximate responses at these locations are calculated as the updated training dataset.

$$\begin{array}{cc}\hfill & \left\{({\mathbf{x}}^1,y_d^1),({\mathbf{x}}^2,y_d^2),\cdots ,({\mathbf{x}}^n,y_d^n)\right\}=\hfill \\ \hfill & \left\{({\mathbf{x}}^1,y^1-{\widehat{y}}_1^1),({\mathbf{x}}^2,y^2-{\widehat{y}}_1^2),\cdots ,({\mathbf{x}}^n,y^n-{\widehat{y}}_1^n)\right\}\hfill \end{array}$$

(3)

Step 3:

By using the updated training dataset in Equation (3), an interpolation-type metamodel ${\widehat{y}}_2\left(\mathbf{x}\right)$ in Equation (4) is constructed to approximate the deviation function $y_d\left(\mathbf{x}\right)$.

$${\widehat{y}}_2\left(\mathbf{x}\right)\approx y_d\left(\mathbf{x}\right)$$

(4)

Step 4:

Finally, the ensemble metamodel ${\widehat{y}}_{ens}\left(\mathbf{x}\right)$ in Equation (5) is constructed by adding the established regression-type metamodel ${\widehat{y}}_1\left(\mathbf{x}\right)$ and interpolation-type metamodel ${\widehat{y}}_2\left(\mathbf{x}\right)$ together. By using Equations (1), (4) and (5), the established ensemble metamodel ${\widehat{y}}_{ens}\left(\mathbf{x}\right)$ can be used to predict the response at any point of interest in the entire design space.

$${\widehat{y}}_{ens}\left(\mathbf{x}\right)={\widehat{y}}_1\left(\mathbf{x}\right)+{\widehat{y}}_2\left(\mathbf{x}\right)\approx {\widehat{y}}_1\left(\mathbf{x}\right)+y_d\left(\mathbf{x}\right)\approx y\left(\mathbf{x}\right)$$

(5)

2.2. Detailed Modeling Process

To clearly illustrate the proposed metamodeling technique, this paper selects two common regression-type metamodels (PRS and SVR) and two popular interpolation-type metamodels, namely RBFM (RBF with multiquadric-form basis function) and RBFI (RBF with inverse multiquadric-form basis function). Accordingly, four types of ensemble metamodels can be obtained, which are PrsRbfm (Ensemble Scheme 1, ensemble of PRS and RBFM), PrsRbfi (Ensemble Scheme 2, ensemble of PRS and RBFI), SvrRbfm (Ensemble Scheme 3, ensemble of SVR and RBFM) and SvrRbfi (Ensemble Scheme 4, ensemble of SVR and RBFI). The detailed modeling processes of these involved metamodels are introduced as follows.

2.2.1. Step 1: Construction of Regression-Type Metamodels

PRS is a general designation of a series of polynomial regression functions, of which the most popular one is the second-order polynomial model. This paper adopts the second-order polynomial model ${\widehat{y}}_{1,prs}\left(\mathbf{x}\right)$, which can be written as

$$\begin{array}{cc}\hfill {\widehat{y}}_{1,prs}\left(\mathbf{x}\right)& ={\mathbf{z}}^{\mathrm{T}}\mathbf{\beta }={\beta }_0+\sum _{i=1}^k{\beta }_i{x}_i+\sum _{i=1}^k\sum _{j=i}^k{\beta }_{\frac{2j+i(2k-i+1)}{2}}{x}_i{x}_j\hfill \end{array}$$

(6)

where $\mathbf{\beta }={({\beta }_0,{\beta }_1,\cdots ,{\beta }_{\frac{k^2+3k}{2}})}^{\mathrm{T}}$ denotes a coefficient vector, $\mathbf{z}=(1,x_1,x_2,\cdots ,$$x_{k-1}{x}_k,x_k{x}_k{)}^{\mathrm{T}}$ denotes a polynomial basis-function vector.

To estimate $\mathbf{\beta }$, the regression problem in Equation (6) can be transformed as follows by using the initial training dataset.

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}{y}^1\\ y^2\\ ⋮\\ y^n\end{array}\right]=\left[\begin{array}{c}{y}_{d,prs}^1\\ y_{d,prs}^2\\ ⋮\\ y_{d,prs}^n\end{array}\right]+\left[\begin{array}{ccccccccc}1& x_1^1& \cdots & x_k^1& \cdots & x_1^1{x}_1^1& \cdots & x_{k-1}^1{x}_k^1& x_k^1{x}_k^1\\ 1& x_1^2& \cdots & x_k^2& \cdots & x_1^2{x}_1^2& \cdots & x_{k-1}^2{x}_k^2& x_k^2{x}_k^2\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮\\ 1& x_1^n& \cdots & x_k^n& \cdots & x_1^n{x}_1^n& \cdots & x_{k-1}^n{x}_k^n& x_k^n{x}_k^n\end{array}\right]\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_{\frac{k^2+3k}{2}}\end{array}\right]\hfill \end{array}$$

(7)

where ${\mathbf{y}}_{d,prs}={(y_{d,prs}^1,y_{d,prs}^2,\cdots ,y_{d,prs}^n)}^{\mathrm{T}}$ denotes the deviation vector.

Equation (7) can be also expressed as

$$\begin{array}{cc}\hfill \mathbf{y}& =\mathbf{X}\mathbf{\beta }+{\mathbf{y}}_{d,prs}\hfill \end{array}$$

(8)

According to the least squares method, $\mathbf{\beta }$ can be calculated as follows.

$$\mathbf{\beta }={\left({\mathbf{X}}^{\mathrm{T}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}$$

(9)

SVR is a regression function ${\widehat{y}}_{1,svr}\left(\mathbf{x}\right)$ in the high-dimensional space, as shown in Equation (10).

$${\widehat{y}}_{1,svr}\left(\mathbf{x}\right)={\mathbf{\omega }}^{\mathrm{T}}\mathbf{\psi }\left(\mathbf{x}\right)+b$$

(10)

where $\mathbf{\omega }$ denotes the weight vector, $\mathbf{\psi }\left(\mathbf{x}\right)$ denotes the mapping function, and b denotes the bias.

To estimate $\mathbf{\omega }$ and b, the regression problem in Equation (10) can be transformed as an optimization problem in Equation (11) by introducing $ϵ$-insensitive loss function.

$$\begin{array}{cccc}\hfill & \min \hfill & \hfill & \frac{1}{2}{\left|\right|\mathbf{\omega }\left|\right|}^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \left\{\begin{array}{c}{\mathbf{\omega }}^{\mathrm{T}}\mathbf{\psi }\left({\mathbf{x}}^i\right)+b-y^i\le ϵ\hfill \\ y^i-{\mathbf{\omega }}^{\mathrm{T}}\mathbf{\psi }\left({\mathbf{x}}^i\right)-b\le ϵ\hfill \\ i=1,\cdots ,n\hfill \end{array}\right.\hfill \end{array}$$

(11)

To solve Equation (11), the regularization parameter, $C\left(>0\right)$, and the slack variables, ${\xi }^{+\left(i\right)}$ and ${\xi }^{-\left(i\right)}$, are introduced. In addition, Equation (12) can be obtained

$$\begin{array}{cccc}\hfill & \min \hfill & \hfill & \frac{1}{2}{\left|\right|\mathbf{\omega }\left|\right|}^2+C\sum _{i=1}^n({\xi }^{+\left(i\right)}+{\xi }^{-\left(i\right)})\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \left\{\begin{array}{c}{\mathbf{\omega }}^{\mathrm{T}}\mathbf{\psi }\left({\mathbf{x}}^i\right)+b-y^i\le ϵ+{\xi }^{+\left(i\right)}\hfill \\ y^i-{\mathbf{\omega }}^{\mathrm{T}}\mathbf{\psi }\left({\mathbf{x}}^i\right)-b\le ϵ+{\xi }^{-\left(i\right)}\hfill \\ {\xi }^{+\left(i\right)},{\xi }^{-\left(i\right)}\ge 0\hfill \\ i=1,\cdots ,n\hfill \end{array}\right.\hfill \end{array}$$

(12)

The Lagrange dual model of Equation (12) can be expressed as

$$\begin{array}{cccc}\hfill & \max \hfill & \hfill & \left\{\begin{array}{c}-\frac{1}{2}\sum _{i,j=1}^n({\alpha }^{+\left(i\right)}-{\alpha }^{-\left(i\right)})({\alpha }^{+\left(j\right)}-{\alpha }^{-\left(j\right)})\hfill \\ k〈{\mathbf{x}}^i,{\mathbf{x}}^j〉+\sum _{i=1}^n({\alpha }^{+\left(i\right)}-{\alpha }^{-\left(i\right)})y^i\hfill \\ -\sum _{i=1}^n({\alpha }^{+\left(i\right)}+{\alpha }^{-\left(i\right)})ϵ\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \left\{\begin{array}{c}\sum _{i=1}^n({\alpha }^{+\left(i\right)}-{\alpha }^{-\left(i\right)})=0\hfill \\ 0\le {\alpha }^{+\left(i\right)},{\alpha }^{-\left(i\right)}\le C\hfill \\ i=1,\cdots ,n\hfill \end{array}\right.\hfill \end{array}$$

(13)

where ${\alpha }^{+\left(i\right)}$ and ${\alpha }^{-\left(i\right)}$ denote the Lagrange multipliers, $k〈{\mathbf{x}}^i,{\mathbf{x}}^j〉$$=\mathbf{\psi }{\left({\mathbf{x}}^i\right)}^{\mathrm{T}}\mathbf{\psi }\left({\mathbf{x}}^j\right)$ denotes a kernel function, which has several different forms. This paper chooses the Gaussian kernel function, which can be expressed as

$$k〈\mathbf{x},{\mathbf{x}}^i〉=exp(-\gamma ||\mathbf{x}-{\mathbf{x}}^i{\left|\right|}^2)$$

(14)

According to Equation (13), ${\alpha }^{+\left(i\right)}$ and ${\alpha }^{-\left(i\right)}$ can be first obtained. According to KKT conditions [43], $\mathbf{\omega }$ and b can be then calculated.

2.2.2. Step 2: Update of Training Dataset

First, $\mathbf{\beta }$ calculated by Equation (9) can be used to substitute the one in Equation (6). Second, the approximate responses of established PRS $({\widehat{y}}_{1,prs}^1,{\widehat{y}}_{1,prs}^2,\cdots ,{\widehat{y}}_{1,prs}^n)$ at the initial sampling locations $({\mathbf{x}}^1,{\mathbf{x}}^2,\cdots ,{\mathbf{x}}^n)$ can be calculated according to Equation (6). Then, the updated training dataset of PRS can be expressed as

$$\begin{array}{cc}\hfill & \left\{({\mathbf{x}}^1,y_{d,prs}^1),({\mathbf{x}}^2,y_{d,prs}^2),\cdots ,({\mathbf{x}}^n,y_{d,prs}^n)\right\}=\hfill \\ \hfill & \left\{({\mathbf{x}}^1,y^1-{\widehat{y}}_{1,prs}^1),({\mathbf{x}}^2,y^2-{\widehat{y}}_{1,prs}^2),\cdots ,({\mathbf{x}}^n,y^n-{\widehat{y}}_{1,prs}^n)\right\}\hfill \end{array}$$

(15)

Similarly, according to Equation (10), the updated training dataset of SVR can be obtained and expressed as

$$\begin{array}{cc}\hfill & \left\{({\mathbf{x}}^1,y_{d,svr}^1),({\mathbf{x}}^2,y_{d,svr}^2),\cdots ,({\mathbf{x}}^n,y_{d,svr}^n)\right\}=\hfill \\ \hfill & \left\{({\mathbf{x}}^1,y^1-{\widehat{y}}_{1,svr}^1),({\mathbf{x}}^2,y^2-{\widehat{y}}_{1,svr}^2),\cdots ,({\mathbf{x}}^n,y^n-{\widehat{y}}_{1,svr}^n)\right\}\hfill \end{array}$$

(16)

2.2.3. Step 3: Construction of Interpolation-Type Metamodels

The general form of RBF can be expressed as

$${\widehat{y}}_{rbf}\left(\mathbf{x}\right)=\sum _{i=1}^n{\lambda }_i\varphi \left(\left|\right|\mathbf{x}-{\mathbf{x}}^i\left|\right|\right)$$

(17)

where ${\lambda }_i$ denotes an interpolation coefficient, $r=\left|\right|\mathbf{x}-{\mathbf{x}}^i\left|\right|=\sqrt{{(\mathbf{x}-{\mathbf{x}}^i)}^{\mathrm{T}}(\mathbf{x}-{\mathbf{x}}^i)}$ denotes the distance between points $\mathbf{x}$ and ${\mathbf{x}}^i$. $\varphi \left(r\right)$ denotes a radially symmetric basis function, which has several different forms, such as:

• Gaussian $\varphi \left(r\right)=e^{(-r^2/c^2)}$
• Multiquadric $\varphi \left(r\right)={(r^2+c^2)}^{\frac{1}{2}}$
• Inverse multiquadric $\varphi \left(r\right)={(r^2+c^2)}^{-\frac{1}{2}}$
• Thin plate spline $\varphi \left(r\right)=\left(r^2\right)log\left(r\right)$

The interpolation coefficient ${\lambda }_i$ can be calculated by using the given training dataset $({\mathbf{x}}^i,y^i)$$(i=1,\cdots ,n)$.

$$\mathbf{\lambda }={\mathbf{A}}^{-1}\mathbf{y}$$

(18)

where

$$\begin{array}{lll}\mathbf{\lambda }& =& {\left[\begin{array}{c}{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\end{array}\right]}^{\mathrm{T}}\\ \mathbf{A}& =& \left[\begin{array}{ccccc}& \varphi \left(\left|\right|{\mathbf{x}}^1-{\mathbf{x}}^1\left|\right|\right),& \varphi \left(\left|\right|{\mathbf{x}}^1-{\mathbf{x}}^2\left|\right|\right)& \cdots & \varphi \left(\left|\right|{\mathbf{x}}^1-{\mathbf{x}}^n\left|\right|\right)\\ & \varphi \left(\left|\right|{\mathbf{x}}^2-{\mathbf{x}}^1\left|\right|\right),& \varphi \left(\left|\right|{\mathbf{x}}^2-{\mathbf{x}}^2\left|\right|\right)& \cdots & \varphi \left(\left|\right|{\mathbf{x}}^2-{\mathbf{x}}^n\left|\right|\right)\\ & ⋮& ⋮& \ddots & ⋮\\ & \varphi \left(\left|\right|{\mathbf{x}}^n-{\mathbf{x}}^1\left|\right|\right),& \varphi \left(\left|\right|{\mathbf{x}}^n-{\mathbf{x}}^2\left|\right|\right)& \cdots & \varphi \left(\left|\right|{\mathbf{x}}^n-{\mathbf{x}}^n\left|\right|\right)\end{array}\right]\end{array}$$

After choosing the multiquadric-form basis function, RBFM (${\widehat{y}}_{rbfm}\left(\mathbf{x}\right)$) can be constructed to approximate the actual model $y\left(\mathbf{x}\right)$ by replacing ${\widehat{y}}_{rbf}\left(\mathbf{x}\right)$ and ${\lambda }_i$ in Equation (17) with ${\widehat{y}}_{rbfm}\left(\mathbf{x}\right)$ and ${\lambda }_{i,rbfm}$. The coefficient ${\lambda }_{i,rbfm}$ can be calculated based on Equation (18). Similarly, after choosing the inverse multiquadric-form basis function, RBFI (${\widehat{y}}_{rbfi}\left(\mathbf{x}\right)$) can be constructed to approximate the actual model $y\left(\mathbf{x}\right)$. The coefficient ${\lambda }_{i,rbfi}$ of ${\widehat{y}}_{rbfi}\left(\mathbf{x}\right)$ can be calculated based on Equation (18).

Additionally, by choosing the multiquadric-form basis function, a model ${\widehat{y}}_{2,rbfm1}\left(\mathbf{x}\right)$ can be constructed to approximate the deviation function of PRS ${\mathbf{y}}_{d,prs}$. By replacing the initial training dataset $({\mathbf{x}}^i,y^i)$$(i=1,\cdots ,n)$ with the updated training dataset of PRS $({\mathbf{x}}^i,y_{d,prs}^i)$$(i=1,\cdots ,n)$, the coefficient ${\lambda }_{i,2rbfm1}$ of ${\widehat{y}}_{2,rbfm1}\left(\mathbf{x}\right)$ can be calculated on the basis of Equation (18). Similarly, by choosing the inverse multiquadric-form basis function, a model ${\widehat{y}}_{2,rbfi1}\left(\mathbf{x}\right)$ can be constructed to approximate the deviation function of PRS ${\mathbf{y}}_{d,prs}$.

Finally, by choosing the multiquadric-form basis function, a model ${\widehat{y}}_{2,rbfm2}\left(\mathbf{x}\right)$ can be constructed to approximate the deviation function of SVR ${\mathbf{y}}_{d,svr}$. By choosing the interpolation-type metamodel, a model ${\widehat{y}}_{2,rbfi2}\left(\mathbf{x}\right)$ can be constructed to approximate the deviation function of SVR ${\mathbf{y}}_{d,svr}$.

2.2.4. Step 4: Construction of Ensemble Metamodels

By adding the established ${\widehat{y}}_{1,prs}\left(\mathbf{x}\right)$ and ${\widehat{y}}_{2,rbfm1}\left(\mathbf{x}\right)$ together, PrsRbfm (${\widehat{y}}_{prsrbfm}\left(\mathbf{x}\right)$) can be subsequently constructed as follows.

$${\widehat{y}}_{prsrbfm}\left(\mathbf{x}\right)={\widehat{y}}_{1,prs}\left(\mathbf{x}\right)+{\widehat{y}}_{2,rbfm1}\left(\mathbf{x}\right)$$

(19)

Being similar to PrsRbfm, PrsRbfi (${\widehat{y}}_{prsrbfi}\left(\mathbf{x}\right)$) can be constructed as follows.

$${\widehat{y}}_{prsrbfi}\left(\mathbf{x}\right)={\widehat{y}}_{1,prs}\left(\mathbf{x}\right)+{\widehat{y}}_{2,rbfi1}\left(\mathbf{x}\right)$$

(20)

SvrRbfm (${\widehat{y}}_{svrrbfm}\left(\mathbf{x}\right)$) can be constructed as follows.

$${\widehat{y}}_{svrrbfm}\left(\mathbf{x}\right)={\widehat{y}}_{1,svr}\left(\mathbf{x}\right)+{\widehat{y}}_{2,rbfm2}\left(\mathbf{x}\right)$$

(21)

SvrRbfi (${\widehat{y}}_{svrrbfi}\left(\mathbf{x}\right)$) can be constructed as follows.

$${\widehat{y}}_{svrrbfi}\left(\mathbf{x}\right)={\widehat{y}}_{1,svr}\left(\mathbf{x}\right)+{\widehat{y}}_{2,rbfi2}\left(\mathbf{x}\right)$$

(22)

The established ensemble metamodels, namely PrsRbfm, PrsRbfi, SvrRbfm, and SvrRbfi, can be used to predict the response at any point of interest in the entire design space by using Equations (19)–(22).

3. Numerical Experiments

3.1. Benchmark Problems

Referred to the website (http://www.sfu.ca/~ssurjano/index.html) and relevant literature [32], six common benchmark problems (BPs) are selected to compare the performance of the individual metamodels (PRS, SVR, RBFM, and RBFI) and the ensemble metamodels (PrsRbfm, PrsRbfi, SvrRbfm, and SvrRbfi).

BP1: 

Goldstein Price Function

$$\begin{array}{cc}\hfill f\left(\mathbf{x}\right)=& \left[1+{(x_1+x_2+1)}^2\times (19-14x_1+3x_1^2-14x_2+6x_1{x}_2+3x_2^2)\right]\times \hfill \\ & \left[30+{(2x_1-3x_2)}^2\times (18-32x_1+12x_1^2+48x_2-36x_1{x}_2+27x_2^2)\right]\hfill \end{array}$$

(23)

where $x_i\in [-2,2]$, for $i=1,2$.

BP2: 

Friedman Function

$$\begin{array}{cc}\hfill f\left(\mathbf{x}\right)=& 10sin\left(\pi x_1{x}_2\right)+20{\left(x_3-0.5\right)}^2+10x_4+5x_5\hfill \end{array}$$

(24)

where $x_i\in [0,1]$, for all $i=1,\cdots ,5$.

BP3: 

Power Sum Function

$$f\left(\mathbf{x}\right)=\sum _{j=1}^6{\left[\left(-\sum _{i=1}^6{x}_i^j\right)-36\right]}^2$$

(25)

where $x_i\in [0,6]$, for all $i=1,\cdots ,6$.

BP4: 

Rosenbrock Function

$$f\left(\mathbf{x}\right)=\sum _{i=1}^6\left[100{\left(x_{i+1}-x_i^2\right)}^2+{\left(x_i-1\right)}^2\right]$$

(26)

where $x_i\in [-5,10]$, for all $i=1,\cdots ,7$.

BP5: 

Zakharov Function

$$f\left(\mathbf{x}\right)=\sum _{i=1}^9{x}_i^2+{\left(\sum _{i=1}^{9}0.5i{x}_i\right)}^2+{\left(\sum _{i=1}^{9}0.5i{x}_i\right)}^4$$

(27)

where $x_i\in [-5,10]$, for all $i=1,\cdots ,9$.

BP6: 

Powell Function

$$\begin{array}{cc}\hfill f\left(\mathbf{x}\right)=& \sum _{i=1}^2\left[{\left(x_{4i-3}+10x_{4i-2}\right)}^2+5{\left(x_{4i-1}-x_{4i}\right)}^2+{\left(x_{4i-2}-2x_{4i-1}\right)}^4+10{\left(x_{4i-3}-x_{4i}\right)}^4\right]\hfill \end{array}$$

(28)

where $x_i\in [-4,5]$, for all $i=1,\cdots ,10$.

3.2. Numerical Setting

For all the benchmark problems, the MATLAB routine “lhsdesign” is used to generate training points and test points. Referred to Jin, Chen, and Simpson [44], $n=\frac{3(k+1)(k+2)}{2}$ training points are selected for a k-dimension problem. Moreover, as many test points as possible should be used in practice, since insufficient test points may increase the uncertainty of the results. This paper selects $n_{tst}=$ 20,000 test points for each benchmark problem. Since the DOE sampling scheme may have an obvious influence on the performance of the metamodels, 100 different training and test sets are selected for each problem. The detailed numerical settings for all the benchmark problems are listed in

Table 1. Detailed numerical settings for the benchmark problems.

Benchmark Problem	NO. of Variables	NO. of Training Points	NO. of Test Points	NO. of Training and Test Sets
BP1	2	18	20,000	100
BP2	5	63	20,000	100
BP3	6	84	20,000	100
BP4	7	108	20,000	100
BP5	9	165	20,000	100
BP6	10	198	20,000	100

. The shape parameters (c) of RBFM and RBFI are both selected as 1 by referring to relevant literature [34,45,46]. The parameters ($ϵ$, C, and $\gamma$) of SVR are selected by using the cross-validation method, which was introduced in detail in the published paper of the authors [47].

3.3. Performance Criteria

The root mean square error (RMSE) and the max absolute error (MAE) are selected as the performance criteria.

RMSE can be expressed as

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^{n_{tst}}{(y^i-{\widehat{y}}^i)}^2}{n_{tst}}}$$

(29)

where $n_{tst}$ denotes the number of test points.

MAE can be expressed as

$$MAE=max|y_i-\widehat{y_i}|,i=1,2,\cdots ,n_{tst}$$

(30)

4. Results and Discussion

4.1. RMSE

Figure 2 shows the boxplots of RMSE of the metamodels over 100 test sets for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points. It can be seen that: (1) for all the benchmark problems, the most accurate ensemble metamodels outperform the most accurate individual metamodels; (2) without exception, the least accurate individual metamodels perform worse than the least accurate ensemble metamodels; (3) for each benchmark problem, the performance differences among the four individual metamodels are greater than that among the four ensemble metamodels.

To provide a better comparison for these metamodels, the error values are normalized with respect to the most accurate individual metamodel for each benchmark problem.

Table 2. Normalized means of RMSE of the metamodels for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points.

BPs	BP1	BP2	BP3	BP4	BP5	BP6
PRS	1.280	1.866	1.113	1.000	1.000	1.000
SVR	1.224	1.000	1.262	1.149	1.006	1.133
RBFM	1.000	1.108	1.000	1.001	1.123	1.385
RBFI	1.133	1.261	1.536	2.175	2.073	2.166
PrsRbfm	0.929	1.043	0.981	0.957	0.889	0.989
PrsRbfi	0.977	1.173	1.062	0.990	0.985	0.994
SvrRbfm	0.968	0.922	1.039	1.006	0.778	1.080
SvrRbfi	1.010	0.937	1.176	1.102	0.910	1.109
Best & Best	−7.1%	−7.8%	−1.9%	−4.3%	−22.2%	−1.1%
Worst & Worst	−21.1%	−37.1%	−23.5%	−49.3%	−52.5%	−48.8%

shows the normalized means of RMSE of the metamodels for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points. The bold values in Table 2 are the most accurate individual/ensemble metamodels, the italic values are the least accurate individual/ensemble metamodels, the underlined values are the ensemble metamodels that perform better than all the individual metamodels, the “Best & Best” values denote the differences between the most accurate ensemble metamodels and individual metamodels, and the “Worst & Worst” values denote the differences between the least accurate ensemble metamodels and individual metamodels. From Table 2, it can be seen that: (1) compared with the most accurate individual metamodels, the means of RMSE of the most accurate ensemble metamodels are reduced, ranging from 1.1% to 22.2%; (2) compared with the least accurate individual metamodels, the means of RMSE of the least accurate ensemble metamodels are reduced, ranging from 21.1% to 52.5%; (3) except for BP3, more than two ensemble metamodels perform better than the most accurate individual metamodels; (4) for BP5, all the four ensemble metamodels perform better than the most accurate individual metamodel.

Table 3. Frequency of the accuracy ranking (using RMSE) of the metamodels for the six benchmark problems with $\frac{3(k+1)(k+2)}{2}$ training points.

Ranking	1st	2nd	3rd	4th	5th	6th	7th	8th
PRS	0	0	2	0	2	0	0	2
SVR	0	0	1	0	0	2	3	0
RBFM	0	1	0	2	1	0	2	0
RBFI	0	0	0	0	0	1	1	4
Total	0	1	3	2	3	3	6	6
PrsRbfm	4	1	0	1	0	0	0	0
PrsRbfi	0	2	1	2	0	1	0	0
SvrRbfm	2	1	1	1	1	0	0	0
SvrRbfi	0	1	1	0	2	2	0	0
Total	6	5	3	4	3	3	0	0

shows the frequency of the accuracy ranking (using RMSE) of the metamodels for the six benchmark problems with $\frac{3(k+1)(k+2)}{2}$ training points. It can be seen that: (1) the frequency of the ensemble metamodels that rank 1st or 2nd is 11, yet the frequency of the individual metamodels is only one; (2) the frequency of the individual metamodels that rank 7th or 8th is 12, yet the frequency of the ensemble metamodels is zero; (3) considered the frequency of the metamodels that rank the top/bottom two, all the ensemble metamodels have better performance than the individual metamodels; (4) PrsRbfm performs best among the four ensemble metamodels, followed by SvrRbfm, PrsRbfi, and SvrRbfi.

To clearly compare the accuracy of each ensemble metamodel with their corresponding individual metamodels, Figure 3 shows the normalized means of RMSE of each ensemble scheme for the six benchmark problems with $\frac{3(k+1)(k+2)}{2}$ training points. It can be seen that: (1) in Scheme 1, PrsRbfm ranks 1st among PRS, RBFM, and PrsRbfm for all the benchmark problems; (2) in Scheme 2, PrsRbfi ranks 1st for all the benchmark problems; (3) in Scheme 3, SvrRbfm ranks 1st for four benchmark problems and 2nd for two benchmark problems; although RBFM ranks 1st for two benchmark problems, it is the worst performer for three benchmark problems; (4) in Scheme 4, without exception, the accuracy of SvrRbfi outperforms that of SVR and RBFI.

Table 4. Normalized standard deviations of RMSE of the metamodels for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points.

BPs	BP1	BP2	BP3	BP4	BP5	BP6
PRS	1.000	1.000	1.646	1.566	1.001	1.072
SVR	2.067	1.542	1.121	10.663	6.188	6.844
RBFM	1.462	1.167	1.000	1.000	1.204	1.488
RBFI	1.660	1.100	1.235	1.117	1.000	1.000
PrsRbfm	1.397	1.185	1.304	1.592	0.744	0.995
PrsRbfi	1.423	1.193	1.493	1.577	0.953	1.055
SvrRbfm	1.581	1.143	1.509	2.021	1.427	3.714
SvrRbfi	1.696	1.199	1.123	7.041	3.991	5.542
Best & Best	39.7%	14.3%	12.3%	57.7%	−25.6%	−0.5%
Worst & Worst	−17.9%	−22.2%	−8.4%	−34.0%	−35.5%	−19.0%

shows the normalized standard deviations of RMSE of the metamodels for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points. It can be seen that: (1) compared with the most accurate individual metamodels, the standard deviations of RMSE of the most accurate ensemble metamodels are reduced for BP5 and BP6, yet the standard deviations are increased for the other four benchmark problems; (2) compared with the least accurate individual metamodels, the standard deviations of RMSE of the least accurate ensemble metamodels are reduced, ranging from 8.4% to 35.5%.

According to the above experimental results, we think the proposed metamodeling approach could reduce the risk of selecting the worst individual metamodel, and the constructed ensemble metamodels perform better than the used individual metamodels in terms of accuracy. In particular, PrsRbfm performs best among the four ensemble metamodels, followed by SvrRbfm, PrsRbfi, and SvrRbfi.

To provide an explicit explanation for the better performance of the proposed approach, a low-dimensional problem (BP1) and an ensemble scheme (ensemble of SVR and RBFM) are selected as examples. Figure 4 shows the contour plot of the actual function and the approximate functions of SVR, RBFM, and SvrRbfm. It can be seen that: (1) SVR has better global trend fitting capacity than RBFM, such as in the red box area; (2) RBFM performs better in the vicinity of the sampling locations, such as in the red ellipse region; (3) SvrRbfm combines the global trend of SVR and the local accuracy of RBFM, such as in the red box area and the red ellipse region.

Therefore, the reason for the better performance of the ensemble metamodels may be that the proposed metamodeling approach combines the advantages of the regression-type and interpolation-type metamodels. The actual model is regarded as the sum of a regression-type model and a deviation function. Some useful information is first extracted by the regression-type metamodel to capture the global trend of the actual model in the entire design space. Then, some other information is extracted from the deviations at the sampling locations by using the interpolation-type metamodel to achieve the local accuracy in the vicinity of sampling locations.

4.2. Effect of Performance Criteria

The choice of different performance criteria may influence the results of the metamodels. To reduce the source of uncertainty in the results as much as possible, the max absolute error (MAE) is selected as another performance criterion.

Figure 5 shows the boxplots of MAE of the metamodels over 100 test sets for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points.

Table 5. Normalized means of MAE of the metamodels for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points.

BPs	BP1	BP2	BP3	BP4	BP5	BP6
PRS	1.055	1.583	1.000	1.000	1.000	1.000
SVR	1.149	1.000	1.405	1.312	1.080	1.325
RBFM	1.000	1.148	1.165	1.174	1.385	1.708
RBFI	1.189	1.200	1.626	2.515	1.939	2.353
PrsRbfm	0.910	1.111	0.952	0.980	0.944	1.015
PrsRbfi	0.965	1.278	0.999	0.999	0.996	1.002
SvrRbfm	0.950	0.933	1.164	1.168	0.957	1.244
SvrRbfi	1.021	0.956	1.338	1.285	1.052	1.302
Best & Best	−9.0%	−6.7%	−4.8%	−2.0%	−5.6%	0.2%
Worst & Worst	−14.2%	−19.3%	−17.7%	−48.9%	−45.7%	−44.7%

shows the normalized means of MAE of the metamodels for each benchmark problem with $\frac{3(k+1)(k+2)}{2}$ training points. From Figure 5 and Table 5, it can be seen that: (1) for each benchmark problem, the performance differences among the four ensemble metamodels are less than that among the four individual metamodels; (2) except for BP6, more than two ensemble metamodels perform better than the most accurate individual metamodels; (3) compared with the most accurate individual metamodels, the means of MAE of the most accurate ensemble metamodels are reduced for five benchmark problems; (4) compared with the least accurate individual metamodels, the means of MAE of the least accurate ensemble metamodels are reduced, ranging from 14.2% to 48.9%.

Table 6. Frequency of the accuracy ranking (using MAE) of the metamodels for the six benchmark problems with $\frac{3(k+1)(k+2)}{2}$ training points.

Ranking	1st	2nd	3rd	4th	5th	6th	7th	8th
PRS	1	0	2	1	0	1	0	1
SVR	0	0	1	0	0	2	3	0
RBFM	0	0	0	1	3	0	2	0
RBFI	0	0	0	0	0	1	0	5
Total	1	0	3	2	3	4	5	6
PrsRbfm	4	0	1	1	0	0	0	0
PrsRbfi	0	3	2	0	0	0	1	0
SvrRbfm	1	2	0	3	0	0	0	0
SvrRbfi	0	1	0	0	3	2	0	0
Total	5	6	3	4	3	2	1	0

shows the frequency of the accuracy ranking (using MAE) of the metamodels for the six benchmark problems with $\frac{3(k+1)(k+2)}{2}$ training points. It can be seen that: (1) considered the frequency of the metamodels that rank the top/bottom two, PrsRbfm, PrsRbfi, and SvrRbfm outperform all the individual metamodels; (2) although SvrRbfi is a little worse than PRS, it still performs better than its corresponding individual metamodels (SVR and RBFI); (3) PrsRbfm is the best performer of the four ensemble metamodels, followed by SvrRbfm, PrsRbfi, and SvrRbfi.

In summary, the choice of the performance criteria influence the results slightly, but the conclusions obtained by the two criteria remain unchanged.

4.3. Effect of Sampling Densities

The choice of different sampling densities may also influence the results of the metamodels. To investigate the effect of the sampling densities, this paper selects another two schemes with different sampling densities, which are $n=\frac{5(k+1)(k+2)}{4}$ and $n=\frac{7(k+1)(k+2)}{4}$.

Table 7. Normalized means of RMSE of the metamodels for each benchmark problem with $\frac{7(k+1)(k+2)}{4}$ training points.

BPs	BP1	BP2	BP3	BP4	BP5	BP6
PRS	1.400	2.282	1.104	1.000	1.342	1.000
SVR	1.209	1.000	1.283	1.154	1.000	1.130
RBFM	1.000	1.186	1.000	1.012	1.496	1.371
RBFI	1.174	1.383	1.555	2.188	2.847	2.195
PrsRbfm	0.921	1.102	0.952	0.946	1.158	0.991
PrsRbfi	0.985	1.271	1.044	0.987	1.316	0.992
SvrRbfm	0.958	0.937	1.040	1.000	0.919	1.076
SvrRbfi	1.014	0.949	1.191	1.102	0.972	1.105
Best & Best	−7.9%	−6.3%	−4.8%	−5.4%	−8.1%	−0.9%
Worst & Worst	−27.6%	−44.3%	−23.4%	−49.6%	−53.8%	−49.7%

shows the normalized means of RMSE of the metamodels for each benchmark problem with $\frac{7(k+1)(k+2)}{4}$ training points. It can be seen that: (1) compared with the most accurate individual metamodels, the means of RMSE of the most accurate ensemble metamodels are reduced, ranging from 0.9% to 8.1%; (2) compared with the least accurate individual metamodels, the means of RMSE of the least accurate ensemble metamodels are reduced, ranging from 23.4% to 53.8%; (3) except for BP3, more than two ensemble metamodels perform better than the most accurate individual metamodels; (4) all the ensemble metamodels perform better than the four individual metamodels; (5) PrsRbfm is the best performer among the four metamodels, while SvrRbfi is the worst performer.

Table 8. Normalized means of RMSE of the metamodels for each benchmark problem with $\frac{5(k+1)(k+2)}{4}$ training points.

BPs	BP1	BP2	BP3	BP4	BP5	BP6
PRS	1.268	1.577	1.115	1.022	1.000	1.000
SVR	1.248	1.000	1.242	1.206	1.027	1.174
RBFM	1.000	1.020	1.000	1.000	1.134	1.379
RBFI	1.100	1.126	1.517	2.151	2.030	2.093
PrsRbfm	0.937	1.008	1.015	0.991	0.918	0.991
PrsRbfi	0.980	1.098	1.079	1.014	0.989	0.995
SvrRbfm	0.965	0.924	1.053	1.035	0.831	1.073
SvrRbfi	1.002	0.938	1.167	1.152	0.951	1.139
Best & Best	−6.3%	−7.6%	1.5%	−0.9%	−16.9%	−0.9%
Worst & Worst	−20.9%	−30.4%	−23.1%	−46.4%	−51.3%	−45.6%

shows the normalized means of RMSE of the metamodels for each benchmark problem with $\frac{5(k+1)(k+2)}{4}$ training points. It can be seen that: (1) compared with the most accurate individual metamodels, the means of RMSE of the most accurate ensemble metamodels are reduced for five benchmark problems, ranging from 0.9% to 16.9%; (2) compared with the least accurate individual metamodels, the means of RMSE of the least accurate ensemble metamodels are reduced, ranging from 20.9% to 51.3%; (3) all the ensemble metamodels have better performance than the four individual metamodels.

In summary, the choice of different sampling densities influences the results slightly, but the conclusions obtained by the three schemes with different sampling densities remain unchanged.

4.4. Significance of Results

The results above have proven the effectiveness of the proposed method to some extent. To further demonstrate the advantages, the proposed method is compared with some other popular ensemble metamodels, which are BPS (Best PRESS surrogate), PWS (PRESS weighted average surrogate), and OWSD (Optimal weighted surrogate using the diagonal elements). The detailed descriptions of these ensemble metamodels can be found in relevant literature [35,37]. Additionally, Kriging with first order polynomial regression function (KRG1) and Kriging with second-order polynomial regression function (KRG2) are also included in the performance comparison. To be noted, the principle and modeling process of Kriging are different from that of the proposed metamodeling approach in this paper.

Figure 6 compares the performance of PrsRbfm, SvrRbfm, KRG1, KRG2, BPS, PWS, and OWSD. It can be seen that: (1) for BP1, PrsRbfm and SvrRbfm perform better than the other five metamodels; (2) for BP2, SvrRbfm and BPS are the best two performers; (3) for BP3, the accuracy of PrsRbfm and BPS are better than that of the other metamodels; (4) for BP4, PrsRbfm and KRG2 are the best two performers; (5) for BP5, SvrRbfm and BPS are more accurate than other metamodels; (6) for BP6, PrsRbfm and KRG2 perform better the other metamodels.

In summary, the proposed metamodeling approach possesses some advantages when compared with KRG1, KRG2, BPS, PWS, and OWSD.

5. Conclusions

This paper proposed a novel metamodeling approach for building ensemble metamodels. Four types of ensemble metamodels, namely PrsRbfm, PrsRbfi, SvrRbfm, and SvrRbfi, were constructed by choosing four individual metamodels, namely PRS, SVR, RBFM, and RBFI. The performance of these metamodels was investigated through six popular benchmark problems. The effects of the performance criteria and sampling densities on the performance of the metamodels were studied. Additionally, the significance of the results was discussed by comparing the proposed method with some other popular ensemble metamodels. According to the results, some findings of this work could be concluded as follows:

(1)

According to the experimental results, the proposed metamodeling approach could reduce the risk of choosing the worst individual metamodel, and the constructed ensemble metamodels perform better than the selected individual metamodels in terms of accuracy.

(2)

The reason for the better performance of the ensemble metamodels may be that the proposed metamodeling approach combines the advantages of the regression-type and interpolation-type metamodels. The ensemble metamodels not only capture the global trend of the actual model in the entire design space, but also achieve the local accuracy in the vicinity of sampling locations.

(3)

The choices of different performance criteria and sampling densities influence the results slightly, but the obtained conclusions remain unchanged.

(4)

The proposed metamodeling approach possesses some advantages when compared with some other popular ensemble metamodels.