An Improved Blind Kriging Surrogate Model for Design Optimization Problems

Abstract

Surrogate modeling techniques are widely employed in solving constrained expensive black-box optimization problems. Therein, Kriging is among the most popular surrogates in which the trend function is considered as a constant mean. However, it also encounters several challenges related to capturing the overall trend with a relatively limited number of function evaluations as well as searching feasible points with complex or discontinuous feasible regions. To address this above issue, this paper presents an improved surrogate blind Kriging (IBK) and a combined infill strategy to find the optimal solution. According to enhancing the prediction accuracy of metamodels of objective and constraints, the high-order effects of regression function in the blind Kriging are identified by promising a variable selection technique. In addition, an infill strategy is developed based on the probability of feasibility, penalization, and constrained expected improvement for updating blind Kriging metamodels of the objective and constraints. At each iteration, two infill sample points are allocated at the positions to achieve improvement in optimality and feasibility. The IBK metamodels are updated by the newly-added infill sample points, which leads the proposed framework search to rapidly converge to the optimal solution. The performance and applicability of the proposed model are tested on several numerical benchmark problems via comparing with other metamodel-based constrained optimization methods. The obtained results indicate that IBK generally has a greater efficiency performance and outperforms the competitors in terms of a limited number of function evaluations. Finally, IBK is successfully applied to structural design optimization. The optimization results show that IBK is able to find the best feasible design with fewer evaluation functions compared with other studies, and this demonstrates the effectiveness and practicality of the proposed model for solving the constrained expensive black-box engineering design optimization problems.

1. Introduction

The optimization process is critical in the engineering design, which requires lower computational cost, robustness, stability, and accuracy. However, most conventional optimization techniques face several challenges in solving the black-box global optimization problem, such as the unavailable expressions, gradient information of objective function or constraints, and the time-consuming effort. One way to circumvent these issues is based on the surrogate model (SM) [1,2]. Here, SM is an invaluable tool that is used to approximate the expensive computational models during the optimization process indirectly or directly. There are commonly used metamodels, such as response surface methodology [3,4], radial basis function (RBF) [5,6], Kriging [7,8], support vector regression (SVR) [9], neural network (NN) [10,11,12,13], inverse distance weighting (IDW) [14], and so on. To achieve the optimal solution, the sampling technique and infill strategy are employed to build and iteratively refine the metamodel in order to ameliorate the solution found during the performance process. Therein, the sampling technique, which is also known as the design of the experiment, creates a set of points over the domain, and the initial metamodel is fitted to the observed points. Various sampling techniques for generating the sample points are developed, including Monte Carlo, random, and Latin hypercubic sampling (LHS) techniques. For more details, interested readers are suggested to refer to ref. [15].

In the past decade, variants of the metamodel-based optimization method have been successfully developed for solving the optimization problems with expensive simulations. For instance, Jones et al. [16] proposed an efficient global optimization (EGO) method where Kriging and expected improvement (EI) are introduced to find the solution. In addition, Huang et al. [17,18] extended the EGO algorithm by using the augmented EI function to determine the next sampling point. Additionally, a new infill strategy with adaptive radius and direction search is developed by Dong et al. [19] to indicate all local optimal values. To reduce the dimension of hyper-parameters, Zhao [20] incorporated the maximal information coefficient into Kriging. Gutmann [21] associated the measure of bumpiness in the RBF model to find the global minimum. Regis and Shoemaker [22,23,24] proposed a stochastic response surface to identify the promising points. More recently, Shepard [14] introduced the IDW, and then Joseph [25] added a linear regression function to improve accuracy. A conjugate between IDW and RBF is suggested by Bemporad [26].

The constrained black-box optimization problem is popular in many practical engineering designs. Hence, it has attracted the remarkable attention of researchers in recent years. They have been successfully applied to address this issue, such as ConstrLMSRBF [27], COBRA [28], RCGO [6], CARS [29], and so on. Among the surrogate models, Kriging has attracted more attention than other models due to estimating the prediction error and capacity to approximate the highly nonlinear functions. Therein, its regression function can be considered as the part trying to catch the general trend and thus the largest variations of the data. In recent times, Li et al. [30] proposed a new Kriging-based constrained global optimization algorithm in which the global optimal solution is obtained by two pivotal phases. An enhanced approach based on the modification of the probability of an improvement algorithm is presented by Carpio et al. Carpio et al. [31]. Additionally, Qian et al. [32] illustrated a new infill strategy, in which the position of the new sampling point is the intersection between the confidence interval and the constraint boundary. A combination of Kriging and the mixture of experts is proposed by Bartoli et al. [33] to improve the model’s accuracy. In addition, Forrester and Keane [34] presented a constrained expected improvement criterion to attain a new sample for updating the surrogate model. Similarly, Shi et al. [35] developed a probability of constrained improvement based on filter technology. However, there are several challenges due to multimodal, nonlinear functions, limited sample size, and a strong overall trend exists. In addition, identifying the right regression function for a set of data with interactions between variables is a difficult task. To address the above issues, researchers have tried to improve the prediction accuracy of Kriging by adjusting either the regression function, the stochastic process, or both. A well-known author, Joseph et al. [36], presented the blind Kriging model, in which the optimal basis functions are estimated by the Bayesian variable selection technique. Additionally, Kersaudy et al. [37] introduced a combination of polynomial chaos expansions and universal Kriging. Zhang et al. [38,39] proposed the regularization method for constructing a trend function and the penalized blind likelihood Kriging. Nevertheless, the above-mentioned models are only established to solve the unconstrained optimization problem. Furthermore, to the best of our knowledge, it has still not been yet utilized for constrained optimization thus far.

