## 2 Multi-Objective Reliability Optimization Method  
### Kriging Surrogate Model  
Kriging is one of the most used estimation methods for spatial data interpolation, and it has been widely used to build surrogate models in automotive engineering in recent years [21; 22; 23]. The Kriging surrogate model covers the global trend and local nonlinearity of the response, and especially has high accuracy for predicting nonlinear response. It can analyze the trend and dynamic characteristics of known information [24; 25; 26].  
The Kriging surrogate model usually includes a polynomial function and a random distribution. Its approximate function expression is:  
\[\hat{y}(\mathbf{x})=\mathbf{f}^{\mathrm{T}}(\mathbf{x})\cdot\mathbf{\beta}+\mathbf{z}(\mathbf{x}) \tag{1}\]  
where \(\mathbf{\beta}=\big{[}\beta_{1},\dots,\beta_{p}\big{]}^{\mathrm{T}}\) represents the regression coefficient vector; \(\mathbf{f}^{\mathrm{T}}(\mathbf{x})\) represents a polynomial with design vector \(\mathbf{x}\); \(\mathbf{z}(\mathbf{x})\) represents a random distribution, which can be expressed as a random function with zero mean and standard deviation of \(\sigma\).  
The covariance of random distribution \(\mathbf{z}(\mathbf{x})\) is:  
\[Cov\big{[}z(\mathbf{x}_{i}),z(\mathbf{x}_{j})\big{]}=\sigma^{2}\mathbf{R} \tag{2}\]  
where \(\mathbf{R}=\big{[}R(\mathbf{x}_{i},\mathbf{x}_{j})\big{]}\) is the spatial correlation equation of any two nondiagonal sample points \(\mathbf{x}_{i}\) and \(\mathbf{x}_{j}\) in the matrix, which plays an absolute role in the simulation accuracy.  
Using Gaussian correlation equation to express \(R\big{(}\mathbf{x}_{i},\mathbf{x}_{j}\big{)}\) is:  
\[R\big{(}\mathbf{x}_{i},\mathbf{x}_{j}\big{)}=EXP\Bigg{(}-\sum_{k=1}^{m}\lambda_{k} \Big{|}x_{ik}-x_{jk}\Big{|}^{2}\Bigg{)} \tag{3}\]  
where \(m\) is the number of design variables; \(\lambda_{k}\) represents the correlation coefficient of the fitted surrogate model; and \(x_{ik}\) and \(x_{jk}\) represent the \(k\)th value of \(x_{i}\) and \(x_{j}\), respectively.  
After determining the correlation function, the estimated value \(\hat{y}(\mathbf{x})\) of the approximate response can be obtained. The estimated value \(\hat{\mathbf{\beta}}\) of regression coefficient \(\mathbf{\beta}\) and the estimated value \(\hat{\sigma}^{2}\) can be expressed as:  
\[\left\{\begin{array}{l}\hat{y}(\mathbf{x})=\mathbf{f}^{\mathrm{T}}(\mathbf{x})\overset{ \mathbf{\wedge}}{\mathbf{\beta}}+\mathbf{r}^{\mathrm{T}}(\mathbf{X}^{*})R^{-1}\bigg{(}\mathbf{y}- \mathbf{F}\overset{\mathbf{\wedge}}{\mathbf{\beta}}\bigg{)}\\ \overset{\mathbf{\wedge}}{\mathbf{\beta}}=\big{(}\mathbf{F}^{\mathrm{T}}R^{-1}\mathbf{F}\big{)} ^{-1}\mathbf{F}^{\mathrm{T}}R^{-1}\mathbf{y}\\ \overset{\mathbf{\wedge}}{\mathbf{\beta}}=\frac{\big{(}\mathbf{y}-\mathbf{F}\overset{\mathbf{\wedge }}{\mathbf{\beta}}\big{)}^{\mathrm{T}}R^{-1}\big{(}\mathbf{y}-\mathbf{F}\overset{\mathbf{\wedge }}{\mathbf{\beta}}\big{)}}{\mathbf{\beta}}\end{array}\right. \tag{4}\]  
where \(\mathbf{r}^{\mathrm{T}}(\mathbf{X}^{*})=[R(\mathbf{X}^{*},\mathbf{x}_{1}),R(\mathbf{X}^{*},\mathbf{x}_{ 2}),\dots,R(\mathbf{X}^{*},\mathbf{x}_{n})]\) is the correlation coefficient vector between prediction point \(\mathbf{X}^{*}\) and sample point \(\mathbf{x}\); \(\mathbf{F}\) is the \((n\times p)\) design matrix.  
The maximum likelihood estimation for parameter \(\lambda_{k}\) is:  
\[\underset{\lambda_{k}\geq 0}{\max}(\lambda_{k})=-\frac{\big{[}n\ln\big{(}\hat{ \sigma}^{2}\big{)}+\ln|\mathbf{R}|\big{]}}{2} \tag{5}\]