A Semi-Parallel Active Learning Method Based on Kriging for Structural Reliability Analysis

Abstract

The reliability analysis system is currently evolving, and reliability analysis efforts are also focusing more on correctness and efficiency. The effectiveness of the active learning Kriging metamodel for the investigation of structural system reliability has been demonstrated. In order to effectively predict failure probability, a semi-parallel active learning method based on Kriging (SPAK) is developed in this study. The process creates a novel learning function called $U_{\mathrm{A}}$, which takes the correlation between training points and samples into account. The $U_{\mathrm{A}}$ function has been developed from the U function but is distinct from it. The $U_{\mathrm{A}}$ function improves the original U function, which pays too much attention to the area near the threshold and the accuracy of the surrogate model is improved. The semi-parallel learning method is then put forth, and it works since $U_{\mathrm{A}}$ and U functions are correlated. One or two training points will be added sparingly during the model learning iteration. It effectively lowers the required training points and iteration durations and increases the effectiveness of model building. Finally, three numerical examples and one engineering application are carried out to show the precision and effectiveness of the suggested method. In application, evaluation efficiency is increased by at least 14.5% and iteration efficiency increased by 35.7%. It can be found that the proposed algorithm is valuable for engineering applications.

1. Introduction

In engineering practice, reliability analysis is crucial in determining the likelihood that a system would fail due to random variables [1,2,3]. A performance function $G\left(x\right)$ characterizes the response of the system, when $G\left(x\right)>0$ the system is safe, otherwise, it is unsafe. The Monte Carlo Simulation (MCS) is still the most used reliability analysis technique [4,5,6]. It is thought to be the most precise and reliable method for system reliability analysis and comprises a random population $n_{MC}$ that is assessed based on the performance function. The failure probability $P_f$ is then calculated as the ratio of $n_{G\le 0}$ to the sampling $n_{MC}$ using Equation (1).

$$P_f=\frac{n_{G\le 0}}{n_{MC}}$$

(1)

In order to obtain a higher accuracy of analysis, the larger population $n_{MC}$ is better, i.e., $n_{MC}=1\times {10}^6$. However, it is computationally expensive and time-consuming to evaluate the performance function for the entire population, especially when performance function is implicit, is a complex system [7,8,9], or is the finite element analysis (FEA) [10,11]. It is highly necessary to study efficient reliability analysis approaches. Aside from numerical simulation methods such as MCS and its advanced tactics, such as importance sample (IS) [12], subset simulation (SS) [13], and line sampling [14], the most commonly utilized approaches are analytical methods and surrogate model methods.

The operational efficacy of analytical techniques, such as the first-order reliability method (FORM) and second-order reliability method (SORM), are considerable [11,15]. The Taylor approximation of the performance function at the most likely point is the basis of these techniques. However, for extremely nonlinear situations or when trying to solve the implicit function problem, the inaccuracy based on the methods could be unacceptable [16]. Due to these difficulties, the simple-to-evaluate surrogate model methods are an essential component of reliability analysis [11]. The surrogate model methods aim at approximating the performance function using the strategic design of experiments (DoE). The method results in a sufficiently accurate prediction of the performance function with a lower computing demand [17,18,19,20]. Over the last few decades, many surrogate model methods have been developed, including response surface methods (RSM) [21,22,23], polynomial chaos expansion (PCE) [24,25,26], support vector machine (SVM) [27,28,29,30], Kriging method [31,32,33,34], artificial neural networks (ANN) [35,36], and others.

Wang et al. [29] proposed a regression-model-adaptive algorithm based on the Bayesian support vector machine. The learning function of this algorithm uses penalty function to select sample points efficiently. Based on the generalized polynomial chaos expansion, Sepahvand et al. [37] performed stochastic structural modal analysis, which includes uncertain parameters. Taking advantage of the global prediction and local interpolation capability of PCE and Kriging, Aghatise et al. [26] proposed a hybrid method, which used the U learning function to drive the active learning function. Since there are multiple failure modes, Fauriat and Gayton [38] developed several Kriging metamodels and emphasized the AK-SYS sampling technique. However, in actual engineering, many studies combine IS and SS with Kriging surrogate model methods to solve the small failure probability problem [39]. Xiao et al. [40] combine a stratified IS method with the adaptive Kriging to estimate reliability.

Of the above methods, the Kriging method shows significant advantages [41,42]. First, Kriging is a reliable interpolation technique. The performance function’s precise value is predicted at a point in DoE [43,44,45,46]. Moreover, it estimates the local variance of the predictions in addition to the expected values without increasing the computational cost [11,41]. Thanks to the variance, the sequential construction of the surrogate model is carried out using the Kriging-based active learning reliability method. After the initial surrogate model is built with a small number of initial training points, more training points are added until the surrogate model accurately represents the performance function. In the process of developing a surrogate model, a learning function is essential for choosing the optimal training points. The criterion for choosing the following optimum training points is based on the learning function. Various Kriging-based surrogate model reliability techniques have been developed by various learning functions. Based on the cumulative distribution function (CDF) and probability density function (PDF), the expected improvement (EI) learning function is constructed [47,48], then the efficient global optimization (EGO) is proposed. On the basis of EGO, the efficient global reliability analysis (EGRA) uses the expected feasibility function to determine training points [18,49]. The active learning reliability method combines Kriging with Monte Carlo Simulation (AK-MCS) and constructs a learning function that can predict the probability of the correct sign of the performance function [4,50]. In addition, other leaning functions are also available, Lv [51] proposed an H learning function based on the information entropy. Sun [52] proposed an LIF learning function, and the studies on learning functions can be found in References.