In this study, an improved surrogate blind Kriging is first presented to handle computational expensive constrained optimization problems. In our work, the higher-order effects of the trend function are estimated by a Bayesian variable selection technique. Simultaneously, a new infill strategy is developed based on the probability of feasibility, penalization, and constrained expected improvement to refine the surrogate models and improve the solution found during the optimization process. According to this scheme, two sample points will be indicated corresponding to the exploration and exploitation stages for each iteration. It helps the IBK metamodel to learn the complex or discontinuous feasible domains. The performance and applicability of the proposed approach are demonstrated through several benchmark optimization problems. The outcomes of the proposed model are compared with other models to evaluate its efficiency and reliability. The obtained results showed that our work is able to find the best feasible design with fewer evaluation functions.

The rest of the paper is organized as follows. Section 2 introduces the IBK model. Next, an infill strategy based on the combination is presented in Section 3. Afterward, several numerical examples are investigated to demonstrate the efficiency of the proposed model in Section 4. Finally, the conclusions are outlined in Section 5.

2. Improved Blind Kriging

2.1. Basics

The primary objective of a metamodel is to search for the approximate function from a set of data points. A brief summary of Kriging and some basic formulas are described in this subsection.

Consider the function $\mathbf{y}\left(\mathbf{x}\right)$ as defined in the domain $\Omega$. Let the n sample points $\mathbf{X}={\left\{{\mathbf{x}}^{\left(1\right)},{\mathbf{x}}^{\left(2\right)},...,{\mathbf{x}}^{\left(n\right)}\right\}}^T$ with $\mathbf{x}\in {\mathbb{R}}^d$ and the true response values $\mathit{y}={\left\{y^{\left(1\right)},y^{\left(2\right)},...,y^{\left(n\right)}\right\}}^T$ with $\mathbf{y}\in \mathbb{R}$, respectively. The universal Kriging (UK) model postulates a combination of the regression function and a stochastic process, as demonstrated in Equation (1)

$$Y\left(\mathbf{x}\right)=f\left(\mathbf{x}\right)+Z\left(\mathbf{x}\right),$$

(1)

where $Y\left(\mathbf{x}\right)$ represents a black-box function; $f\left(\mathbf{x}\right)=\nu {\left(\mathbf{x}\right)}^T{\mathit{\beta }}_m$ denotes the regression function; $Z\left(\mathbf{x}\right)$ is characterized by Gaussian process with zero mean and stationary covariance. For the blind Kriging, the trend function utilizes a regression function, as defined as Equation (2)

$$f\left(\mathbf{x}\right)=\nu {\left(\mathbf{x}\right)}^T{\mathit{\beta }}_m,$$

(2)

in which $\nu \left(\mathbf{x}\right)={\left\{1,{\nu }_1\left(\mathbf{x}\right),...,{\nu }_m\left(\mathbf{x}\right)\right\}}^T$ is the vector that contains known basic functions; m is the number of basis functions, and ${\mathit{\beta }}_m={\left\{{\mathit{\beta }}_0,{\mathit{\beta }}_1,...,{\mathit{\beta }}_m\right\}}^T$ denotes the unknown coefficients [40]. All of them are identified through the feature selection methods in Section 2.2.

Let $\Psi$ be the correlation matrix of the samples, $\mathit{\psi }\left(\mathbf{x}\right)$ the correlation vector between a point $\mathbf{x}$ and the sample point in the data set, and ${\mathit{F}}_m$ the model matrix of the sample points. Then the predicted value at any point x by the blind Kriging model is given as,

$$\widehat{y}\left(\mathbf{x}\right)=\nu {\left(\mathbf{x}\right)}^T{\widehat{\mathit{\beta }}}_m^{}+\mathit{\psi }{\left(\mathbf{x}\right)}^T{\Psi }^{-1}\left(\mathbf{y}-{\mathit{F}}_m^{}{\widehat{\mathit{\beta }}}_m^{}\right),$$

(3)

where the variance ${\widehat{\sigma }}_m^2$ and mean value ${\widehat{\mathit{\beta }}}_m^{}$ of variance can be calculated as follows

$$\begin{array}{c}{\widehat{\mathit{\beta }}}_m^{}=\frac{{\mathit{F}}_m^T{\Psi }^{-1}\mathbf{y}}{{\mathit{F}}_m^T{\Psi }^{-1}{\mathit{F}}_m}\hfill \\ {\widehat{\sigma }}_m^2=\frac{{\left(\mathbf{y}-{\mathit{F}}_m^{}{\widehat{\mathit{\beta }}}_m^{}\right)}^T{\Psi }^{-1}\left(\mathbf{y}-{\mathit{F}}_m^{}{\widehat{\mathit{\beta }}}_m^{}\right)}{n}\hfill \end{array}.$$

(4)

To indicate the optimal hyper-parameters $\theta$, in this paper, the differential evolution (DE) algorithm is employed to maximize the likelihood function. Their searching space is restricted to $\left({10}^{-3},{10}^2\right)$ [34,41]. It should be noted that $\theta$ is estimated two times. The first time, they are found at m = 0, and after choosing m, the hyper-parameters will be estimated again. The interested reader is referred to [36,39,40] for more details about deriving them.

Mean square error (MSE) at points in the design domain is given by [42]

