Ensemble of surrogates with hybrid method using global and local measures for engineering design

Abstract

Surrogate models are usually used as a time-saving approach to reduce the computational burden of expensive computer simulations for engineering design. However, it is difficult to choose an appropriate model for an unknown design space. To tackle this problem, an effective method is forming an ensemble model that combines several surrogate models. Many efforts were made to determine the weight factors of ensemble, which include global and local measures. This article investigates the characteristics of global and local measures, and presents a new ensemble model which combines the advantages of these two measures. In the proposed method, the design space is divided into two parts, and different strategies are introduced to evaluate the weight factors in these two parts respectively. The results from numerical and engineering design cases show that the proposed ensemble model has satisfactory robustness and accuracy (it performs best for most cases tested in this article), while spending almost the equivalent modeling time (the additional cost is not more than 6.7% for any case tested in this article) compared with the combined global and local ensemble models.

Appendices

Appendix A: Description of the selected surrogate models

In this appendix, a brief overview of the mathematical formulations of PRS, RBF and KRG surrogate models is provided.

1.1 A.1 Polynomial response surface (PRS)

The PRS approximation is one of the most well-established surrogate models. The most commonly used PRS model is the following second-order form

$$ \widehat{f}(x)={\beta}_0+\sum \limits_{i=1}^d{\beta}_i{x}_i+\sum \limits_{i=1}^d\sum \limits_{j\ge i}^d{\beta}_{ij}{x}_i{x}_j $$

(27)

where d is the number of design variables, β ₀, β _i and β _ij are the unknown coefficients to be determined by the least squares technique. Here we run the MATLAB® routine “pinv” to obtain the Moore-Penrose generalized inverse matrix of the unknown coefficients, which was proved to be the optimal least square solution by Penrose (1955).

1.2 A.2 Radial basis function (RBF)

RBF models were originally developed to approximate multivariate functions. The general form of the RBF approximation can be expressed as

$$ \widehat{f}(x)=\sum \limits_{i=1}^n{w}_i\varphi \left(\left\Vert x-{x}_i\right\Vert \right) $$

(28)

where n denotes the number of sample points, w _i are the unknown coefficients to be determined, ‖·‖ represents the Euclidean norm and φ(·) is the so-called basis function. Powell (1987) suggested several forms of the basis function φ(·):

• Gaussian $\varphi (r)={e}^{-\frac{r^2}{2{\sigma}^2}}$

• ultiquadric $\varphi (r)=\sqrt{r^2+{\sigma}^2}$

• nverse Multiquadric $\varphi (r)=1/\sqrt{r^2+{\sigma}^2}$

• Thin-Plate Spline φ(r) = r ² ln(r)

where σ ≥ 0. In this article we use the multiquadric basis function with σ = 1 (suggested by Acar and Rais-Rohani (2009)), for its prediction accuracy and convergence ability with increased sample points. In order to obtain the unknown coefficients w _i , we substitute the n sample points into the (28) to form an equation as

$$ y=\varPhi \cdot w $$

(29)

where y is the vector of sample responses and Φ is an n × n matrix of basis functions. The coefficients vector w is obtained by solving (29).

1.3 A.3 Kriging (KRG)

The basic assumption of KRG is the estimation of the response in the form

$$ Y(x)=\mu (x)+Z(x) $$

(30)

where the response Y consists of a known polynomial μ(x) which globally approximates the trend of the function and a stochastic component Z(x) which generates deviations, so that the Kriging model interpolates the sample points. The correlation between the random variables Y(x⁽ⁱ⁾) and Y(x^(j)) is given by

$$ Corr\left[Y\left({\mathrm{x}}^{(i)}\right),Y\left({\mathrm{x}}^{(j)}\right)\right]=\exp \left(-\sum \limits_l^d{\theta}_l{\left|{x}_l^{(i)}-{x}_l^{(j)}\right|}^{p_l}\right) $$

(31)

where d is the number of design variables, θ _l and P _l (l = 1, ⋯, d) are unknown parameters to be estimated. Here we only consider a constant term to represent the mean of the overall surface (the Ordinary Kriging) and fix the parameters P _l  = 2(l = 1, ⋯, d) (the stationary Gaussian correlation function case). Then we search the optimal θ _l in the range of [10⁻³, 10²] (suggested by Forrester et al. (2008)) with the GA (Genetic Algorithm) toolbox of MATLAB®.

