Ensemble of surrogate based global optimization methods using hierarchical design space reduction

Abstract

For computationally expensive black-box problems, surrogate models are widely employed to reduce the needed computation time and efforts during the search of the global optimum. However, the construction of an effective surrogate model over a large design space remains a challenge in many cases. In this work, a new global optimization method using an ensemble of surrogates and hierarchical design space reduction is proposed to deal with the optimization problems with computation-intensive, black-box objective functions. During the search, an ensemble of three representative surrogate techniques with optimized weight factors is used for selecting promising sample points, narrowing down space exploration and identifying the global optimum. The design space is classified into: Original Global Space (OGS), Promising Joint Space (PJS), Important Local Space (ILS), using the newly proposed hierarchical design space reduction (HSR). Tested using eighteen representative benchmark and two engineering design optimization problems, the newly proposed global optimization method shows improved capability in identifying promising search area and reducing design space, and superior search efficiency and robustness in identifying the global optimum.

Appendices

Appendix 1: Description of selected surrogate modeling techniques

1.1 Polynomial response surface, PRS

Polynomial response surface, also called response surface method (RSM) is one of the most popular surrogate modeling techniques. It first generates a number of design points, then uses either a second-order model or a higher order model to replace the unknown system. The most commonly used PRS model is the second-order model in the form as:

$$ \tilde{y}\left(\boldsymbol{x}\right)={\beta}_0+\sum \limits_{i=1}^n{\beta}_i{\boldsymbol{x}}_i+\sum \limits_{i=1}^n{\beta}_{ii}{\boldsymbol{x}}_i^2+\sum \limits_{i<j}\sum \limits_{i,j=1}^n{\beta}_{ij}{\boldsymbol{x}}_i{\boldsymbol{x}}_j $$

(19)

where n is the number of variables in the input vector x, $\tilde{y}\left(\boldsymbol{x}\right)$ is the response surface approximation of the actual function. β₀, β_i, β_ii, β_ij represent the unknown regression coefficients which can be computed by the least squares method.

1.2 Radial basis function, RBF

Radial basis function is originally developed by Hardy in 1971 to fit irregular topographic contours of geographical data. It has been known tested and verified for several decades and many positive properties have been identified.

For a data set consisting of the values of design variables and response values at N sample points, the actual function can be approximated as

$$ \tilde{y}\left(\boldsymbol{x}\right)=\sum \limits_{i=1}^N{\lambda}_i\phi \left(\left\Vert \boldsymbol{x}\hbox{-} {\boldsymbol{x}}_i\right\Vert \right) $$

(20)

where x is the vector of design variables, x_i is the vector of design variables at the ith point, $\left\Vert \boldsymbol{x}\hbox{-} {\boldsymbol{x}}_i\right\Vert =\sqrt{\left(\boldsymbol{x}\hbox{-} {\boldsymbol{x}}_i\right){\left(\boldsymbol{x}\hbox{-} {\boldsymbol{x}}_i\right)}^T}$ is the Euclidean norm, ϕ is a basis function, and λ_i is the coefficient for the ith basis function. The approximation function $\tilde{y}\left(\boldsymbol{x}\right)$ is a linear combination of some basis functions with weight coefficients λ_i. The most commonly used basis functions are listed in Table 10. The Multiquadric is selected as basis function due to its prediction accuracy and simple structure in this paper. Many literatures find that c₁ = 1 is suitable for most function approximations (Mullur and Messac 2006; McDonald et al. 2007; Kim et al. 2009; Pan et al. 2014b). Thus, c₁ is set to 1 in this work.

Table 10 Commonly used basis functions

1.3 Kriging, KRG

Kriging estimates the value of a function as a combination of a known function like a linear model such as a polynomial trend and departures which represent the low and high frequency variation components, respectively. The form of KRG is usually expressed by

$$ \tilde{f}\left(\mathrm{x}\right)=\sum \limits_{i=\boldsymbol{1}}^m{\beta}_i{y}_i\left(\mathrm{x}\right)+Z\left(\mathrm{x}\right) $$