$$\begin{array}{c}{\widehat{s}}^2\left(\mathbf{x}\right)={\widehat{\sigma }}_m^2\left[1-\mathit{\psi }{\left(\mathbf{x}\right)}^T{\Psi }^{-1}\mathit{\psi }\left(\mathbf{x}\right)+\left({\mathit{F}}_m^T{\Psi }^{-1}\mathit{\psi }\left(\mathbf{x}\right)\right.\right.\hfill \\ {\left.-\nu \left(\mathbf{x}\right)\right)}^T{\left({\mathit{F}}_m^T{\Psi }^{-1}{\mathit{F}}_m^{}\right)}^{-1}\left({\mathit{F}}_m^T{\Psi }^{-1}\mathit{\psi }\left(\mathbf{x}\right)-\nu \left(\mathbf{x}\right)\right).\hfill \end{array}$$

(5)

To provide insight, a single variable test function from five initial samples, which is assumed as a black-box objective function, is investigated here, as shown in Figure 1. In which the solid black line represents the true function, the blue dashed line shows the prediction value, the red dots denote the set of sample points, and the solid red line shows the root mean squared error. It is easily seen that the predicted value at sample points is equal to the real value, and its MSE is zero. Note that MSE is used as the measure of uncertainty in each prediction and to estimate the approximate accuracy of the BK model.

2.2. Variable Selection

As mentioned by Couckuyt [40], BK can capture the most variance in the data set by data analysis methods. A collection of candidate functions is considered for selection by using cross-validation prediction error (CVPE). Now, we consider a trend function that combines the regression function with the set of candidate functions to fit the data in the following linear model

$$f\left(\mathbf{x}\right)=\sum _{i=0}^m{\mathit{\beta }}_i{\nu }_i\left(\mathbf{x}\right)+\sum _{i=0}^t{\alpha }_i{u}_i\left(\mathbf{x}\right),$$

(6)

where t is the number of potential functions; $u\left(\mathbf{x}\right)=\left\{1,u_1\left(\mathbf{x}\right),...,u_t\left(\mathbf{x}\right)\right\}$ denotes the set of candidate functions, and $\alpha$ is the vector of corresponding coefficients. Note that $\mathit{\beta }$ has already been determined independently $\alpha$. When the transfer function f can be highly nonlinear, the number of candidate features to approximate are larger than the number of sample points. Therefore, we cannot determine all the coefficients ${\alpha }_i$. In this study, a function prior to all parameters in the linear model is utilized to overcome this issue. Additionally, the trend function is improved by considering the high-order effect of the Bayesian variable selection technique, including linear, quadratic, cubic, quartic, and two-factor interaction effects. Consequently, the total of candidate variables and the mean term is $t=4\left(d^2+2{\displaystyle \sum _{i=1}^{d-1}}i\right)$. Let us consider five-level, spaced factors with possible levels ${\mathit{l}}^{\left(1\right)},{\mathit{l}}^{\left(2\right)},{\mathit{l}}^{\left(3\right)},{\mathit{l}}^{\left(4\right)}$, and ${\mathit{l}}^{\left(5\right)}$, which are defined as follows

$$\begin{array}{c}{\mathit{l}}^{\left(1\right)}=min\left(\mathbf{X}\right),\\ {\mathit{l}}^{\left(2\right)}=\frac{mean\left(\mathbf{X}\right)}{2},\\ {\mathit{l}}^{\left(3\right)}=mean\left(\mathbf{X}\right),\\ {\mathit{l}}^{\left(4\right)}=\frac{3mean\left(\mathbf{X}\right)}{2},\\ {\mathit{l}}^{\left(5\right)}=max\left(\mathbf{X}\right).\end{array}$$

(7)

Therefore, the model matrix using orthogonal polynomial coding, which is the column lengths of $\sqrt{5}$, is

$${\mathit{U}}_j=\left[\begin{array}{ccccc}1& -\sqrt{2}& \sqrt{\frac{10}{7}}& -\sqrt{\frac{1}{2}}& \sqrt{\frac{1}{14}}\\ 1& -\sqrt{\frac{1}{2}}& -\sqrt{\frac{5}{14}}& \sqrt{2}& -\sqrt{\frac{8}{7}}\\ 1& 0& -\sqrt{\frac{10}{7}}& 0& \sqrt{\frac{18}{7}}\\ 1& \sqrt{\frac{1}{2}}& -\sqrt{\frac{5}{14}}& -\sqrt{2}& -\sqrt{\frac{8}{7}}\\ 1& \sqrt{2}& \sqrt{\frac{10}{7}}& \sqrt{\frac{1}{2}}& \sqrt{\frac{1}{14}}\end{array}\right].$$

(8)

The corresponding correlation matrix is

$${\Psi }_j=\left[\begin{array}{ccccc}1& {\psi }_j\left(l_j^{\left(2\right)}\right)& {\psi }_j\left(l_j^{\left(3\right)}\right)& {\psi }_j\left(l_j^{\left(4\right)}\right)& {\psi }_j\left(l_j^{\left(5\right)}\right)\\ {\psi }_j\left(l_j^{\left(2\right)}\right)& 1& {\psi }_j\left(l_j^{\left(2\right)}\right)& {\psi }_j\left(l_j^{\left(3\right)}\right)& {\psi }_j\left(l_j^{\left(4\right)}\right)\\ {\psi }_j\left(l_j^{\left(3\right)}\right)& {\psi }_j\left(l_j^{\left(2\right)}\right)& 1& {\psi }_j\left(l_j^{\left(2\right)}\right)& {\psi }_j\left(l_j^{\left(3\right)}\right)\\ {\psi }_j\left(l_j^{\left(4\right)}\right)& {\psi }_j\left(l_j^{\left(3\right)}\right)& {\psi }_j\left(l_j^{\left(2\right)}\right)& 1& {\psi }_j\left(l_j^{\left(2\right)}\right)\\ {\psi }_j\left(l_j^{\left(5\right)}\right)& {\psi }_j\left(l_j^{\left(4\right)}\right)& {\psi }_j\left(l_j^{\left(3\right)}\right)& {\psi }_j\left(l_j^{\left(2\right)}\right)& 1\end{array}\right],$$

(9)

and the ${\mathit{R}}_j\left(5\times 5\right)$ variance-covariance matrix is calculated as follows

$${\mathit{R}}_j={\mathit{U}}_j^{-1}{\Psi }_j{\left({\mathit{U}}_j^{-1}\right)}^T.$$

(10)

Note that we only use the diagonal elements of ${\mathit{R}}_j$ [43]. Hence, the ith element of the diagonal matrix $\mathit{R}\left(\left(t+1\right)\times \left(t+1\right)\right)$ can be written as

$${\mathit{R}}_{i,i}=\prod _{j=1}^d{r}_{jl}^{l_{ij}}{r}_{jq}^{q_{ij}}{r}_{jc}^{c_{ij}}{r}_{j{q}_r}^{{q_r}_{ij}},$$

(11)

in which $l_{ij}=1$ if ${\alpha }_i$ includes the factor’s linear effect j and 0 otherwise. Similarly, $q_{ij},c_{ij},{q_r}_{ij}$ are assigned to the quadratic, cubic, and quartic effects, respectively.

$$\begin{array}{c}{r}_{jl}^{}=\frac{{\mathit{R}}_j\left(2,2\right)}{{\mathit{R}}_j\left(1,1\right)},\hfill \\ r_{jq}^{}=\frac{{\mathit{R}}_j\left(3,3\right)}{{\mathit{R}}_j\left(1,1\right)},\hfill \\ r_{jc}^{}=\frac{{\mathit{R}}_j\left(4,4\right)}{{\mathit{R}}_j\left(1,1\right)},\hfill \\ r_{j{q}_r}^{}=\frac{{\mathit{R}}_j\left(5,5\right)}{{\mathit{R}}_j\left(1,1\right)},\hfill \\ 0⩽l_{ij}+q_{ij}+c_{ij}+{q_r}_{ij}⩽2.\end{array}$$

(12)

The sample data will be normalized to the interval [0, 1]. The encoded samples for the linear, quadratic, cubic, and quartic effects are expressed as follows

$$\begin{array}{c}{x}_{j,l}=\frac{1}{\sqrt{2}}\left(\frac{x_j}{l_j^{\left(5\right)}}\right),\hfill \\ x_{j,q}=\sqrt{\frac{5}{14}}\left[{\left(\frac{x_j}{l_j^{\left(5\right)}}\right)}^2-2\right],\hfill \\ x_{j,c}=\frac{5}{6\sqrt{2}}\left[{\left(\frac{x_j}{l_j^{\left(5\right)}}\right)}^3-\frac{17}{5}\left(\frac{x_j}{l_j^{\left(5\right)}}\right)\right],\hfill \\ x_{j,q_r}=\frac{35}{12\sqrt{14}}\left[{\left(\frac{x_j}{l_j^{\left(5\right)}}\right)}^4-\frac{31}{7}{\left(\frac{x_j}{l_j^{\left(5\right)}}\right)}^2+\frac{72}{35}\right],\hfill \end{array}$$

(13)

where $x_j$ defines the jth column X. The two-factor interaction terms can be built from these basic effects. After the completion of the construction of R as shown above, the posterior mean of $\mathbf{\alpha }$ is defined as

$$\widehat{\mathbf{\alpha }}={\left(\frac{{\tau }_m}{{\sigma }_m}\right)}^2\mathit{R}{{\mathit{F}}_c}^{\prime }{\Psi }^{-1}\left(\mathbf{y}-{\mathit{F}}_m{\mathit{\beta }}_m\right),$$

(14)

$$var\left(\widehat{\mathbf{\alpha }}\right)={\tau }_m^2\left(\mathit{R}-{\left(\frac{{\tau }_m}{{\sigma }_m}\right)}^2\mathit{R}{{\mathit{F}}_c}^{\prime }{\Psi }^{-1}{\mathit{F}}_c\mathit{R}\right),$$

(15)

where ${\mathit{F}}_c$ is the matrix of all candidate variables. As shown in refs. [36,40,43], the absolute of $\widehat{\alpha }$ is used instead of the standardized coefficient for variable selection. Hence, variables selected at each step correspond to the largest value $\widehat{\alpha }$. It should be noted that $\frac{{\tau }_m}{{\sigma }_m}$ is a constant, which is set to 1 for simplicity of the computation. To be specific, for the best value of m, the blind Kriging model has to go through five stages, which are described in detail as follows.

Step 1: The LHS technique is used to create initial sample points X and the obtained response values, respectively.

Step 2: Construct the ordinary Kriging surrogate model, and the leave-one-out cross-validation prediction error $\left(CVP{E}_{(m=0)}\right)$ is estimated to measure accuracy.

Step 3: Determine the coefficients $\left|\widehat{\alpha }\right|$ corresponding to each promising feature, and arrange $\left|{\widehat{\alpha }}_i\right|$ from largest to smallest.

Step 4: If improving the accuracy of prediction has not been satisfied,

Step 4.1: Add a candidate function corresponding to each $\left|{\widehat{\alpha }}_i\right|$ coefficient to the regression function.

Step 4.2: Construct the intermediate Kriging surrogate model with the new regression function.

Step 4.3: $CVP{E}_{\left(m\right)}$ is estimated to measure the accuracy of the model,

$$\begin{array}{c}CVP{E}_{\left(m\right)}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n}c{v}_i^2}{n}},\hfill \\ c{v}_i^{}=y_i\left({\mathbf{x}}_i\right)-{\widehat{y}}_i\left({\mathbf{x}}_i\right).\end{array}$$