Once the correlation function has been selected, the response is predicted as

$$ \widehat{y}=\widehat{\mu}+{r}^T{R}^{-1}\left(y-\mathbf{1}\widehat{\mu}\right) $$

(32)

where the matrix R ⁻¹ is the inverse of the correlation matrix R whose element R _ij is equal to the (31), y is the vector of sample responses and 1 represents an n × 1 vector of ones. The estimated value of $\widehat{\mu}$ and the expressions of r ^T are

$$ \widehat{\mu}=\frac{{\mathbf{1}}^T{R}^{-1}y}{{\mathbf{1}}^T{R}^{-1}\mathbf{1}} $$

(33)

$$ {r}^T=\left( Corr\left[Y\left({x}^{(1)}\right),Y(x)\right]\kern0.5em \cdots \kern0.5em Corr\left[Y\left({x}^{(n)}\right),Y(x)\right]\right) $$

(34)

Detailed derivation of Kriging can be found in Jones (2001) and Forrester et al. (2008).

Appendix B: Description of the numerical test functions

In this appendix, the description of six numerical test functions is provided. The landscapes of two-variable functions are depicted in Figs. 6 and 7.

Fig. 6

The landscape of Branin-Hoo function

Fig. 7

The landscape of Camelback function

Fig. 8

Boxplots of RMSE for six numerical examples

Fig. 9

Boxplots of AAE for six numerical examples

Fig. 10

Boxplots of MAE for six numerical examples

1.1 B.1 Branin-Hoo function

$$ f(x)={\left({x}_2-\frac{5.1{x}_1^2}{4{\pi}^2}+\frac{5{x}_1}{\pi }-6\right)}^2+10\left(1-\frac{1}{8\pi}\right)\cos \left({x}_1\right)+10 $$

(35)

where x ₁ ∈ [−5, 10] and x ₂ ∈ [0, 15]_.

1.2 B.2 Camelback function

$$ f(x)=\left(4-2.1{x}_1^2+\frac{x_1^4}{3}\right){x}_1^2+{x}_1{x}_2+\left(-4+4{x}_2^2\right){x}_2^2 $$

(36)

where x ₁ ∈ [−2, 2] and x ₂ ∈ [−2, 2]_.

1.3 B.3 and B.4 Hartman functions

$$ f(x)=-\sum \limits_{i=1}^4{c}_i\exp \left[-\sum \limits_{j=1}^n{a}_{ij}{\left({x}_j-{p}_{ij}\right)}^2\right] $$

(37)

where x _i  ∈ [0, 1]. Two types of Hartman functions are given based on different number of input variables: (1) Hartman-3 with three input variables (test function 3), and (2) Hartman-6 with six input variables (test function 4). While the parameter c in each function is the same vector ${\left[1\kern0.5em 1.2\kern0.5em 3\kern0.5em 3.2\right]}^T$, the other two parameters a and p are shown in Table 5 and 6.

Table 5 Parameters used in Hartman-3 function, j = 1, 2, 3

Table 6 Parameters used in Hartman-6 function, j = 1, 2, ⋯, 6

1.4 B.5 Extended-Rosenbrock function

$$ f(x)=\sum \limits_{i=1}^{m-1}\left[{\left(1-{x}_i\right)}^2+100{\left({x}_{i+1}-{x}_i^2\right)}^2\right] $$

(38)

where x _i  ∈ [−5, 10], i = 1, 2, ⋯, m = 9.

1.5 B.6 Dixon-Price function

$$ f(x)={\left({x}_1-1\right)}^2+\sum \limits_{i=2}^mi{\left[2{x}_i^2-{x}_{i-1}\right]}^2 $$

(39)

$$ {x}_i\in \left[-10,10\right],i=1,2,\cdots, m=12. $$

(where)

Appendix C: Test results for determining the form of ES-HGL

In this appendix, the test results referenced in Section 3.3 are provided in Table 7, 8 and 9.

Table 7 Test results of different region radius forms

Table 8 Test results of different HWF forms

Table 9 Test result for evaluating the effect of the hybrid method

1.1 C.1 Test result for determining the form of the region radius

1.2 C.2 Test result for determining the form of the HybridWeight Factor (HWF)

1.3 C.3 Test result for evaluating the effect of the hybrid method

Appendix D: Boxplots for six numerical examples (Figures 8, 9, and 10)