The core of all the methods is the strategies of surrogate model construction and training points selection. EGO focuses on the expected improvement of model error, but the stopping criterion of learning process has no uniform standard, which may be difficult to converge for some problems. EGRA and LIF are similar to EGO. The theoretical accuracy of AK-MCS is 95.45%, and the stopping criterion is relatively stable. However, its learning process mainly focuses on the accuracy of the classification of test points, and lacks the accuracy of the overall model.

To further improve the efficiency of reliability methods, a semi-parallel active learning method based on Kriging is proposed in this paper. It takes into account the relationship between training points and sample points comprehensively but seldom increases the computational amount of the learning process and can consider more sample points to participate in the learning. Furthermore, the new training points are added based on the new learning process semi-parallel to update the surrogate model. Because of these innovations, it involves fewer iterations and is more efficient than traditional parallel learning or sequential learning.

The remainder of this paper is organized as follows: in Section 2, the Kriging method and current Kriging-based reliability techniques are discussed. The new method and the implemental procedure are described in Section 3. Three examples followed in Section 4 and an engineering application in Section 5 verifies the effectiveness and practicability of the proposed method. Conclusions are made in Section 6.

2. Review

2.1. Kriging

In place of the performance function, which was initially designed for forecasts in mining engineering and geostatistics, Kriging is beneficial for creating intricate and highly nonlinear models [53]. The approximated performance function is a realization of a statistical Gaussian process $G\left(x\right)$ as Equation (2):

$$G\left(x\right)=f{\left(x\right)}^T\beta +Z\left(x\right)$$

(2)

In the equation, $f\left(x\right)={\left[f_1\left(x\right),\dots ,f_P\left(x\right)\right]}^T$ there is a vector of the basic functions. $\beta$ is a vector of regression coefficients. $f{\left(x\right)}^T\beta$ forms an approximation of the response of the performance function in mean. $Z\left(x\right)$ is the stationary Gaussian process with zero mean and covariance between two points $x_i$ and $x_j$ calculated by Equation (3):

$$Cov\left[Z\left(x_i\right),Z\left(x_j\right)\right]={\sigma }_Z^2·R\left(x_i, x_j,\theta \right)$$

(3)

${\sigma }_Z^2$ is the variance of $Z\left(x\right)$. $R\left(x_i, x_j,\theta \right)$ is the stochastic process correlation function with respect to parameter $\theta$. The correlation function can be defined using a variety of models, including exponential, linear, spline, spherical, Gaussian, and so on. The Gaussian model is the most popular among them. The model can be formulated by Equation (4):

$$R\left(x_i, x_j,\theta \right)={\displaystyle \prod }_{m=1}^M\exp \left[-{\theta }_m{\left(x_{i,m}-x_{j,m}\right)}^2\right]$$

(4)

where $x_{i,m}$, $x_{j,m}$ and ${\theta }_m$ are the mth dimension of $x_i$, $x_j$ and $\theta$. $\left|x_{i,m}-x_{j,m}\right|$ represents the distance in the dimension $x_m$. ${\theta }_m$, which is a scale parameter that indicates the correlation between the points within dimension mth. The larger ${\theta }_m$ and $\left|x_{i,m}-x_{j,m}\right|$ are representative of a shorter correlation length. M is the number of input variables.

Based on a set of training points (i.e., [$x_1,x_2,\dots ,x_{\mathrm{n}}$]) with $x_i\in R^{\mathrm{n}}$ the ith experiment, and $Y$ with $Y_i=G\left(x_i\right)\in R$, $\beta$ and ${\sigma }^2$ are estimated according to Equation (5):

$$\widehat{\beta }={\left(F^{\mathrm{T}}{R}^{-1}F\right)}^{-1}{F}^{\mathrm{T}}{R}^{-1}Y$$

(5)

where $\widehat{\beta }$ is the generalized least square estimation of $\beta$, $F$ is a matrix with each row $f{\left(x\right)}^{\mathrm{T}}$, $R=R\left(x_i, x_j,\theta \right)$ is the matrix of correlation between each pair of points of the DoE.

Based on Equation (2) to Equation (5), the prediction of the Kriging model of the performance function with untried points $x$ can be described by Equation (6):

$$\widehat{G}\left(x\right)=f{\left(x\right)}^{\mathrm{T}}\widehat{\beta }+r{\left(x\right)}^{\mathrm{T}}{R}^{-1}\left(Y-F\widehat{\beta }\right)$$

(6)

$r\left(x\right)$ is the correlation vector between the untried points $x$ and each of the training points, which is estimated by Equation (7):

$$r\left(x\right)={\left[R\left(x, x_1,\theta \right),R\left(x, x_2,\theta \right),\dots ,R\left(x, x_{\mathrm{n}},\theta \right)\right]}^{\mathrm{T}}$$