(16)

Step 5: The best set of features to minimize $\left(CVP{E}_{\left(m\right)}\right)$, will be chosen to build the final BK model. Finally, the new hyper-parameters are indicated for the new trend function [40] of the surrogate model.

3. Infill Strategy

The infill strategy is employed to iteratively refine the metamodel and guide the algorithm toward a promising region to improve the solution found. Its trade-off between exploration and exploitation decides the success of the implementation. Although there are several available infill strategies, such as expected improvement (EI), expected violation, probability of improvement function (PI), lower confidence bounding, etc., these procedures may be difficult to find feasible points or feasible disconnected regions [6,28]. In this section, a new combined infill strategy, which consists of the probability of feasibility, penalization, and a constrained expected improvement, is introduced to overcome the above limitations.

The mathematical formulation of the constrained expensive black-box global optimization problem

$$\begin{array}{c}\underset{\mathbf{x}\in {\mathbb{R}}^d}{min}f\left(\mathbf{x}\right),\hfill \\ s.t.g_j\left(\mathbf{x}\right)\le 0,j=1,2,...,m,\hfill \\ {\mathbf{x}}_{lb}\le \mathbf{x}\le {\mathbf{x}}_{ub},\hfill \end{array}$$

(17)

where f and $g_j$ denote the objective function and inequality constraints, respectively. Normally, if there is one equality constraint, it is converted to two inequality constraints using a very small relaxation factor. ${\mathbf{x}}_{ub}$ and ${\mathbf{x}}_{lb}$ are the upper and lower bounds of the design variables, respectively. It should be noted that the above problem belongs to the class of derivative-free optimization, so its derivatives are not available. In order to solve Equation (17), the constrained EI is usually employed as the infill criterion and expressed as follows