(21)

where β_i is an unknown constant estimated, Z (x) is assumed to be a realization of a stochastic process with zero mean and a nonzero covariance. The i, jth element of covariance matrix of Z(x) is formulized as

$$ Cov\left[Z\left({\mathrm{x}}^i\right),Z\left({\mathrm{x}}^j\right)\right]={\sigma}_z^{\boldsymbol{2}}{R}_{ij} $$

(22)

where ${\sigma}_z^{\boldsymbol{2}}$ is the process variance, and R_ij is the correlation function between the ith and jth data points. The Gaussian function is used as the correlation function in this study, given by

$$ R\left({\mathrm{x}}^i,{\mathrm{x}}^j\right)={R}_{ij}=\exp \left\{-\sum \limits_{k=1}^n{\theta}_k{\left({x}_k^i-{x}_k^j\right)}^{\boldsymbol{2}}\right\} $$

(23)

where θ_k is distinct for each dimension, and these unknown parameters are generally determined by solving a nonlinear optimization problem.

Appendix 2: List of benchmark optimization problems

• (1) Six-hump Camel-Back function (SC) with n = 2

$$ f(x)=4{x}_1^2-2.1{x}_1^4+{x}_1^6/3+{x}_1{x}_2-4{x}_2^2+4{x}_2^4\kern0.2em $$

(24)

• (2) Branin function (BR) with n = 2

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

(25)

• (3) Generalized polynomial function (GF) with n = 2

$$ f(x)={\left(1.5-{x}_1\left(1-{x}_2\right)\right)}^2+{\left(2.25-{x}_1\left(1-{x}_2^2\right)\right)}^2+{\left(2.625-{x}_1\left(1-{x}_2^3\right)\right)}^2 $$

(26)

• (4) Goldstein and Price function (GP) with n = 2

$$ {\displaystyle \begin{array}{l}f(x)=\left[1+{\left({x}_1+{x}_2+1\right)}^2\left(19-14{x}_1+3{x}_1^2-14{x}_2+6{x}_1{x}_2+3{x}_2^2\right)\right]\ast \\ {}\left[30+{\left(2{x}_1-3{x}_2\right)}^2\left(18-32{x}_1+12{x}_1^2+48{x}_2-36{x}_1{x}_2+27{x}_2^2\right)\right]\end{array}} $$

(27)

• (5) Rosenbrock function (RB) with n = 2

$$ f(x)=100{\left({x}_2-{x}_1^2\right)}^2+{\left({x}_1-1\right)}^2 $$

(28)

• (6) Himmelblau function (HM) with n = 2

$$ f(x)={\left({x}_1^2+{x}_2-11\right)}^2+{\left({x}_1+{x}_2^2-7\right)}^2 $$

(29)

• (7) Cross-IN-TRAY Function (CT) with n = 2

$$ f(x)=-0.0001{\left(\left|\sin \left({x}_1\right)\sin \left({x}_2\right)\exp \left(\left|100-\frac{\sqrt{x_1^2+{x}_2^2}}{\pi}\right|\right)\right|+1\right)}^{0.1} $$

(30)

• (8) Drop-Wave Function (DW) with n = 2

$$ f(x)=-\frac{1+\cos \left(12\sqrt{x_1^2+{x}_2^2}\right)}{0.5\left({x}_1^2+{x}_2^2\right)+2} $$

(31)

• (9) Levy Function (LF3) with n = 3

$$ {\displaystyle \begin{array}{c}f(x)={\sin}^2\left(\pi {w}_1\right)+\sum \limits_{i=1}^{n-1}{\left({w}_i-1\right)}^2\left[1+10{\sin}^2\left(\pi {w}_i+1\right)\right]+{\left({w}_n-1\right)}^2\left[1+{\sin}^2\left(2\pi {w}_n\right)\right],\\ {}\kern0.5em {w}_i=1+\frac{x_i-1}{4}\end{array}} $$

(32)

• (10) Ackley Function (AF4) with n = 4