(7)

The Kriging predicted variance is evaluated by Equation (8):

$${\widehat{\sigma }}_{\widehat{G}}^2\left(x\right)={\widehat{\sigma }}_Z^2[1+u^{\mathrm{T}}{\left(F^{\mathrm{T}}{R}^{-1}F\right)}^{-1}u-r{\left(x\right)}^{\mathrm{T}}{R}^{-1}r\left(x\right)]$$

(8)

where $u$ and ${\widehat{\sigma }}_Z^2$ can be estimated respectively by the following

$$u\left(x\right)=F^{\mathrm{T}}{R}^{-1}r\left(x\right)-f\left(x\right)$$

(9)

$${\widehat{\sigma }}_Z^2\left(x\right)=\frac{1}{\mathrm{n}}{\left(Y-F\widehat{\beta }\right)}^T{R}^{-1}\left(Y-F\widehat{\beta }\right)$$

(10)

Based on the above analysis, the Kriging prediction at an untried point $x$ is normally distributed [50] as Equation (11):

$$G\left(x\right)~N\left({\mu }_{\widehat{G}}\left(x\right),{\widehat{\mathsf{\sigma }}}_{\widehat{G}}^2\left(x\right)\right)$$

(11)

The prediction $\widehat{G}\left(x_i\right)$ in the training point $x_i$ of the design of experiments is exact, i.e., $\widehat{G}\left(x_i\right)=G\left(x_i\right)$. The variance is zero in these points. Therefore, it is an exact interpolation method, while in unexplored areas, the Kriging variance is predicted based on Equation (8). It enables the quantification of the uncertainty of local predictions. This property is important in reliability studies, which relate to the index to improve the surrogate model.

2.2. Kriging-Based Active Learning Algorithms

The accuracy of the Kriging surrogate model depends on the selection of training points. The same number of training points may build entirely distinct models with different training points. Additionally, because the initial training points are frequently chosen at random, they may not accurately reflect the properties of the performance functions. The Kriging-based active learning algorithm is researched to construct the surrogate model intelligently, accurately and efficiently. The primary active learning methods in this respect are the EI learning function [47,48], expected feasibility function (EFF), improved sequential Kriging reliability analysis (ISKRA), and U learning function [4]. These fundamental techniques are refined by numerous innovative active learning algorithms. The EI and U learning functions are briefly introduced in this section.

2.2.1. EI learning Function

The EI learning function is a method first presented by Jones, also known as the efficient global optimization (EGO) method [47]. It is similar to other learning functions that are used to decide where to add the new training point. It is described as the expected improvement that one point in the search space will provide a better solution than the current surrogate model. The expected values and variances predicted on the interpolating point are subjected to a normal distribution based on the aforementioned Equations (6) and (8). Assuming that the minimal performance function of the training points is $G_{min}$, the expected values and variances predicted on the interpolating point are $\widehat{G}\left(x\right)$ and $\widehat{\sigma }\left(x\right)$, respectively. The improvement of the objective function for the minimization issue is given by Equation (12):

$$I\left(x\right)=max(G_{min}-\widehat{G}\left(x\right),0)$$

(12)

Then, the expectation of $I\left(x\right)$ is described as Equation (13):