$$E{I}_C\left(\mathbf{x}\right)=E{I}_{C_i}\left(\mathbf{x}\right)+E{I}_{C_r}\left(\mathbf{x}\right),$$

(18)

with

$$\begin{array}{cc}\hfill & \begin{array}{c}E{I}_{C_i}\left(\mathbf{x}\right)=\left(f\left({\mathbf{x}}^{*}\right)-{\mu }_n\left(\mathbf{x}\right)\right)\Phi \left(\frac{f\left({\mathbf{x}}^{*}\right)-{\mu }_n\left(\mathbf{x}\right)}{{\sigma }_n^{}\left(\mathbf{x}\right)}\right)PF\left(\mathbf{x}\right),\hfill \\ E{I}_{C_r}\left(\mathbf{x}\right)={\sigma }_n^{}\left(\mathbf{x}\right)\varphi \left(\frac{f\left({\mathbf{x}}^{*}\right)-{\mu }_n\left(\mathbf{x}\right)}{{\sigma }_n^{}\left(x\right)}\right)PF\left(\mathbf{x}\right),\hfill \\ PF\left(\mathbf{x}\right)=\prod _{j=1}^m{Pr}_j\left(g_j\left(\mathbf{x}\right)\le 0\right),\hfill \\ {Pr}_j\left(g_j\left(\mathbf{x}\right)\le 0\right)=\Phi \left(\frac{-{\mu }_{g_j}\left(\mathbf{x}\right)}{{\sigma }_{g_j}^{}\left(\mathbf{x}\right)}\right),\hfill \end{array}\hfill \end{array}$$

where ${\mu }_{g_j}\left(\mathbf{x}\right)$ and ${\sigma }_{g_j}^{}\left(\mathbf{x}\right)$ are the posterior mean and covariance function of the jth constraint, respectively.

As indicated by Haftka et al. [44], $E{I}_{C_i}\left(\mathbf{x}\right)$ indicates the exploitation sample points where the posterior mean of the objective function is small and has low uncertainty, and the constraints are likely to be near the boundaries of feasible operation. In contrast, $E{I}_{C_r}\left(\mathbf{x}\right)$ aims to identify exploration sample points, which tend to be biased towards regions that are more likely to satisfy the constraints and high uncertainty of the objective function. However, the standard EI${}_C$ often fails when there is no initial feasible point in the data set [44]. To circumvent this, $PF\left(\mathbf{x}\right)$ is one of the most used methods to favor the feasible region. Note that the value of $PF\left(\mathbf{x}\right)$ is easy to be zero when one of the constraints violates the design specification [44]. Consequently, the high-priority data points are missed or overlooked.

To tackle these limitations, a combined computation strategy is first introduced in this study to achieve an initial feasible point as well as find the optimal solution. This scheme helps to reduce the number of constraint violations and explore other feasible regions, as follows:

$$PFC\left(\mathbf{x}\right)={\sigma }_n^{}\left(\mathbf{x}\right)PF\left(\mathbf{x}\right),$$

(19)

$$PIC\left(\mathbf{x}\right)=-\sum _{j=1}^{m}max\left(0,{\mu }_{g_j}\left(\mathbf{x}\right)\right),$$

(20)

$$IS{C}_{14}\left(\mathbf{x}\right)=\prod _{j=1}^m{Pr}_j\left(g_j\left(\mathbf{x}\right)\le 0\right).D_j,$$

(21)

