## 2 Background  
### Kriging Model  
The Kriging model was originally proposed by [PERSON] [38] to predict the location of a mine hole in a geostatistical community. Then, it was extended by [PERSON] et al. [6] for modeling an experiment of a computer. The Kriging model is also called the Gaussian process model, which is a kind of interpolative model. The Kriging model can be expressed as  
\[\hat{y}(\mathbf{x})=\beta+Z(\mathbf{x}) \tag{1}\]  
where \(\mathbf{x}\) represents the vector of the design variables, which is a \(d\)-dimensional vector \(\mathbf{x}=\{x_{1},x_{2},\ldots,x_{d}\}\), \(\beta\) is an unknown parameter which denotes the global tendency, \(Z()\) is a static Gaussian process with zero mean and non-zero variance \(\sigma^{2}\), which represents the local deviation. In the static Gaussian process, spatial correlation is used to organize the relationship between any two samples. Generally, the squared exponential function is utilized, which can be expressed as  
\[R(\mathbf{x}_{i},\mathbf{x}_{j};\mathbf{\theta})=\exp(-\sum\limits_{k=1}^{d}\theta_{k}(x_ {i}^{k}-x_{j}^{k})^{p_{j}}) \tag{2}\]  
where \(\mathbf{\theta}\) and \(\mathbf{P}\) are the hyper-parameters used to control the smooth and the correlation between two sample points. Generally, the hyper-parameter vector \(\mathbf{P}\) is set to be \(p_{i}=2;i=1,2,\ldots,d\)[39]. The core point of the modeling process of the Kriging model is to determine the unknown parameters. Because the responses obey the multivariable Gaussian distribution, the unknown parameter can be obtained by maximum likelihood estimation (MLE) [14]. The likelihood function can be organized as  
\[L\big{(}y(\mathbf{x}_{1}),y(\mathbf{x}_{2}),\ldots,(x_{N})\big{|}\sigma,\beta,\mathbf{ \theta}\big{)}=\frac{1}{(2\pi\sigma^{2})^{\frac{N}{2}}}\exp[-\frac{\sum\limits _{i=1}^{N}\big{(}\mathbf{y}(\mathbf{x}_{i})-\beta\big{)}^{2}}{2\sigma^{2}}] \tag{3}\]  
where \(N\) is the number of samples. Then, Equation (3) can be simplified by taking the natural logarithm,  
\[\ln(L)=-\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln\big{(}\sigma^{2}\big{)}-\frac{1}{ 2}\ln(|\mathbf{R}|)-\frac{(\mathbf{y}-\mathbf{1}\beta)^{T}\mathbf{R}^{-1}(\mathbf{y}-\mathbf{1}\beta)} {2\sigma^{2}} \tag{4}\]  
where \(\mathbf{y}\) is an \(N\)-dimensional vector that consists of the real responses, \(1\) is an \(N\)-dimensional vector that consists of \(1\),  
The values of \(\beta\) and \(\sigma^{2}\) can be obtained by setting the derivatives of Equation (4) concerning \(\beta\) and \(\sigma^{2}\) to be \(0\),  
\[\hat{\beta}=\frac{\mathbf{f}^{T}\mathbf{R}^{-1}\mathbf{y}}{\mathbf{f}^{T}\mathbf{R}^{-1}\mathbf{f}} \tag{5}\]\[\partial^{2}=\frac{\left(\mathbf{y}-\mathbf{1}\hat{\beta}\right)^{T}\!\mathbf{R}^{-1}\!\left( \mathbf{y}-\mathbf{1}\hat{\beta}\right)}{N} \tag{6}\]  
Then, substituting Equations (5) and (6) into Equation (4), and remove the constant terms, Equation (4) yields the concentrated \(\ln\)-likelihood function  
\[\ln(L)=-\frac{N}{2}\ln\!\left(\partial^{2}\right)-\frac{1}{2}\ln\!\left(|R|\right) \tag{7}\]  
It is difficult to obtain an analytical solution of \(\Theta\) because of high non-linearity and non-differentiality. Therefore, a numerical solution is obtained instead. The optimization algorithm, such as the genetic algorithm (GA) [40] and particle swarm optimization algorithm (PSO) [41], can be used to find the optimized values of \(\Theta\). The Kriging model is widely adopted in surrogate model-based engineering optimization because it can provide both the predicted value and variance [42]. The predicted value of an un- sampled point can be determined by minimizing the mean square error.