$$E\left[I\left(x\right)\right]=\{\begin{array}{cc}\left(G_{min}-\widehat{G}\left(x\right)\right)·\varnothing \left(\frac{G_{min}-\widehat{G}\left(x\right)}{\widehat{\sigma }\left(x\right)}\right)+\widehat{\sigma }·\phi \left(\frac{G_{min}-\widehat{G}\left(x\right)}{\widehat{\sigma }\left(x\right)}\right)& \widehat{\sigma }>0\\ 0& \widehat{\sigma }=0\end{array}$$

(13)

where $\varnothing$ and $\phi$ are the CDF and PDF of standard normal distribution, respectively. The new training point where $E\left[I\left(x\right)\right]$ is maximal, as shown by Equation (14), is the most effective new training point to increase the accuracy of the surrogate model.

$$x_{new}=arg max\left(E\left[I\left(x\right)\right]\right)$$

(14)

The Kriging surrogate model will be recreated with the new training point included. Additionally, the EI function will be rebuilt. The training will continue until the EI value at its maximum point falls below the set tolerance, which is typically $1\times {10}^{-6}$. Based on Equations (13) and (14), the new training point is chosen, where $\widehat{y}\left(x\right)$ is close to the minimum or $\widehat{y}\left(x\right)$ has a larger variance $\widehat{\sigma }\left(x\right)$ with a high possibility. It is efficient to obtain the global minimum of the performance function, but it is not efficient to find a specified threshold value of the reliability analysis. Then, EFF is acquired to provide an instruction of how well the true value of the response is expected to satisfy the equality constraint $\widehat{y}\left(x\right)=\overline{z}$. The details of the EFF can be obtained from refs [48,54].

2.2.2. U Learning Function

The sign of the prediction in MCS is crucial for the reliability analysis. The point with a high potential risk of crossing the threshold $G\left(x\right)=0$ should be added to the DoE. The potentially crossing points can present the properties of: (1) Being near $G\left(x\right)=0$ or the threshold value; (2) Being significantly uncertain or having high Kriging variance; (3) Both happen simultaneously. As a result, Equation (15) is proposed as the learning function for AK based on prediction $\widehat{G}\left(x\right)$ and variance $\widehat{\sigma }\left(x\right)$ [4].

$$U\left(x\right)=\frac{\left|\widehat{G}\left(x\right)\right|}{\widehat{\sigma }\left(x\right)}$$

(15)

Here, it is assumed that the prediction $\widehat{G}\left(x\right)$ has a mean of zero based on Equation (11), then $U\left(x\right)$ is obeyed to normal distribution with mean zero and covariance $\widehat{\sigma }\left(x\right)$. So, $\varnothing \left(U\left(x\right)\right)$ is the probability of the points $x$ correctly classified according to the sign of $\widehat{G}\left(x\right)$. When $U\left(x\right)$ is smaller, the probability of being correctly classified is lower. The new training point is discovered to be the one with the minimal $U\left(x\right)$. This process will continued until the stopping criterion $min\left(U\left(x\right)\right)\ge \gamma$ is satisfied. Commonly, $\gamma$ = 2 is suitable, with the probability of correct classification being equal to 95.45%. AK-MCS has shown effectiveness in a variety of challenging instances. However, for the global design space modeling and small failure probability problem, the population required for AK-MCS will be very large. To overcome this problem, some methods, for example combined with IS, can be implemented, referring to study [55] for details.

To approximate the performance function more accurately at the extreme region and critical region where $G\left(x\right)$ is close to zero, adding the EI learning point and U function learning point to be training points simultaneously has been considered. It requires two completely independent learning functions, as well as two iterating stopping criteria, which will lengthen the calculation time and perhaps result in non-convergence.

3. Proposed Method

3.1. New Learning Function

In this study, an advanced learning function is proposed based on U function. It is a novel method that can add new training points according to the requirements that further consider the relationship between training points and sample points. The learning function shown as Equation (16):

$$U_{\mathrm{A}}\left(x\right)=\frac{\left|\widehat{y}\left(x\right)\right|}{\widehat{\sigma }\left(x\right)}·r_m{\left(x\right)}^p$$

(16)

where, $r_m\left(x\right)$, maximal correlation coefficient, is the maximum of correlation vector $r\left(x\right)$ in Equation (7), $r_m\left(x\right)\in \left(0,1\right]$. $p$ the power exponent of $r_m\left(x\right)$ reflects the importance of the correlation coefficient, it can be any non-negative rational number, such as 0, 2, 4 and so on. In particular, when $p=0$, Equation (16) turns to the traditional U function. If $p$ is increased, the MCS sample points with small correlation coefficients that are far from the current training points and the critical region have a higher probability of being the new training points. $U_{\mathrm{A}}\left(x\right)$ pays more attention to the region where the points are uncertain. As a matter of experience, $p=2$ is usually appropriate. The stopping criterion of $U_{\mathrm{A}}\left(x\right)$ is identical to the U function shown as Equation (17):

$$min U_A\left(x\right)\ge \gamma$$

(17)

However, there is the feature that $U_{\mathrm{A}}\left(x\right)$ is also less than or equal to $U\left(x\right)$ for the same point because of $r_m{\left(x\right)}^p\le 1$. The probability of being correctly classified is equal to $2\varnothing \left(\frac{\gamma }{r_m{\left(x\right)}^p}\right)-1$; in other words, it has a higher probability of accurate classification than the U function. It is a remarkable effect on improving learning accuracy.

3.2. A Semi-Parallel Active Learning Method Based on Kriging

Based on the $U_{\mathrm{A}}\left(x\right)$ aforementioned, a new semi-parallel active learning method based on Kriging is proposed to balance the estimation accuracy and the computational cost. In SPAK, the new training points are added semi-parallel to update the surrogate model. The semi-parallel indicates that the added new training points may be one point or two points during one learning process, which depends on the different positions of learning functions. So SPAK is substantially different from the traditional parallel learning and sequential learning method. The flowchart of SPAK is shown in Figure 1 and the procedure is described as follows:

Step 1. Initialization. Generate a population of $n_{\mathrm{MC}}$ sample points set $S=\left\{x_i:i=1,2,\dots ,n_{\mathrm{MC}}\right\}$ based on PDF in the design space. It presents the candidate training points of the active learning process. The population of $n_{\mathrm{MC}}$ will be changed greater if the condition of updating is satisfied.

Step 2. Define original DoE. In order to build an initial Kriging surrogate model, a DoE contains $n_{\mathrm{TP}}$ samples that are stochastically selected from $n_{\mathrm{MC}}$ or the Latin hypercube sampled in the truncated domain as [μ_x ± 5σ_x]. Then, the performance function is invoked to calculate the responses of these samples. The size of samples in the original DoE usually is small. In the cited paper [4], the number is larger with the number of random variables increased.

Step 3. Construct or refine the Kriging model. The Kriging model is constructed or refined based on the current DoE and responses of these training points. Here, the toolbox of DACE is recommended.

Step 4. Predict using the Kriging model. Kriging predictions $\widehat{G}\left(x_i\right)$, variances ${\widehat{\sigma }}_{\widehat{G}}{}^2\left(x_i\right)$, and correlation coefficient vector $r\left(x_i\right)$ for $i=1,2,\dots ,n_{\mathrm{MC}}$ with Equations (6)–(8) are obtained. Then, the probability of failure is obtained based on Equation (1).

Step 5. Identify new training points in semi-parallel. In this step, the two learning functions in Equations (15) and (16) are calculated for all $n_{\mathrm{MC}}$ sample points, simultaneously. The new potential training points are identified, where $U\left(x\right)$ and $U_{\mathrm{A}}\left(x\right)$ reach the minimum. For most situations, learning functions are primarily based on predictions $\widehat{G}\left(x_i\right)$ and variances ${\widehat{\sigma }}_{\widehat{G}}{}^2\left(x_i\right)$, they are the same point for high probability. So it will only add one new point to the current training points. Oppositely, a point with a low correlation $r\left(x_i\right)$ to the current training points is most likely different from the point based on the U function. Then, both of the two points will be added to the current training points parallelly.

Step 6. The stopping criterion of the learning process. If the best next points are identified with the learning functions, the values of the corresponding learning functions are compared to the stopping criterion in Equation (18).

$$min\left(U_{min}, U_{Amin}\right)\ge \gamma$$

(18)

Step 7. Update DoE with new points. If the stopping criterion in Step 6 is not satisfied, the learning carries on and the new points are evaluated on the performance function. The learning process turns back to Step 3 to refine the new surrogate model with the updated DoE composed of $n_{\mathrm{TP}}+1$ or $n_{\mathrm{TP}}+2$ points.

Step 8. Compute $COV\left({\widehat{P}}_{\mathrm{f}}\right)$. If the stopping criterion in Step 6 is achieved, then the learning is finished and the Kriging surrogate model is precise enough. As a result, it remains to be seen if $n_{\mathrm{MC}}$ is an appropriate size for the MCS. If the size is too small, the probability of failure is trustless. Here, it is acceptable to have a coefficient of variation below 2%. It is estimated as Equation (19):

$$COV\left({\widehat{P}}_{\mathrm{f}}\right)=\sqrt{\frac{1-{\widehat{P}}_{\mathrm{f}}}{{\widehat{P}}_{\mathrm{f}}\times n_{\mathrm{MC}}}}\le 2\%$$

(19)

Step 9. Update Monte Carlo population. If $COV\left({\widehat{P}}_{\mathrm{f}}\right)$ is too large, the Monte Carlo population $n_{\mathrm{MC}}$ will be added, which is generated like in Step 1. Then, SPAK turns to Step 4 to forecast the added new population and the active learning method carries on until the stopping criterion is met again. In order to improve efficiency, the evaluations of previous performance function are continually used.

Step 10. End of SPAK. If $COV\left({\widehat{P}}_{\mathrm{f}}\right)$ is lower than the target, SPAK stops and the last estimation of the failure probability is considered as the result of the performance function.

The new learning function and semi-parallel active learning method based on Kriging are illustrated above, and the learning method can be programed in Matlab 2016b by MathWorks Inc, Natick, MA, USA.

4. Numerical Examples

The performance of a complex system is usually affected by various uncertainties, such as variant loads and material degradation [56]. It is not practicable to directly evaluate the system using brute force MCS because it would take too much time and calculation to study the reliability of a complex system with uncertain factors. Therefore, we usually adopt other more efficient methods. A surrogate model is an effective method to replace complex engineering problems in reliability analysis.

To facilitate comparison with existing methods, in this section, three examples are tested to illustrate the advantage of the SPAK method. These are general examples of reliability analysis from references [5,57,58]. First, a highly nonlinear curve function is tested to observe the process of different methods. Next, the six-hump camel-back (SC) function with two dimensions is evaluated. Low-dimensional examples are easy to display graphically, so they can show the learning processes of training points for each iteration, and reveal the different characteristics of the learning functions. Lastly, a nonlinear and high dimensions analytical function with the dynamic response of a nonlinear oscillator is tested. This example demonstrates the effectiveness of high dimensional nonlinear problems. The failure probabilities of the three examples are also significantly reduced, which increases the difficulty of reliability analysis. Through comparison, it is found that the proposed SPAK has higher efficiency, accuracy, and better practicability. The examples are experimented by Matlab.

4.1. Example 1: Highly Nonlinear Curve Function

The first mathematical example, designed to find the extreme response of $G\left(x\right)$, is employed to demonstrate the difference of the new training points among the typical Kriging surrogate model methods. The highly nonlinear performance function with one-dimensional space, cited in paper [57], is provided as Equation (20).

$$G\left(x\right)=1.0417x^5-13.25x^4+59.792x^3-112.75x^2+75.167x+5$$

(20)

where variable x is within the interval [0, 5]. The global minimum occurs at x = 4.515 with G(4.515) = −2.5036, whereas the local minimum is located at x = 1.957 with G(1.957) = 3.9824, shown in Figure 2.

The performance function is firstly evaluated at initial samples $x_1~x_5$: [0, 1, 2, 2.5, 5], and the obtained performance function values are $G\left(x_1\right)~G\left(x_5\right)$: [5, 15, 4.004, 6.630, 10.148]. With these initial sample points, the one-dimensional surrogate model can be built to approximate the limit state function $G\left(x\right)$. In order to clarify the difference of the modeling process among the methods, the random Kriging method without the learning process is built by 15 points. Based on the learning functions and stopping criterions of different methods, the learning processes are shown as

Table 1. Results of adding training points for different methods for Example 1.

No. of Iteration (NoI)	Kriging	EGO	AK	SPAK
0	0	0	0	0
	1	1	1	1
	2	2	2	2
	2.500	2.500	2.5000	2.500
	5	5	5	5
1	0.4877	3.469	3.6818	3.6818 3.7500
2	0.6349	4.238	3.0951	4.2338 4.3749
3	1.3924	2.975	1.6927	4.6184
4	2.7344	1.541	4.2485	4.7260
5	3.1617	0.457	4.1366	4.7512 3.0909
6	4.0736	4.588	4.4860	4.2054
7	4.5289	3.864	4.6726	/
8	4.5668	2.245	4.2030	/
9	4.7875	4.516	4.7538	/
10	4.8244	/	4.2054	/
Pf	0.10996	0.10996	0.10996	0.10996
errmax (x ∈ [4, 5])	4.4123 × 10−4	26.747 × 10−4	7.1643 × 10−4	4.6021 × 10−4
σ (x ∈ [4, 5])	2.3423 × 10−4	11.6403 × 10−4	2.7163 × 10−4	1.5857 × 10−4

and Figure 3. The EGO achieves convergence in nine iterations, the AK and SPAK require 10 and six iterations, respectively. It should be noted that SPAK adds two training points at the first, second and fifth iterations; thus, reduces the number of iterations. So, SPAK is the most efficient.

In Figure 3, the EGO adds the new points near extreme point with a high probability. It seems more efficient than the AK and SPAK on a global space to identify the extreme value of the performance function. However, for the reliability analysis, under the premise of global modeling, more attention should be paid to the threshold region to determine a critical boundary. Based on the learning function in Equation (15), the AK selects a new point that may be crossed, as shown in Figure 3c. The SPAK considers the correlation between sample points and training points based on Equation (16), the added training points are shown clearly in Figure 3d. Meanwhile, SPAK may add the two new training points in once iterated parallelly; the training points are shown as a plus sign (+) and an asterisk (*) in Figure 3d.

The result of the performance function near the critical region has the most direct effect on the failure probability. Therefore, we choose the extreme point and the vicinity of the critical region as the main area of error analysis. Furthermore, in the interval $x\in \left[4,5\right]$, the performance function occurs twice out-crossing. For the sake of comparison, the residual error $y_{res}$ is determined between the actual and predicted values, as described in Equation (21), and the residual errors of several surrogate model approaches are illustrated in Figure 4.

$$y_{res}\left(x\right)=\mathrm{abs}\left(G\left(x\right)-\widehat{G}\left(x\right)\right)$$

(21)

The EGO has the maximum residual error $y_{res}=26.747\times {10}^{-4}$ and the standard deviation $\sigma =11.6403\times {10}^{-4}$. The AK follows with $7.1643\times {10}^{-4}$ and $2.7163\times {10}^{-4}$. The Kriging and SPAK are similar. However, when the sample points change, the calculation results of the Kriging method are unstable without a learning process. The SPAK performs even better with $4.6021\times {10}^{-4}$ and $1.5857\times {10}^{-4}$. The surrogate model based on SPAK near the critical region is accurate enough. It is verified that SPAK has a good balance between efficiency and accuracy by semi-parallel learning.

4.2. Example 2: Six-Hump Camel Back Function

The six-hump camel back function is used for the second example [58]. The design space of the SC function in this paper is smaller compared to the original function, as shown in Equation (22). The design variables are uniform distribution in the interval $-2\le x_1\le 2,-0.5\le x_2\le 1$. This function has six minimums, one of which is global minimum, which occurs at $x=[-0.09029, 0.7126]$ with $G\left(x\right)=-0.0316$, as indicated by the asterisk in Figure 5. To simulate a small probability of failure, the function is shifted up by one unit.

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

(22)

Fifteen initial training points are generated by uniform random sampling. The surrogate models $\widehat{G}\left(x\right)$ are constructed and refined by different training methods. The training process is shown in Figure 6. For the random Kriging surrogate model, 60 points are used for random sample training, which is shown in Figure 6a. Due to the randomness, the training points are uncontrollable, and the accuracy near the threshold value of the surrogate model without learning process is not guaranteed. The EGO adds the 20 training points mainly near the area around the two extreme points in Figure 6b, and the other areas of training are disregarded. In interval $x_1\in \left[-0.5, 0.5\right],x_2\in \left[0.5, 1\right]$, the AK and SPAK train the surrogate models around the threshold value, AK adds seven points in Figure 6c, and SPAK adds six points in Figure 6d. It is important to emphasize that SPAK adds the training points semi-parallel, the triangle training points are generated based on the U function in sequential, while the plus training points based on the $U_A$ and the asterisk points based on the U function are in parallel. Therefore, SPAK only iterates 14 times with 26 training points. It has significantly fewer iterations than the AK.

Furthermore, to eliminate the contingency of the results caused by the initial samples, 20 groups of random samples are selected for the experiment, and the average values of 20 results are finally taken for comparisons.

Table 2. Results of different methods for Example 2.

Method	NoF	NoI	Pf	Relative err (%)	yres (x1 ∈ [−0.5, 0.5] x2 ∈ [0.5, 1])	σ (x1 ∈ [−0.5, 0.5]) x2 ∈ [0.5, 1])
MCS	1 × 106	/	0.002971	/	/	/
Kriging	60	/	0.002625	11.63	0.27934	0.06587
EGO	35.45	20.45	0.002881	3.036	0.34872	0.09747
AK	38.85	23.85	0.002955	0.5301	0.14960	0.03135
SPAK	39.1	14.2	0.002972	0.0471	0.10461	0.02629