where $PFC\left(\mathbf{x}\right)$ is combined by the covariance function and the probability of feasibility function; $PIC\left(\mathbf{x}\right)$ is the penalization function corresponding to the constraint set; $IS{C}_{14}\left(\mathbf{x}\right)$ is the infill sampling criterion for disconnected feasible regions; D is the distance to the nearest feasible point, and given by:

$$D=\underset{{\mathbf{x}}_{feas}}{min}\left(\frac{∥{\mathbf{x}}_{feas}-\mathbf{x}∥}{range}\right),$$

(22)

in which range is the lag distance at which it reaches the sill from the variogram model.

According to the proposed strategy, if the feasible point is not found in the data, $PFC\left(\mathbf{x}\right)$ and $PIC\left(\mathbf{x}\right)$ will be used as an alternative for obtaining two infill sample points. In this case, providing that there exists a feasible region, two feasible sample points will be identified. More specifically, the first point is located in the sparsely sampled area by maximizing $PFC\left(\mathbf{x}\right)$, so ${\sigma }_n^{}\left(\mathbf{x}\right)$ characterizes the sample density of the objective function in the design space. The other point indicates the corresponding minimum value of the total constraints $PIC\left(\mathbf{x}\right)$. When the models fit poorly and do not have a feasible region, $PIC\left(\mathbf{x}\right)$ is utilized to improve the locations that violated constraints. On the contrary, once an initial feasible point is found, Equations (18) and (21) are employed to determine infill sample points. $IS{C}_{14}\left(\mathbf{x}\right)$ aims to effectively explore other feasible regions. The algorithm framework is provided in Algorithm 1.

Algorithm 1 The combined infill strategy for the expensive constrained optimization problem

Input:  d        : number of design variables       m        : number of constraints       $N_{max}$: maximum expensive evaluation number       $\left[{\mathbf{x}}_{lb};{\mathbf{x}}_{ub}\right]$: design space Output:  $[{\mathbf{x}}_{best},{\mathbf{f}}_{best}]$: the optimum solution 1: Generate $n_0$ sample points using LHS from the design space, evaluate the fitness values, and collect a set of initial observations ${\mathcal{D}}_{n_0}^{}=\left\{{\mathbf{x}}_i,y_i,{\mathbf{g}}_i\right\}\forall i=1,...,n_0$ 2: Determine ${\mathcal{D}}_{}^{fes}$ from ${\mathcal{D}}_{n_0}^{}$ 3: Set $n=n_0$ 4: while$n\le N_{max}$do 5:    Build surrogate models for the objective and each of the constraint functions 6:    if ${\mathcal{D}}_n^{fes}=\varnothing$ then 7:      Find ${\mathbf{x}}_{n+1}$ by maximizing Equation (19) 8:      Find ${\mathbf{x}}_{n+2}$ by maximizing Equation (20) 9:    else 10:      Find ${\mathbf{x}}_{n+1}$ by maximizing Equation (18) 11:      Find ${\mathbf{x}}_{n+2}$ by maximizing Equation (21) 12:    end if 13:    Evaluate function values at ${\mathbf{x}}_{n+1,2}$, append ${\mathcal{D}}_{n+2}^{}={\mathcal{D}}_n^{}\cup \left\{\left({\mathbf{x}}_{n+1,2},y_{n+1,2},{\mathbf{G}}_{n+1,2}\right)\right\}$, and update ${\mathcal{D}}_{n+1,2}^{fes}$ 14:    $n=n+2$ 15:    Update $[{\mathbf{x}}_{best},{\mathbf{f}}_{best}]$ 16: end while

4. Numerical Examples

4.1. Study in High-Order Effects of the Trend Function

To evaluate the high-order efficiency of the trend function, a planar truss structure for maximum passive vibration isolation [45], as shown in Figure 2, is introduced. The authors of Keane and Bright [46] achieved results for this problem through analysis and experiment. The structures consist of 42 members with the same material properties. The right-hand end nodes are fixed. The unit force excitation is subjected to node 11 over a frequency range 100–200Hz. The x- and y-coordinates of nodes seventh) and eighteenth are considered as design variables, and the other nodes are fixed as per the regular structure. The objective function is stress at the left-hand end node. The interested reader is referred to [45,46] for more details. LHS is used to obtain data for building metamodels that differ in the number of sample points, and 100 validation runs were performed to validate the results. In this case, the average Euclidean error (AEE) is used to measure the prediction error using

$$AEE(\widehat{y},y)=\frac{1}{k}\sum _{i=1}^k\sqrt{{\left(\widehat{y}-y\right)}^2}.$$

(23)

As shown in Figure 3, it is easily seen that the BK model improves compared to the Kriging model. Clearly, the CVPE is smaller than the Kriging model. It is easy to explain that the linear, quadratic effect causes a significant rise in accuracy. On the other hand, the difference in AEE scores between the models was an insignificant amount with the number of sample points less than 30. This is explained by the fact that if the sample size is small, then the high-order effect of the regression function will not be much. However, the high-order BK model achieved a vast improvement over another model when the data size increased. Then, the regression function with the high-order effects is represented as Equations (24) and (25) for 100 and 200 sample points, respectively. Clearly, IBK can capture the overall trend to the largest variations of the data.

$$\begin{array}{c}1+{\mathrm{x}}_1{\mathrm{x}}_3+{\mathrm{x}}_1+{\mathrm{x}}_4+{\mathrm{x}}_1{\mathrm{x}}_4+{\mathrm{x}}_3{\mathrm{x}}_4+{\mathrm{x}}_2+{\mathrm{x}}_1{\mathrm{x}}_2^3{+\mathrm{x}}_1^4{\mathrm{x}}_2^4+{\mathrm{x}}_2{\mathrm{x}}_4{+\mathrm{x}}_1^3{\mathrm{x}}_2{+\mathrm{x}}_1^4{\mathrm{x}}_3{+\mathrm{x}}_1^2{\mathrm{x}}_3{+\mathrm{x}}_1^2{\mathrm{x}}_4.\hfill \end{array}$$