$$ f(x)=-20\exp \left(-0.2\sqrt{\frac{1}{n}\sum \limits_{i=1}^n{x}_i^2}\right)-\exp \left(\frac{1}{n}\sum \limits_{i=1}^n\cos \left(2\pi {x}_i\right)\right)+20+\exp (1) $$

(33)

• (11) Hartman function (HN6) with n = 6

$$ {\displaystyle \begin{array}{l}f(x)=-\sum \limits_{i=1}^4{c}_i\exp \left[-\sum \limits_{j=1}^n{\alpha}_{ij}{\left({x}_j-{p}_{ij}\right)}^2\right],\kern0.4em i=1,2,\dots, n.\kern0.6em {c}_i=\left[1\kern0.5em 1.2\kern0.5em 3\kern0.5em 3.2\right]\\ {}\left[{a}_{ij}\right]=\left[\begin{array}{cccccc}10& 3& 17& 3.5& 1.7& 8\\ {}0.05& 10& 17& 0.1& 8& 14\\ {}3& 3.5& 1.7& 10& 17& 8\\ {}17& 8& 0.05& 10& 0.1& 14\end{array}\right],\left[{p}_{ij}\right]={\left[\begin{array}{cccccc}1312& 1696& 5569& 124& 8283& 5886\\ {}2329& 4135& 8307& 3736& 1004& 9991\\ {}2348& 1451& 3522& 2883& 3047& 6650\\ {}4047& 8828& 8732& 5743& 1091& 381\end{array}\right]}^{\ast }1{0}^{-4}\end{array}} $$

(34)

• (12) Trid function (TR6 and TR10) with n = 6, 10

$$ f(x)=\sum \limits_{i=1}^n{\left({x}_i-1\right)}^2-\sum \limits_{i=2}^n{x}_i{x}_{i-1} $$

(35)

• (13) Sum squares function (SF) with n = 12

$$ f(x)=\sum \limits_{i=1}^n{ix}_i^2\kern0.1em $$

(36)

• (14) Griewank function (GN4, GN8 and GN12) with n = 4, 8, 12

$$ f(x)=\frac{1}{4000}\sum \limits_{i=1}^n{x}_i^2-\prod \limits_{i=1}^n\cos \left(\frac{x_i}{\sqrt{i}}\right)+1 $$

(37)

• (15) A function of 16 variables (F16) with n = 16

$$ {\displaystyle \begin{array}{c}f(x)=\sum \limits_{i=1}^n\sum \limits_{j=1}^n{\alpha}_{ij}\left({x}_i^2+{x}_i+1\right)\left({x}_j^2+{x}_j+1\right),\kern0.4em i,j=1,\kern0.3em 2,\kern0.3em \dots, n\\ {}{a}_{ij\left(\mathrm{row}1-8\right)}=\left[\begin{array}{cccccccccccccccc}1& 0& 0& 1& 0& 0& 1& 1& 0& 0& 0& 0& 0& 0& 0& 1\\ {}0& 1& 1& 0& 0& 0& 1& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ {}0& 0& 1& 0& 0& 0& 1& 0& 1& 1& 0& 0& 0& 1& 0& 0\\ {}0& 0& 0& 1& 0& 0& 1& 0& 0& 0& 1& 0& 0& 0& 1& 0\\ {}0& 0& 0& 0& 1& 1& 0& 0& 0& 1& 0& 1& 0& 0& 0& 1\\ {}0& 0& 0& 0& 0& 1& 0& 1& 0& 0& 0& 0& 0& 0& 1& 0\\ {}0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 1& 0& 1& 0& 0& 0\\ {}0& 0& 0& 0& 0& 0& 0& 1& 0& 1& 0& 0& 0& 0& 1& 0\end{array}\right]\kern0.2em \\ {}{a}_{ij\left(\mathrm{row}\ 9-16\right)}=\left[\begin{array}{cccccccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 1& 0& 0& 0& 1\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 1& 0& 0& 0\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 1& 0& 0\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& 0& 0\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ {}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]\kern0.1em \end{array}} $$

(38)