shows the results among random Kriging, AK, EGO, SPAK, and brute MCS. The true failure probability is calculated as 0.002971 based on MCS with 1 × 10⁶ sampling evaluations. The results obtained by the four methods are quite different. The random Kriging method has the maximal relative error 11.63%. Obviously, this method is no longer applicable to the problem of small failure probability. The EGO has error 3.036%, better than the random Kriging, and only needs 35.45 number of function evaluations (NoF) and 20.45 NoI. The NoF based on AK and SPAK is close to the former, 38.85, and the latter, 39.1; however, SPAK has fewer NoI than AK, which only has 14.2 iterations. Iteration determines the number of cycles. There is no doubt that SPAK will be more efficient. In terms of accuracy, SPAK is also the best, with a relative error 0.0471%. The absolute error $y_{res}$ and $\sigma$ of SPAK around the threshold value is reduced, which makes SPAK more accurate. The surrogate model is closer to the real performance function.

4.3. Example 3: Stochastic Linear Oscillation System

The third example, which is a rare failure probability, deals with a stochastic oscillation system [5]. It consists of a nonlinear undamped single degree of freedom (SDoF) system in Figure 7. The performance function is defined as:

$$G\left(x\right)=3\mathrm{S}-\left|\frac{2F}{M{\omega }^2}\sin \left(\frac{\omega t}{2}\right)\right|$$