(24)

$$\begin{array}{c}1+{\mathrm{x}}_1{\mathrm{x}}_4+{\mathrm{x}}_1+{\mathrm{x}}_1{\mathrm{x}}_3{+\mathrm{x}}_1^4{\mathrm{x}}_3+{\mathrm{x}}_1{\mathrm{x}}_2+{\mathrm{x}}_3{\mathrm{x}}_4{+\mathrm{x}}_1^4{\mathrm{x}}_2^4+{\mathrm{x}}_2{+\mathrm{x}}_1^4{\mathrm{x}}_3^2{+\mathrm{x}}_2^4{\mathrm{x}}_4+{\mathrm{x}}_1{\mathrm{x}}_3^2{+\mathrm{x}}_2^4{\mathrm{x}}_4^2.\hfill \end{array}$$

(25)

4.2. Synthetic Test Problems

In this section, nine well-known benchmark problems, in which the number of variables, constraints, and the best know solutions are given

Table 1. Constrained optimization benchmark problems.

Problems	Best Known Value	Dimensionality	No. Constraints
G24	−5.508	2	2
G8	−0.0958	2	2
Two_bars	0.0309	3	2
Four_bars	1400	4	1
Ibeam	0.0131	4	2
G5MOD	5126.5	4	5
G4	−30,665.539	5	6
Hesse	−310	6	6
SR	2994.42	7	11

, are investigated to evaluate the efficiency of the proposed model. More information about these problems can be seen in Appendix A and in ref. [28]. Although these optimization problems are not expensive to evaluate, they are treated as all computationally expensive functions to conduct meaningful comparisons of the performance of the alternative methods. The obtained results will be compared with the result of RCGO, FLT-AKM, COBRA-local, COBRA-global, and ConstrLMSRBF. Therein, the two algorithms RCGO and FLT-AKM are implemented by the authors following Wu et al. [6] and Shi et al. [35]. In the RCGO method, the distance coefficient in the first phase is $\Pi =[0.001,0.005,0.01,0.05,0.1]$, and $\Gamma =[0,1,2,3,4]$ is the value of the exponent in the second phase. In order to obtain a fair comparison between the different methods, all the tests are performed in a Matlab^TM environment using an Intel^® Core(TM) i5-8500 CPU 3.0 GHz desktop machine (Microsoft, New York, NY, USA). Each metamodel is run ten times on each test problem to reduce the effect of random error. Otherwise, the searching process of the model is terminated when either the maximum number of evaluation functions is reached or the best result does not improve after 10 iterations. LHS created ten different initial data sets. The differential evolution algorithm (DE) is employed to find hyper-parameters and the next points in the infill strategy.

The obtained results are summarized in

Table 2. Comparison of the obtained results for the test problems.

		Test Problems
		G24	G8	Four Bars	Two Bars	Ibeam	G5MOD	G4	Hesse	SR
IBK	Best	−5.4669	−0.095	1400	0.031	0.013	5126.5	−30,665.5	−310	3042.2
	Median	−5.4075	−0.091	1400	0.031	0.014	5135.2	−30,665.5	−310	3051.3
	Mean	−5.3832	−0.089	1400	0.031	0.014	5167.1	−30,665.5	−310	3053.3
	NEF	33.2	46.8	17.4	36.2	60	74.0	24.2	38.6	68.4
FLT-AKM	Best	−5.46	−0.095	1400	0.031	0.013	5739.9	−30,611.7	−306.553	3024.4
	Median	−5.4173	−0.084	1400	0.031	0.015	6157.9	−30,416.8	−297.34	3061
	Mean	−5.415	−0.085	1400	0.031	0.018	6187.8	−30,411.2	−297.093	3058.5
	NEF	34	54.8	17.5	37.8	74.8	59.2	35.2	58.6	59.6
RCGO	Best		−0.096	1400	0.031	0.013	5323.5	−30,628.8	−306.428
	Median		−0.096	1400	0.031	0.014	6099.8	−30,389.5	−294.812
	Mean	N/A(10)	−0.073	1400	0.033	0.016	5992.9	−30,395.3	−296.239	N/A(10)
	NEF		51.7	17.2	40.4	67.8	65.4(3)	35	73.2
COBRA-Local	Best	-	−0.1	-	-	-	5126.5	−30,665.5	−309.94	2994.4
	Median		−0.1				5126.51	−30,665.2	−297.87	2994.7
[28]	Mean		−0.09				5126.51	−30,665.1	−296.25	2994.7
	NEF		50				50	50	50	50
COBRA-Global	Best	-	−0.1	-	-	-	5126.5	−30,665.4	−309.97	2994.7
	Median		−0.1				5126.51	−30,664.9	−297.87	2994.7
[28]	Mean		−0.09				5126.62	−30,664.9	−296.25	2994.7
	NEF		50				50	50	50	50

, which provides the average number of evaluation functions (NEF), objective values of the best solution, median solutions, and mean solutions found by the six metamodels. Firstly, it is easily seen that RCGO has failed to indicate the feasible solution for G24, SR in all trials, and G5MOD in three trials. These are mainly due to the small feasible region, all the initial points are infeasible, or the incoherent distance coefficients. Clearly, IBK outperformed other algorithms for G24, G8, Two-bars, G4, Ibeam, and Hesse with the best solution. In addition, the IBK requires fewer evaluations than other metamodels. Figure 4 provides a graphical illustration of the overall process of identifying the infill samples of the G24 problem. Firstly, the feasible domain is discontinuous and includes one global and two local minimums corresponding to peaks of the feasible domain. All initial sample points are quite analogous to the metamodels for a fair comparison. Clearly, the infill samples obtained by IBK are distributed near the peaks of the feasible region as well as the boundary of the constraint set. On the contrary, the obtained result from FLT-AKM divided two groups of the infill samples, and the first group concentrates on the infeasible region far from the constrained boundary. This can easily be explained by the fact that the order of constraints is polynomials of degree four. Hence, IBK can easily capture the high-order effects compared with Kriging of the FLT-AKM. Furthermore, our infill strategy allows us to reduce the number of constraint violations and explore other feasible domains. More specifically, our model requires only 32 evaluation functions, while FLT-AKM requires 45 analyses for the convergence performance. This shows that the features of IBK and the combination infill strategy outperform other well-known existing algorithms.

For the Four-bars test, the obtained results of all these models are similar, and there are not many differences between the models. It is interesting that, among the problems, it requires few evaluation functions. Hence, IBK will not show much effect with small sample points or a simple, feasible domain. It can be observed that the optimal results found by the COBRA provide the best solution to the G5MOD and SR problems. However, IBK performs better than FLT-AKM and RCGO in terms of median, mean and best values of feasible optimum. Note that COBRA and RCGO require the initial data that contain the feasible points. However, it is difficult to find a feasible point for highly nonlinear problems, which then leads to failure. Clearly, RCGO could not find a feasible point in all trials for the SR problem. Although the result achieved from SR by IBK is close to FLT-AKM, IBK shows stability when the median and mean are close to the best value.

4.3. Structural Design Optimization

Finally, the proposed model is applied to perform the shape optimization of truss structures. This benchmark problem has been previously analyzed by Shao [48] and Suprayitno [49] using the clustering-based surrogate model and Kriging, respectively. A 21-bars planar truss shown in Figure 5 is investigated for design optimization. Coordinates of all free nodes are treated as continuous design variables. The design aim is to maximize the truss performance, which is defined as the maximum load ratio ($F_{max}$) that the truss can carry to the self-weight of the truss (W). The structure is subjected to stress and stability ratio constraints. The stress limitations of all members are the same for tension and compression and set to 200 MPa. The cross-sectional areas of all bars are 25 cm${}^2$. The number of initial sample points is set at 14 [48]. Here, the critical buckling stress $S_{cr,i}={\pi }^{2}E/{\left(L_{e,i}/\rho \right)}^2$. Whereas E and $\rho$ are Young’s modulus of the truss material and the radius of gyration of the truss bar, respectively. It can be formulated as follows

$$\begin{array}{c}\\ Maximize\frac{F_{max}}{W}=f\left(x_1,x_2,x_3,y_3,x_4,y_4,y_6\right),\hfill \\ subjectedto\left\{{\sigma }_t,{\sigma }_c\right\}⩽{\sigma }_a,\hfill \\ \frac{{\sigma }_{c,i}}{S_{cr,i}}⩽1,\hfill \\ 0.8⩽x_1,x_3,y_3,y_4,y_6⩽3m,\hfill \\ 3.3⩽x_2,x_4⩽5.5m.\hfill \end{array}$$

(26)

The number of maximum evaluation functions is set to 100 for all models. The obtained optimal results are listed in

Table 3. Comparison of the obtained results for the 21-bars planar truss.

Design Variables	IBK	FLT-AKM	RCGO	CMLS [48]	EORKS [49]
$x_1$	3	3	2.98	2.37	3
$x_2$	4.31	4.251	4.19	4.174	4.209
$x_3$	1.907	1.91	1.91	1.877	1.902
$y_3$	2.267	2.275	2.278	1.721	2.285
$x_4$	4.218	4.218	4.237	4.024	4.222
$y_4$	1.706	1.708	1.705	1.438	1.709
$y_6$	2.723	2.702	2.624	2.186	2.644
Best	25.85	25.46	24.72	22.79	25.66
Median	25.27	25.21	23.46	-	-
Mean	25.36	25.33	23.75	-	-
NFE	55.6	65.7	72.3	95	87

and Figure 6. It can easily be seen that IBK found the best feasible solution, where the objective value is 25.8541, while the maximum tensile stress and maximum compressive stress are 94.21 and -115.701 Mpa, respectively. In addition, IBK significantly decreased the number of simulations. These have demonstrated that IBK can solve expensive black-box structure design optimization problems effectively.

5. Conclusions

In this study, an improved surrogate blind Kriging combined with the new infill strategy is introduced to solve the constrained expensive black-box optimization problems. According to that core idea, the high-order effects of the trend function are identified from the training data using the Bayesian variable selection method. Simultaneously, the novel infill strategy based on the probability of feasibility, penalization, and constrained expected improvement is constructed to reduce the number of constraint violations, explore other feasible domains, and handle the discontinuous design domain. The efficiency of the proposed approach is demonstrated through several numerical examples for synthetic test problems and shapes the optimization of the truss structure. The obtained results have indicated that the optimum solution obtained by this work is in good agreement with ten of the eleven tests. The proposed model saves a number of evaluation functions in almost all problems in comparison with other metamodels. In addition, it is also shown that the IBK predictor is simpler to interpret and is more robust than an ordinary Kriging predictor. This proved to be an efficient model for solving real-world engineering problems. Hence, it promises to be robust and effective at resolving complex problems in the expensive optimization of a complex design domain.