(23)

where $\omega =\sqrt{\frac{K_1+K_2}{M}}$, and input random vector is given as $x={\left[M,K_1,K_2,S,F,t\right]}^{\mathrm{T}}$. The probabilistic characteristics are summarized in

Table 3. Random variables of the SDoF oscillation system.

Variable	Mean-Value	Variable Coefficient	Distribution
M	1	0.05	Normal
K1	1	0.1	Normal
K2	0.1	0.1	Normal
S	0.5	0.1	Normal
F	0.8	0.15	Normal
t	1	0.15	Normal

.

The results of the example depending on different methods are summarized in

Table 4. Results of different methods for Example 3.

Method	NoF	NoI	Pf	Relative Error (%)
MCS	7 × 106	/	3.9897 × 10−4	/
Kriging	200	/	3.9290 × 10−4	1.52
EGO	24.3	9.3	5.4754 × 10−4	10.42
AK	46.95	31.95	4.0403 × 10−4	1.27
SPAK	43.55	25.5	3.9959 × 10−4	0.154

. The MCS method is used to calculate the failure probability with 7 × 10⁶ evaluations. The result based on EGO is inaccurate, with relative errors of 10.42%, though only requires 24.3 evaluations and 9.3 iterations. The SPAK is able to obtain precise results of the rare failure probability, with relative errors that are only 0.154%, far lower than other methods. The SPAK, adding the training points in semi-parallel with the fewer NoF and NoI than AK, is efficient for the reliability analysis of the SDoF oscillation system. It also shows that the increase of random variables has little influence on the project.

Through the above three examples, it is verified that the proposed method has more efficiency and accuracy than the comparison method for nonlinear and high-dimensional reliability analysis problems. It also indicates that the method will have a wide application value.

5. Application

In addition to general numerical examples, the effectiveness of the proposed method can be verified by practical engineering problems. One aquarium is going to be built in the second basement floor of a commercial building. The acrylic board, which is mostly made of polymethyl methacrylate (PMMA), is used to build the aquarium. To ensure safety in case of fire, according to the national standard Code for Design of Building Fire Protection (GB 50016-2014), it is necessary to conduct a safety and reliability analysis. When PMMA receives heat from the flame, the surface temperature of the board will rise. The board loses its mechanical properties and even combusts, thus, collapses, as shown in Figure 8.

The effective thickness model of the board under a different radiation intensity was obtained by the cone calorimeter experiment under 2 h as shown in Figure 9. The expression is shown as Equation (24). The combustion performance and fire resistance that limit the time of building components with various fire resistance grades are specified in national standard GB 50016-2014.

$$b^{\prime }=b_0-\left(-0.0661q^2+8.4396q-89.654\right)$$

(24)

where, b₀ is the initial thickness of the board. It equals to 200 mm. q is the radiation intensity of flame, and has a mean value of 40 kW/m²; the variable coefficient is 0.15. The other random variables include depth of water, dimensions of aquarium structure, and material properties of acrylic, which are truncated random variables with μ ± 2σ and they are shown in Figure 10 and

Table 5. Random variables of the aquarium.

Variable	Mean-Value	Variable Coefficient	Distribution	Description
AB (mm)	7482	0.01	Truncated Normal	Maximum arc length of acrylic board
R (mm)	19,950	0.1	Truncated Normal	Radius of acrylic board
H (mm)	6050	0.1	Truncated Normal	Depth of water
Q (kW/m2)	40	0.15	Truncated Normal	Intensity of thermal radiation
E (MPa)	3130	0.1	Truncated Normal	Modulus of elasticity
[σb] (MPa)	69.5	0.1	Truncated Normal	Allowable strength limit

.

By establishing a finite element simulation model based on Ansys, the maximum stress ${\sigma }_{\mathrm{m}}$ with different conditions is calculated. Ansys carries out a finite element simulation and returns the calculation results to Matlab to build and update the surrogate model based on SPAK. The similar finite element simulation and software call method are also used in reference [59]. Figure 11 shows the finite element model of the aquarium and the stress distribution of the result at one point of input variables. The performance function can be simplified as Equation (25), when $G>0$, it is safe, otherwise, it is unsafe.

$$G=\left[{\sigma }_b\right]-{\sigma }_{\mathrm{m}}$$

(25)

Through DoE and the active learning process, the stress prediction model is constructed and the stress values ${\sigma }_{\mathrm{m}}$ under the fluctuation of design variables are predicted. Then, the failure probability is calculated by Equations (1) and (25). Because of time constraints, the brute MCS method is not available.

Table 6. Results of different methods of the aquarium.

Method	NoF	NoI	Pf	Pr
Kriging	300	/	5.784 × 10−5	0.99994216
EGO	77	57	3.510 × 10−4	0.99964900
AK	62	42	5.522 × 10−5	0.99994478
SPAK	53	27	5.453 × 10−5	0.99994547

shows the results of the comparison between Kriging, EGO, AK, and SPAK, without brute MCS. Therefore, the true failure probability cannot be obtained by MCS. The feasible solution is to analyze accuracy by comparing the results among these methods. The result shows that the proposed method is more available for reliability analysis. Firstly, the failure probability based on the proposed SPAK is 5.453 × 10⁻⁵, which is close to AK-MCS and Kriging. It indicates the correctness of SPAK. EGO with a large error is no more applicable. Secondly, the NoF and NoI of SPAK are reduced to 53 and 27; they are both smaller than other methods. The evaluation efficiency is increased by at least 14.5% and iteration efficiency is increased by 35.7%. Based on these two observations, the effectiveness of the proposed method is finally testified.

6. Conclusions

Currently, the reliability analysis system is becoming more complex, and the efficiency and accuracy of reliability analysis are also increasing. This paper presents a new algorithm called SPAK that uses Kriging models for small failure probabilities. There are two innovations:

(1)

The method proposes a novel learning function that considers the relationship between the training points and samples. It is called the $U_{\mathrm{A}}$ function. Additionally, it is more advanced than the U function, and while it is based on the U function, it is also different from it. The learning function pays attention to the correlation function, predictions and variance, so it overcomes the shortcoming that the training points are always around the area near the threshold. Based on the $U_{\mathrm{A}}$ function, the other essential areas have been highlighted as supplementary. Therefore, the accuracy of the surrogate model based on SPAK is improved.

(2)

To enhance the efficiency of iteration, the method proposes the semi-parallel learning process. It will add one or two training points intelligently, which is only based on one stopping criterion. This learning process reduces the number of iterations and enhances efficiency.

The new learning function and intelligent semi-parallel added points are the most significant advantages compared to other parallel learning methods and sequential learning methods.

To test this new method, three numerical examples and a nonlinear and a multidimensional application are studied. We obtain the following results:

(1)

The failure probability based on SPAK is more accurate than the Kriging, AK and EGO methods in Example 1 to 3. Especially for Examples 2 and 3, the relative error of the failure probability based on SPAK is only one tenth of the other methods. It verifies the accuracy of this method for the failure probability analysis of the engineering application because the real value of failure probability cannot be obtained for practical engineering problems.

(2)

Due to the semi-parallel learning mechanism, the number of iterations required to construct the surrogate model is significantly reduced. In the application of the reliability analysis of the aquarium, evaluation efficiency is increased by at least 14.5% and the iteration efficiency increased by 35.7%. It shows that the semi-parallel learning method has more advantages in computational efficiency, especially in iterative computation.

(3)

Due to SPAK being an approximate model construction method, it focuses on the analysis of the accuracy around the model threshold ($G=0$). Therefore, this method is more suitable for dealing with critical value problems, such as reliability analysis, but may not be the most suitable for solving extreme value problems.

Based on this work, it can be concluded that the proposed algorithm is valuable for constructing reliability the surrogate model. Further, we can combine IS sampling, subset and other methods with SPAK to improve its application value for more rare failure probability analyses. This research group will also carry out relevant research in the follow-up work.