# A Kriging-Based Active Learning Algorithm for Mechanical Reliability Analysis with Time-Consuming and Nonlinear Response

Authors: Cao Tong; Jian Wang; Jinguo Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

When the reliability analysis of the mechanical products with high nonlinearity and time-consuming response is carried out, there will be the problems of low precision and huge computation using the traditional reliability methods. To solve these issues, the active learning reliability methods have been paid much attention in recent years. It is the key to choose an efficient learning function (such as U, EFF, and ERF). The aim of this study is to further decrease the computation and improve the accuracy of the reliability analysis. Inspired from these learning functions, a new point-selected learning function (called HPF) is proposed to update DOE, and a new point is sequentially added step by step to the DOE. The proposed learning function can consider the features like the sampling density, the probability to be wrongly predicted, and the local and global uncertainty close to the limit state. Based on the stochastic property of the Kriging model, the analytic expression of HPF is deduced by averaging a hybrid indicator throughout the real space. The efficiency of the proposed method is validated by two explicit examples. Finally, the proposed method is applied to the mechanical reliability analysis (involving time-consuming and nonlinear response). By comparing with traditional mechanical reliability methods, the results show that the proposed method can solve the problems of large computation and low precision.


---

## Body

## 1. Introduction

In the mechanical reliability analysis (MRA), when the joint probability density functionfX(x) of the random variable vector X = [X1, X2, …, Xn]T is known, the failure probability can be expressed as(1)Pf=Pgx≤0=∫FfXxdx=∫Ig≤0xfXxdx,where gx is the performance function of the mechanical system. If gx > 0, the system is in the safety state; otherwise, it is in the failure state. Ig≤0x=1,gx≤0,0,gx>0, and the failure domain F and the safety domain S of input random variables are F=xgx≤0 and S=xgx>0, respectively.However, the performance functiongx of the MRA is obtained by solving the numerical model (such as finite element method and finite different method) which is usually time-consuming and highly nonlinear. For this situation, the moment reliability methods [1] are not applicable because the calculated reliability’s accuracy is very low due to ignorance of the Taylor expansion’s high-order items. While the random simulation methods [2, 3] (such as Monte Carlo simulation (MCS), important sampling (IS), subset simulation (SS), line sampling (LS), and Latin hypercube sampling) are unaffordable because their calculation computation is too huge for calling a large of numerical models. It also takes too much time to calculate the reliability by the workstation computer. Therefore, in MRA, it is the core problem to minimize the number of calls to the numerical calculation under the premise of ensuring the calculated accuracy.To solve the problem of large computation in MRA, several surrogate models (including polynomial response surfaces [1], sparse polynomials [4], neural networks [5], support vector machines [6], and Kriging [7]) are usually replaced with the real performance function. Since the speed of calculating the surrogate function is much faster than the real performance function, the surrogate model techniques are widely applied to MRA. However, there are two weaknesses in these surrogate model techniques: firstly, the sample points of the design of experiment (DOE) must be arranged in advance before the surrogate model of the performance function is constructed. This fixed DOE’s scheme will cause excessive sampling of areas where design variables are not important and then results in unnecessary computational waste. Secondly, the traditional surrogate-based reliability methods only guarantee that the fitting performance function is as accurate as possible, while the accuracy of the reliability analysis cannot be guaranteed.Kriging, as a surrogate model, shows better fitting accuracy than other surrogate models [4, 8, 9]. Kriging not only is a interpolation model but also provides the mean and variance at an unknown point. Because of the statistical information of Kriging, the Kriging-based active learning reliability method (KARM) has received much attention in the recent years. The basic idea of KARM is to sequentially refresh the DOE, and the constructed surrogate model gradually reaches to the real performance function. Inspired by the famous efficient global optimization (EGO) method [10], Bichon et al. [11, 12] proposed an efficient global reliability analysis (EGRA) method to search the DOE’s points in the vicinity of the limit state. Echard et al. proposed an active learning reliability method (AK-MCS/AK-IS) [13, 14] by combing MCS/IS and Krigng. Yang et al. [15] proposed a reliability mixed algorithm (ALK-HRA) based on probability and convex set hybrid reliability analysis. Wen et al. [16] proposed a dynamic sampling area strategy (SKRA) to avoid selecting samples in areas with low probability density, to improve the efficiency of selecting the best point. Among these methods, AK-MCS [13] is one of the most typical KARMs, whose learning function U is an innovative point-selected learning function that measures the probability of classifying a point into the wrong domain. Based on AK-MCS, there are many scholars who have extended the study of U in order to solve some special conditions (such as small failure probability, multiple dimension, multiple failure regions, and systems) of reliability analysis. For example, Echard et al. [14] combined the AK-MCS algorithm idea with IS (AK-IS) to make it suitable for solving the reliability problem with small failure probability. Tong et al. [17] combined the AK-MCS idea with the SSIS method for the small failure probability problem, to improve the robustness of multidimensional random variables. Cadini et al. [18] improved the AK-IS to solve the reliability problem of multiple areas and small failure probability. Fauriat and Gayton [19] presented the AK-SYS method which is adaptation to system problems.Many studies [20–23] indicated that it is the key problem to choose a good learning function for guiding the DOE’s strategy of KARM. Some efficient learning functions have been presented in recent years. Yang et al. [20] defines the expected risk function (ERF), which provides an indication of how large the predicted performance function value is if the sign of the predicted performance function is wrongly predicted. Lv et al. [21] proposed a new learning function based on information entropy theory and applied it to the reliability algorithm in combination with the LS method. Sun et al. [22] proposed a learning function that considers the degree of improvement in prediction accuracy and probability density based on a dynamic point selection update strategy and established a relatively high-precision Kriging prediction model relatively quickly. Wang and Sun [23] roughly determined the region where gx = 0 according to Kriging prediction statistics and gradually increased the number of sample points in the region to maximize a certain accuracy index of the Kriging model.The aim of this paper is to further decrease the computation and improve the accuracy of reliability analysis. A Kriging-based active learning reliability method is presented, in which an efficient point-selected learning function (called HPF) is proposed to update DOE. The uncertainty of HPF could measure the features like the sampling density, the probability to be wrongly predicted, and the local and global uncertainties close to the limit state. The new point-selected learning function can help select the next point effectively, and only a small number of real performance function evaluations are required to build an accurate Kriging performance function. Finally, the gear reliability analysis is taken as the object of this study, whose performance function response’s behavior is nonlinear and time-consuming. It is a difficult task that gear reliability is evaluated accurately and efficiently by traditional mechanical reliability methods (such as crude MCS, FORM, and RSM). The proposed method is applied to the reliability analysis of gear. This study will provide reference for the reliability design of the gear vibration.This article is organized as follows: Section2 introduces the structural reliability method based on the Kriging surrogate model and Monte Carlo simulation. Section 3.1 presents a new hybrid point-selected learning function and then describes the analytic expression of HPF. Section 3.2 describes the procedures of the proposed Kriging-based active reliability method. Section 4 illustrates the efficiency and correctness of the proposed method through two academic examples. In Section 5, the proposed method is applied to analyze the reliability analysis of gear. Section 6 is the conclusion.


## 2. Kriging Surrogate Model for Structural Reliability Analysis

Following [4], the Kriging response gx can be expressed as(2)gx=∑ipβikix+zx=kxTβ+zx,where β is the vector of regression coefficient; k(x) is the vector of polynomial, where constant polynomials are selected as the basic functions; and z(x) is the random process, and the covariance of z(x) is(3)Ezwzx=σ2Rθ,w,x,where σ2 is the variance of the process and R(θ, w, x) is the correlation function between z(w) and z(x), where a Gaussian correlation function is selected as(4)Rθ,w,x=∏knRkθk,wk,xk=exp−∑k=1nθkwk−xk2.Under the consideration of linear predictor, the predictor at the unobserved pointx is a normal distribution variable before calculating the true response. And the predictor follows a Gaussian distribution:(5)gx∼Ng^x,sx,where the Kriging predictor value g^x and standard variance sx can be expressed as(6)μx=g^x=kxTβ∗+rxTR−1Y−Fβ∗,(7)s2x=σZ1+uxTFTR−1F−1ux−rxTR−1rx,where β∗ is the regression parameters of generalized least squares solution; R is the matrix of stochastic-process correlations between design sites; and R = [R · ij] m × m, where R · ij = R(θ, si, θj), i, j = 1, ⋅⋅⋅, m; rx=Rθ,s1,xRθ,s2,x…Rθ,sm,xT; ux=FTR−1r−k; σz=Y−Fβ∗TY−Fβ∗/m.The optimal coefficientsθ∗ corresponding to maximum likelihood estimation should solve(8)minθψθ≡R1/mσ2,where |R| is the determinant of R.Whenθ∗ is found, a Kriging prediction model can be established.Then, the Kriging surrogate expressiong^x will be substituted for the performance function gx to calculate the failure probability, and it is used to calculate the failure probability:(9)pf≈p^f,g^=∫Ig^≤0xfXxdx.Becauseg^x is explicit expression, the computer can complete the calculation of the performance function gx at millions of sample points in a short time. Therefore, Monte Carlo simulation, as the most common and robust method, is applied to reliability analysis in this paper. The failure probability can be expressed as(10)pf=1nMCS∑i=1nMCSIgxi=ng≤0nMCS,where x = [x1, x2, ⋅⋅⋅, xn]T is the vector of input random variables with n-dimension in the domain; fX(x) is the joint probability density function of x; nMCS is the number of random samples generated by MCS; xi, (i = 1, 2, … , nMCS) is the ith sample in the Monte Carlo population; Igx is a indicator function, and it can be defined as(11)Igx=0,whengx≤0,1,whengx>0.The coefficient of variation ofp^f,g is calculated as(12)δ=1−p^f,gnMCSp^f,g.


## 3. Kriging-Based Active Learning Reliability Method

### 3.1. Point-Selected Learning Function

The idea of the Kriging-based active learning reliability method (KARM) is to perform a relatively small DOE and sequentially add a new point step by step to the DOE, and then the constructed Kriging surrogate model gradually approaches to the real performance function. So the surrogate model is to be replaced with the real performance function when reliability analysis is conducted. Studies show that the key of the Kriging-based active learning reliability method is to choose an efficient point-selected learning function, where some widely used learning functions consist of such as EFF [11], U [13], and ERF [15]. Inspired by the own advantages of these efficient learning functions, combination of some efficient features is proposed and deduced to further improve the efficiency and accuracy of the KARM. The proposed learning function (called “HPF”) could consider the features like the sampling density, the probability to be wrongly predicted, and the local and global uncertainty close to the limit state. And, the derivation of HPF is elaborated below.Figure 1
Probability density function ofgx.Firstly, based on equation (5) of Kriging theory in Section 2, gx is assumed to a random variable and gx∼Ng^x,sx under the Kriging framework. As shown in Figure 1, fggx is the probability density and g^x and sxare the Kriging predictor value and standard variance, to equation (1), the failure probability pf is only related to the positive or negative of the state function gx instead of the absolute value of gx in the X domain. The accuracy of the Kriging-based reliability analysis can be interpreted as the accuracy of g^x = 0 relative to gx = 0. As shown in Figure 2, g^x=0 and gx = 0 divide the X space into four parts: the signs of g^x and gx in the domain D1 and D3 are the same and the signs of g^x and gx in the domain D2 and D4 are converse. The probability that the sign of g^x is opposite to the sign of gx can be expressed as(13)Error=∫D2∪D4fXxdx=∫Ig≤0x−Ig^≤0xfXxdx.Figure 2
Illustration of the accuracy ofg^x=0.In the Kriging framework, the predicted sign at the pointx is a random variable Ig≤0x, and Ig≤0x obeys the distribution of 0-1 with the parameter Ugx=g^x/sx. According to [24], the accuracy of p^f can be measured as(14)Eg=EError=∫Φ−Ugx⋅fXxdx.It is shown from equation (14) that the larger the value of Φ−Ugx⋅fXx is, the larger is the accuracy of probability of predicting gx to be wrong, where Ugx is the learning function of AK-MCS [13] and fXx indicates the sampling points’ density; so Φ−Ugx⋅fXx will be regarded as one part of the proposed indication function.Secondly, as shown in Figure3(a), the larger the value of RFgx=maxgx−0,0,g^x<0,max0−gx,0,g^x≥0, [15] is, the higher the uncertainty of the corresponding sample point in the global area is. According to [11], the point with the larger value of EFgx=maxε−gx,0 is to be the higher uncertainty in local region −ε,ε, which is shown in Figure 3(b). Therefore, RFgx⋅EFgx could balance the local and global uncertainty of gx, and then it measures the uncertainty more reasonably.Figure 3
Uncertainty of different Kriging-based indications: (a) RF(gx); (b) EF(gx).
(a)
(b)However,RFgx⋅EFgx is a random variable, so the expectation of RFgx⋅EFgx is used to measure this uncertainity. Just like EFF [11] and ERF [15], we average RFgx⋅EFgx throughout the real space. So ERFgx⋅EFgx is expressed as(15)ERFgx⋅EFgx=∫−∞+∞RFgx⋅EFgxfgxgxdgx=∫−ε0maxgx−0,0⋅maxε−gx,0fgxgxdgx∫0εmax0−gx,0⋅maxε−gx,0fgxgxdgx.This integral is expressed analytically as (the derivation process can be seen in Appendix)(16)ERFgx⋅EFgx=−signg^xg^x−εsigng^x−g^xsx⋅Φ−εsigng^x−g^xsx−Φ−g^xsx+g^xsxφ−εsigng^x−g^xsx+sx−εsigng^x−g^xφ−g^xsx,where sign⋅ is the sign function and signg^x=1,g^x≥0,−1,g^x<0.Above all, a hybrid point-selected function (HPF) is defined as(17)HPFx=−signg^xfXxΦ−g^xsxgx−εsigng^x−g^xsxΦ−εsigng^x−g^xsx−Φ−g^xsx+g^xsxφ−εsigng^x−g^xsx+sx−εsigng^x−g^xφ−g^xsx,where g^x and sx are the Kriging predictor value g^x and standard variance sx seen in equations (6) and (7), respectively; Φ(⋅) and φ(⋅) are the cumulative distribution function and the probability density function of the standard normal distribution, respectively; ε=k⋅sx and k = 3 is selected in this paper; sign⋅ is the sign function, and signg^x=1,g^x≥0,−1,g^x<0. Note that equation (17) is hereafter referred to as ‘HPF’ in this study because HPF is the abbreviation of the hybrid point-selected function.


### 3.2. Stopping Criterion

The stopping criterion directly determines the number of calls to the real performance function (referred to as ‘ncall’). If ncall is too small, the accuracy of the constructed surrogate model is too low; while if ncall is too large, the computation is too large. Therefore, it is the key point to choose a precise and efficient criterion in the active learning reliability methods.In AK-MCS [13], the stopping criterion can be expressed as(18)minx∈SUx≥2,where U(x) is the learning function of AK-MCS [13], Ux=g^x/sx, S is the sampling points of MCS, and S=x1,x2,...,xnMCS.The stopping criterionminx∈SUx≥2 only considers the sampling point with the worst predicting accuracy. It does not consider the prediction accuracy of all sampling points of MCS. This stopping criterion is too rough and not rigorous. In order to more accurately predict the sign of all sampling points, a stopping criterion considering global samples is presented in the literature [25], namely,(19)PC≥0.99,where the event C denotes all the points of the MCS population are correctly classified, and PC=∏i=1nMCSΦg^xi/sxi.The stopping criterionPC≤0.99 means that the probability that all the points are correctly classified is greater than 0.99. It requires that the prediction accuracy of every point of the MCS population is more accurate. For example, if the prediction accuracy of every point of the 105 MCS population is assumed as P(Ci) = 0.99999, the probability that all points are correctly classified should be PC=∏i=11000000.99999=0.36788. Or there is only one sampling point with P(Ci) < 0.99; the active learning process will not stop. Therefore, this stopping criterion is so strict that the accuracy of the calculated failure probability is very high, but it also requires more iterations. Especially when the real performance function is time-consuming, the calculation of the active learning reliability method is too large.In order to avoid above problems of the classical stopping criteria, an improved stopping criterion is adopted in this paper. This criterion is based on the fact that the signs of all sampling points are independent, and it could consider the global sampling points as well. This stopping criterion is not affected by a few local points with the low prediction accuracy (such asP(Ci) < 0.99). Therefore, it can get a relatively acceptable results faster than others. The improved stopping condition is presented as (the derivation process and more details can be found in the literature [17])(20)if any ofNck≥PthresholdnMCSPcontrolk,k=1,2,…,nc,then iteration stops,where Pthreshold is the threshold probability of classifying the right sign on S; if Pthreshold = 0.998, it means that the probability of predicting a right sign on population S is 0.998 at least; Pcontrolk(k = 1, 2, ..., nc) are uniformly distributed in [Pthreshold,1), and the larger the nc is, the more accurately the stopping condition (20) will judge the state of population (nc equals to 100 in this paper); Nck=the number ofPi>PcontrolkPi∈P, where P=P1,P2,…,PnMCS is the vector corresponding to S=x1,x2,…,xnMCS, e.g., Pi=ΦUxi.To better express how to use the stopping criterion, the judgment process of the stopping criterion is illustrated as follows. After each updating the Kriging prediction model, the following process is conducted.Firstly, the probabilities of making a right sign at all sampling points are calculated by the current Kriging prediction model, andP=P1,P2,…,PnMCS, Pi = Φ(U(xi))(i = 1, 2, ..., nMCS).Secondly,Pcontrolk(k = 1, 2, ..., nc) are uniformly selected from the interval [Pthreshold, 1), and they are Pcontrol1=Pthreshold, Pcontrol2=Pthreshold+1×1−Pthreshold/nc,…,Pcontrolnc=Pthreshold+nc−1×1−Pthreshold/nc, respectively.Thirdly,Nckk=1,2,…,nc are calculated. For example, in case of k = 1, the number of sample points whose Pi is greater than Pcontrol1, where P=P1,P2,…,PnMCS is the vector corresponding to S=x1,x2,…,xnMCS, e.g., Pi = Φ(xi) (i = 1, 2, ..., nMCS). In the similar way, when k = 2, 3, ..., nc, Nc2, Nc3, ... Ncnc are calculated, respectively.Finally, on the basis ofPcontrolk(k = 1, 2, ..., nc) and Nckk=1,2,…,nc, respectively, judge whether Nc1≥PthresholdnMCS/Pcontrol1?, Nc2≥PthresholdnMCS/Pcontrol2?, ..., Ncnc≥PthresholdnMCS/Pcontrolnc? If any of Nck≥PthresholdnMCS/Pcontrolk,k=1,2,…,nc is meet, then the iteration would stop.


### 3.3. Procedures in the Present Study

The procedure (Figure4) of the proposed active learning reliability method can be described as follows:Step 1. Generate a Monte Carlo population S=x1,x2,…,xnMCS in the design space. At this stage, none of them is evaluated on the real performance function gx, and they are regarded as the candidate points to be evaluated on the real performance if KARM requires it.Step 2. Define the initial DOE SDOE=p1,…,pNT with small size (N≥n+1n+2/2, in which n is the number of random variables over the bounds ± five standard deviations, and calculate YDOE=y1,…,yNT corresponding to SDOE with the real performance function gx.Step 3. Construct an Kriging predicted function model g^x from the SDOE and YDOE.Step 4. Identify the best point x∗ by the proposed HPF which is shown in equation (19). The HPF criterion of finding the best point is expressed as(21)x∗=argmaxx∈SHPFx.Step 5. Predict the performance function value g^xi,i=1,…,nMCS, and estimate the failure probability p^f,g^ and its coefficient of variation δ according to equations (10) and (12), respectively.Step 6. If the stopping condition (equation (20)) is satisfied, stop iteration and go to Step 7; otherwise, N = N + 1, and update the previous DOE with the best point x∗, and go to Step 3.Step 7. If δ≤0.03, stop iteration and go to Step 8; otherwise, the population S enlarges and go to Step 3.Step 8. Output the results.Figure 4
Procedure of the proposed reliability method.


## 3.1. Point-Selected Learning Function

The idea of the Kriging-based active learning reliability method (KARM) is to perform a relatively small DOE and sequentially add a new point step by step to the DOE, and then the constructed Kriging surrogate model gradually approaches to the real performance function. So the surrogate model is to be replaced with the real performance function when reliability analysis is conducted. Studies show that the key of the Kriging-based active learning reliability method is to choose an efficient point-selected learning function, where some widely used learning functions consist of such as EFF [11], U [13], and ERF [15]. Inspired by the own advantages of these efficient learning functions, combination of some efficient features is proposed and deduced to further improve the efficiency and accuracy of the KARM. The proposed learning function (called “HPF”) could consider the features like the sampling density, the probability to be wrongly predicted, and the local and global uncertainty close to the limit state. And, the derivation of HPF is elaborated below.Figure 1
Probability density function ofgx.Firstly, based on equation (5) of Kriging theory in Section 2, gx is assumed to a random variable and gx∼Ng^x,sx under the Kriging framework. As shown in Figure 1, fggx is the probability density and g^x and sxare the Kriging predictor value and standard variance, to equation (1), the failure probability pf is only related to the positive or negative of the state function gx instead of the absolute value of gx in the X domain. The accuracy of the Kriging-based reliability analysis can be interpreted as the accuracy of g^x = 0 relative to gx = 0. As shown in Figure 2, g^x=0 and gx = 0 divide the X space into four parts: the signs of g^x and gx in the domain D1 and D3 are the same and the signs of g^x and gx in the domain D2 and D4 are converse. The probability that the sign of g^x is opposite to the sign of gx can be expressed as(13)Error=∫D2∪D4fXxdx=∫Ig≤0x−Ig^≤0xfXxdx.Figure 2
Illustration of the accuracy ofg^x=0.In the Kriging framework, the predicted sign at the pointx is a random variable Ig≤0x, and Ig≤0x obeys the distribution of 0-1 with the parameter Ugx=g^x/sx. According to [24], the accuracy of p^f can be measured as(14)Eg=EError=∫Φ−Ugx⋅fXxdx.It is shown from equation (14) that the larger the value of Φ−Ugx⋅fXx is, the larger is the accuracy of probability of predicting gx to be wrong, where Ugx is the learning function of AK-MCS [13] and fXx indicates the sampling points’ density; so Φ−Ugx⋅fXx will be regarded as one part of the proposed indication function.Secondly, as shown in Figure3(a), the larger the value of RFgx=maxgx−0,0,g^x<0,max0−gx,0,g^x≥0, [15] is, the higher the uncertainty of the corresponding sample point in the global area is. According to [11], the point with the larger value of EFgx=maxε−gx,0 is to be the higher uncertainty in local region −ε,ε, which is shown in Figure 3(b). Therefore, RFgx⋅EFgx could balance the local and global uncertainty of gx, and then it measures the uncertainty more reasonably.Figure 3
Uncertainty of different Kriging-based indications: (a) RF(gx); (b) EF(gx).
(a)
(b)However,RFgx⋅EFgx is a random variable, so the expectation of RFgx⋅EFgx is used to measure this uncertainity. Just like EFF [11] and ERF [15], we average RFgx⋅EFgx throughout the real space. So ERFgx⋅EFgx is expressed as(15)ERFgx⋅EFgx=∫−∞+∞RFgx⋅EFgxfgxgxdgx=∫−ε0maxgx−0,0⋅maxε−gx,0fgxgxdgx∫0εmax0−gx,0⋅maxε−gx,0fgxgxdgx.This integral is expressed analytically as (the derivation process can be seen in Appendix)(16)ERFgx⋅EFgx=−signg^xg^x−εsigng^x−g^xsx⋅Φ−εsigng^x−g^xsx−Φ−g^xsx+g^xsxφ−εsigng^x−g^xsx+sx−εsigng^x−g^xφ−g^xsx,where sign⋅ is the sign function and signg^x=1,g^x≥0,−1,g^x<0.Above all, a hybrid point-selected function (HPF) is defined as(17)HPFx=−signg^xfXxΦ−g^xsxgx−εsigng^x−g^xsxΦ−εsigng^x−g^xsx−Φ−g^xsx+g^xsxφ−εsigng^x−g^xsx+sx−εsigng^x−g^xφ−g^xsx,where g^x and sx are the Kriging predictor value g^x and standard variance sx seen in equations (6) and (7), respectively; Φ(⋅) and φ(⋅) are the cumulative distribution function and the probability density function of the standard normal distribution, respectively; ε=k⋅sx and k = 3 is selected in this paper; sign⋅ is the sign function, and signg^x=1,g^x≥0,−1,g^x<0. Note that equation (17) is hereafter referred to as ‘HPF’ in this study because HPF is the abbreviation of the hybrid point-selected function.


## 3.2. Stopping Criterion

The stopping criterion directly determines the number of calls to the real performance function (referred to as ‘ncall’). If ncall is too small, the accuracy of the constructed surrogate model is too low; while if ncall is too large, the computation is too large. Therefore, it is the key point to choose a precise and efficient criterion in the active learning reliability methods.In AK-MCS [13], the stopping criterion can be expressed as(18)minx∈SUx≥2,where U(x) is the learning function of AK-MCS [13], Ux=g^x/sx, S is the sampling points of MCS, and S=x1,x2,...,xnMCS.The stopping criterionminx∈SUx≥2 only considers the sampling point with the worst predicting accuracy. It does not consider the prediction accuracy of all sampling points of MCS. This stopping criterion is too rough and not rigorous. In order to more accurately predict the sign of all sampling points, a stopping criterion considering global samples is presented in the literature [25], namely,(19)PC≥0.99,where the event C denotes all the points of the MCS population are correctly classified, and PC=∏i=1nMCSΦg^xi/sxi.The stopping criterionPC≤0.99 means that the probability that all the points are correctly classified is greater than 0.99. It requires that the prediction accuracy of every point of the MCS population is more accurate. For example, if the prediction accuracy of every point of the 105 MCS population is assumed as P(Ci) = 0.99999, the probability that all points are correctly classified should be PC=∏i=11000000.99999=0.36788. Or there is only one sampling point with P(Ci) < 0.99; the active learning process will not stop. Therefore, this stopping criterion is so strict that the accuracy of the calculated failure probability is very high, but it also requires more iterations. Especially when the real performance function is time-consuming, the calculation of the active learning reliability method is too large.In order to avoid above problems of the classical stopping criteria, an improved stopping criterion is adopted in this paper. This criterion is based on the fact that the signs of all sampling points are independent, and it could consider the global sampling points as well. This stopping criterion is not affected by a few local points with the low prediction accuracy (such asP(Ci) < 0.99). Therefore, it can get a relatively acceptable results faster than others. The improved stopping condition is presented as (the derivation process and more details can be found in the literature [17])(20)if any ofNck≥PthresholdnMCSPcontrolk,k=1,2,…,nc,then iteration stops,where Pthreshold is the threshold probability of classifying the right sign on S; if Pthreshold = 0.998, it means that the probability of predicting a right sign on population S is 0.998 at least; Pcontrolk(k = 1, 2, ..., nc) are uniformly distributed in [Pthreshold,1), and the larger the nc is, the more accurately the stopping condition (20) will judge the state of population (nc equals to 100 in this paper); Nck=the number ofPi>PcontrolkPi∈P, where P=P1,P2,…,PnMCS is the vector corresponding to S=x1,x2,…,xnMCS, e.g., Pi=ΦUxi.To better express how to use the stopping criterion, the judgment process of the stopping criterion is illustrated as follows. After each updating the Kriging prediction model, the following process is conducted.Firstly, the probabilities of making a right sign at all sampling points are calculated by the current Kriging prediction model, andP=P1,P2,…,PnMCS, Pi = Φ(U(xi))(i = 1, 2, ..., nMCS).Secondly,Pcontrolk(k = 1, 2, ..., nc) are uniformly selected from the interval [Pthreshold, 1), and they are Pcontrol1=Pthreshold, Pcontrol2=Pthreshold+1×1−Pthreshold/nc,…,Pcontrolnc=Pthreshold+nc−1×1−Pthreshold/nc, respectively.Thirdly,Nckk=1,2,…,nc are calculated. For example, in case of k = 1, the number of sample points whose Pi is greater than Pcontrol1, where P=P1,P2,…,PnMCS is the vector corresponding to S=x1,x2,…,xnMCS, e.g., Pi = Φ(xi) (i = 1, 2, ..., nMCS). In the similar way, when k = 2, 3, ..., nc, Nc2, Nc3, ... Ncnc are calculated, respectively.Finally, on the basis ofPcontrolk(k = 1, 2, ..., nc) and Nckk=1,2,…,nc, respectively, judge whether Nc1≥PthresholdnMCS/Pcontrol1?, Nc2≥PthresholdnMCS/Pcontrol2?, ..., Ncnc≥PthresholdnMCS/Pcontrolnc? If any of Nck≥PthresholdnMCS/Pcontrolk,k=1,2,…,nc is meet, then the iteration would stop.


## 3.3. Procedures in the Present Study

The procedure (Figure4) of the proposed active learning reliability method can be described as follows:Step 1. Generate a Monte Carlo population S=x1,x2,…,xnMCS in the design space. At this stage, none of them is evaluated on the real performance function gx, and they are regarded as the candidate points to be evaluated on the real performance if KARM requires it.Step 2. Define the initial DOE SDOE=p1,…,pNT with small size (N≥n+1n+2/2, in which n is the number of random variables over the bounds ± five standard deviations, and calculate YDOE=y1,…,yNT corresponding to SDOE with the real performance function gx.Step 3. Construct an Kriging predicted function model g^x from the SDOE and YDOE.Step 4. Identify the best point x∗ by the proposed HPF which is shown in equation (19). The HPF criterion of finding the best point is expressed as(21)x∗=argmaxx∈SHPFx.Step 5. Predict the performance function value g^xi,i=1,…,nMCS, and estimate the failure probability p^f,g^ and its coefficient of variation δ according to equations (10) and (12), respectively.Step 6. If the stopping condition (equation (20)) is satisfied, stop iteration and go to Step 7; otherwise, N = N + 1, and update the previous DOE with the best point x∗, and go to Step 3.Step 7. If δ≤0.03, stop iteration and go to Step 8; otherwise, the population S enlarges and go to Step 3.Step 8. Output the results.Figure 4
Procedure of the proposed reliability method.


## 4. Validation

### 4.1. Case 1 (Modified Rastrigin Function)

This explicit example taken from [13, 18] is used to validate the efficiency and correctness of the proposed method. And, its performance function is a modified Rastrigin function, which is highly nonlinear involving nonconvex and nonconnex domains of failure. And, the performance function can be expressed as(22)gx=10−∑i=12xi2−5cos2πxi,where x=x1,x2T is the vector of input random variables and x1 and x2 are standard normal distributed random variables, and they are independent.This example is performed by the active learning reliability method (which can be seen in Section3.3). Firstly, 10 points are sampled by the Latin hypercube sampling method as the initial DOE, and 105 points (called population S) generated by MCS are regarded as the candidate points. The failure probability estimated by crude MCS, which is calling the real performance function, is regarded as the “true” solution (called it p^f,g in this paper). Figure 5 shows the iterative process of the evaluated p^f,g^. As the number of DOE’s points increases, the estimated failure probability p^f,g^ converges to the true one p^f,g. It takes about 281 points to obtain the failure probability p^f,g^ with relative high accuracy. Figure 6 shows the initial DOE’s points, iterative DOE’s points, and the Kriging predicted function g^x after the iteration of the proposed method run 281 times. From Figure 6, the Kriging predicted performance function g^x is almost consistent with gx in the main domain, and most of the sample points’ signs are right (Figure 7).Figure 5
Graph of the failure probability.Figure 6
Predicted performance function and DOE points.Figure 7
The sample points by MCS.In order to demonstrate the efficiency of HPF, the active reliability method is tested with the different learning functions (such as EFF, ERF, U, and HPF). To avoid the effect of other random factors, EFF, ERF, U, and HPF are performed on the same Monte Carlo candidate points and initial DOE. Figure8 shows the iterative graphs of p^f,g^/p^f,g corresponding to ncall by the different point-selected learning functions. By comparison, the convergence speed of HPF is faster than that of other learning functions, and HPF requires the least number of calls to the real performance function than other point-selected learning functions at the same precision level.Figure 8
Iterative graphs of thep^f with different point-selected learning functions.To demonstrate the influences of the different stopping criteria on the reliability results,minx∈SUx≥2, PC≥0.99 and the improved criterion equation (20) (Pthreshold = 0.998) are, respectively, applied in the active learning reliability methods. The results of different methods are listed in Table 1, where p^f is the estimated failure probability; δ is the coefficient of variation for the estimated failure probability; ε is the relative percentage error in comparison with the reference failure probability p^f,g and it can be calculated by ε=p^f−p^f,g/p^f,g; and ncall is the number of calls to the real performance function.Table 1
Comparison of thep^f values with different methods.
Reliability method
n
call
p
^
f(10−2)
δ
ε
Learning function
Stopping criterion
MCS
1 × 105
7.279
0.01128
—
HPF
min
x
∈
S
U
x
≥
2
418
7.279
0.01129
0.00000
P(C) ≥ 0.99
438
7.279
0.01129
0.00000
Improved criterion
380
7.283
0.01128
0.00055
EFF
min
x
∈
S
U
x
≥
2
447
7.278
0.01129
0.00014
P(C) ≥ 0.99
457
7.279
0.01129
0.00000
Improved criterion
430
7.277
0.01129
0.00028
ERF
min
x
∈
S
U
x
≥
2
431
7.279
0.01129
0.00000
P(C) ≥ 0.99
439
7.279
0.01129
0.00000
Improved criterion
418
7.280
0.01129
0.00014
U
min
x
∈
S
U
x
≥
2
433
7.279
0.01129
0.00000
P(C) ≥ 0.99
438
7.279
0.01129
0.00000
Improved criterion
400
7.279
0.01129
0.00000It can be seen from Table1 that the failure probabilities obtained by three iteration termination criteria are all quite accurate. The result of the global stopping criterion P(C) ≥ 0.99 is the most accurate one, and its relative error is almost 0. This is because P(C) ≥ 0.99 is the most strict one. However, the improved stopping criterion requires the least number of calls to the real performance function among three criteria, whose relative errors are all less than 0.0028 (ε≤0.00028). Therefore, the improved stopping condition could get a relatively accurate result faster than others.


### 4.2. Case 2 (Multidimensional Performance Function)

The second example is a multidimensional performance function, which is a classic illustration in [13, 25]. The explicit expression of the performance function is given as(23)gc1,c1,m,r,t1,F1=3r−2F1mω02sinω02t12,where ω0=c1+c2/m. The means and standard deviations of the six random variables are listed in Table 2.Table 2
Mean and standard deviation of all random variables.
Variable
Distribution type
Mean
Standard deviation
M
GAUSS
1
0.05
c
1
GAUSS
1
0.1
c
2
GAUSS
0.1
0.01
R
GAUSS
0.5
0.05
F
1
GAUSS
1
0.2
t
1
GAUSS
1
0.2In order to demonstrate the efficiency of HPF, the active reliability method is tested with the different learning functions (such as EFF, ERF, U, and HPF). EFF, ERF, U, and HPF are performed on the same Monte Carlo candidate points and initial DOE. Figure9 shows the iterative graphs of p^f,g^/p^f,g corresponding to ncall by the different point-selected learning functions. By comparison, the convergence speed of HPF is faster than that of other learning functions, and HPF requires the least number of calls to the real performance function than other point-selected learning functions at the same precision level.Figure 9
Iterative graphs of thep^f with different point-selected learning functions.To demonstrate the influences of the different stopping criteria on the reliability results,minx∈SUx≥2, PC≥0.99 and the improved criterion equation (20) (Pthreshold = 0.998) are, respectively, applied in active learning reliability methods. The results of different methods are listed in Table 3, where p^f is the estimated failure probability; δ is the coefficient of variation for the estimated failure probability; ε is the relative percentage error in comparison with the reference failure probability p^f,g and it can be calculated by ε=pf−p^f,g/p^f,g; and ncall is the number of calls to the real performance function.Table 3
Comparison of thep^f with different methods.
Reliability method
n
call
p
^
f(10−2)
δ
ε
Learning function
Stopping criterion
MCS
1 × 105
2.873
0.01839
—
HPF
min
x
∈
S
U
x
≥
2
58
2.872
0.01839
0.00035
P(C) ≥ 0.99
65
2.873
0.01839
0.00000
Improved criterion
36
2874
0.01838
0.00035
EFF
min
x
∈
S
U
x
≥
2
64
2.870
0.01840
0.00104
P(C) ≥ 0.99
73
2.873
0.01839
0.00000
Improved criterion
42
2.892
0.01832
0.00661
ERF
min
x
∈
S
U
x
≥
2
61
2.871
0.01839
0.00070
P(C) ≥ 0.99
69
2.873
0.01839
0.00000
Improved criterion
41
2.847
0.01847
0.00905
U
min
x
∈
S
U
x
≥
2
68
2.873
0.01839
0.00000
P(C) ≥ 0.99
68
2.873
0.01839
0.00000
Improved criterion
40
2.873
0.01839
0.00000It can be seen from Table3 that the failure probabilities obtained by three iteration termination criteria are all quite accurate. The result of the global stopping criterion P(C) ≥ 0.99 is the most accurate one, and its relative error is almost 0. This is because P(C) ≥ 0.99 is the most strict one. However, the improved stopping criterion requires the least number of calls to the real performance function among three criteria, whose relative errors are all less than 0.00104 (ε≤0.00104). Therefore, the improved stopping condition could get a relatively accurate result faster than others.


## 4.1. Case 1 (Modified Rastrigin Function)

This explicit example taken from [13, 18] is used to validate the efficiency and correctness of the proposed method. And, its performance function is a modified Rastrigin function, which is highly nonlinear involving nonconvex and nonconnex domains of failure. And, the performance function can be expressed as(22)gx=10−∑i=12xi2−5cos2πxi,where x=x1,x2T is the vector of input random variables and x1 and x2 are standard normal distributed random variables, and they are independent.This example is performed by the active learning reliability method (which can be seen in Section3.3). Firstly, 10 points are sampled by the Latin hypercube sampling method as the initial DOE, and 105 points (called population S) generated by MCS are regarded as the candidate points. The failure probability estimated by crude MCS, which is calling the real performance function, is regarded as the “true” solution (called it p^f,g in this paper). Figure 5 shows the iterative process of the evaluated p^f,g^. As the number of DOE’s points increases, the estimated failure probability p^f,g^ converges to the true one p^f,g. It takes about 281 points to obtain the failure probability p^f,g^ with relative high accuracy. Figure 6 shows the initial DOE’s points, iterative DOE’s points, and the Kriging predicted function g^x after the iteration of the proposed method run 281 times. From Figure 6, the Kriging predicted performance function g^x is almost consistent with gx in the main domain, and most of the sample points’ signs are right (Figure 7).Figure 5
Graph of the failure probability.Figure 6
Predicted performance function and DOE points.Figure 7
The sample points by MCS.In order to demonstrate the efficiency of HPF, the active reliability method is tested with the different learning functions (such as EFF, ERF, U, and HPF). To avoid the effect of other random factors, EFF, ERF, U, and HPF are performed on the same Monte Carlo candidate points and initial DOE. Figure8 shows the iterative graphs of p^f,g^/p^f,g corresponding to ncall by the different point-selected learning functions. By comparison, the convergence speed of HPF is faster than that of other learning functions, and HPF requires the least number of calls to the real performance function than other point-selected learning functions at the same precision level.Figure 8
Iterative graphs of thep^f with different point-selected learning functions.To demonstrate the influences of the different stopping criteria on the reliability results,minx∈SUx≥2, PC≥0.99 and the improved criterion equation (20) (Pthreshold = 0.998) are, respectively, applied in the active learning reliability methods. The results of different methods are listed in Table 1, where p^f is the estimated failure probability; δ is the coefficient of variation for the estimated failure probability; ε is the relative percentage error in comparison with the reference failure probability p^f,g and it can be calculated by ε=p^f−p^f,g/p^f,g; and ncall is the number of calls to the real performance function.Table 1
Comparison of thep^f values with different methods.
Reliability method
n
call
p
^
f(10−2)
δ
ε
Learning function
Stopping criterion
MCS
1 × 105
7.279
0.01128
—
HPF
min
x
∈
S
U
x
≥
2
418
7.279
0.01129
0.00000
P(C) ≥ 0.99
438
7.279
0.01129
0.00000
Improved criterion
380
7.283
0.01128
0.00055
EFF
min
x
∈
S
U
x
≥
2
447
7.278
0.01129
0.00014
P(C) ≥ 0.99
457
7.279
0.01129
0.00000
Improved criterion
430
7.277
0.01129
0.00028
ERF
min
x
∈
S
U
x
≥
2
431
7.279
0.01129
0.00000
P(C) ≥ 0.99
439
7.279
0.01129
0.00000
Improved criterion
418
7.280
0.01129
0.00014
U
min
x
∈
S
U
x
≥
2
433
7.279
0.01129
0.00000
P(C) ≥ 0.99
438
7.279
0.01129
0.00000
Improved criterion
400
7.279
0.01129
0.00000It can be seen from Table1 that the failure probabilities obtained by three iteration termination criteria are all quite accurate. The result of the global stopping criterion P(C) ≥ 0.99 is the most accurate one, and its relative error is almost 0. This is because P(C) ≥ 0.99 is the most strict one. However, the improved stopping criterion requires the least number of calls to the real performance function among three criteria, whose relative errors are all less than 0.0028 (ε≤0.00028). Therefore, the improved stopping condition could get a relatively accurate result faster than others.


## 4.2. Case 2 (Multidimensional Performance Function)

The second example is a multidimensional performance function, which is a classic illustration in [13, 25]. The explicit expression of the performance function is given as(23)gc1,c1,m,r,t1,F1=3r−2F1mω02sinω02t12,where ω0=c1+c2/m. The means and standard deviations of the six random variables are listed in Table 2.Table 2
Mean and standard deviation of all random variables.
Variable
Distribution type
Mean
Standard deviation
M
GAUSS
1
0.05
c
1
GAUSS
1
0.1
c
2
GAUSS
0.1
0.01
R
GAUSS
0.5
0.05
F
1
GAUSS
1
0.2
t
1
GAUSS
1
0.2In order to demonstrate the efficiency of HPF, the active reliability method is tested with the different learning functions (such as EFF, ERF, U, and HPF). EFF, ERF, U, and HPF are performed on the same Monte Carlo candidate points and initial DOE. Figure9 shows the iterative graphs of p^f,g^/p^f,g corresponding to ncall by the different point-selected learning functions. By comparison, the convergence speed of HPF is faster than that of other learning functions, and HPF requires the least number of calls to the real performance function than other point-selected learning functions at the same precision level.Figure 9
Iterative graphs of thep^f with different point-selected learning functions.To demonstrate the influences of the different stopping criteria on the reliability results,minx∈SUx≥2, PC≥0.99 and the improved criterion equation (20) (Pthreshold = 0.998) are, respectively, applied in active learning reliability methods. The results of different methods are listed in Table 3, where p^f is the estimated failure probability; δ is the coefficient of variation for the estimated failure probability; ε is the relative percentage error in comparison with the reference failure probability p^f,g and it can be calculated by ε=pf−p^f,g/p^f,g; and ncall is the number of calls to the real performance function.Table 3
Comparison of thep^f with different methods.
Reliability method
n
call
p
^
f(10−2)
δ
ε
Learning function
Stopping criterion
MCS
1 × 105
2.873
0.01839
—
HPF
min
x
∈
S
U
x
≥
2
58
2.872
0.01839
0.00035
P(C) ≥ 0.99
65
2.873
0.01839
0.00000
Improved criterion
36
2874
0.01838
0.00035
EFF
min
x
∈
S
U
x
≥
2
64
2.870
0.01840
0.00104
P(C) ≥ 0.99
73
2.873
0.01839
0.00000
Improved criterion
42
2.892
0.01832
0.00661
ERF
min
x
∈
S
U
x
≥
2
61
2.871
0.01839
0.00070
P(C) ≥ 0.99
69
2.873
0.01839
0.00000
Improved criterion
41
2.847
0.01847
0.00905
U
min
x
∈
S
U
x
≥
2
68
2.873
0.01839
0.00000
P(C) ≥ 0.99
68
2.873
0.01839
0.00000
Improved criterion
40
2.873
0.01839
0.00000It can be seen from Table3 that the failure probabilities obtained by three iteration termination criteria are all quite accurate. The result of the global stopping criterion P(C) ≥ 0.99 is the most accurate one, and its relative error is almost 0. This is because P(C) ≥ 0.99 is the most strict one. However, the improved stopping criterion requires the least number of calls to the real performance function among three criteria, whose relative errors are all less than 0.00104 (ε≤0.00104). Therefore, the improved stopping condition could get a relatively accurate result faster than others.


## 5. Application of Mechanical Reliability Analysis

Gear is one of the most important parts of the mechanical equipment. When the gears are engaged, the dynamic transmission error of the gear pair will be generated due to the random errors (such as manufacturing error and assembly error). Too large dynamic transmission error will cause the mechanical system to produce strong impact, vibration, and noise [26], so it is regarded as a failure mode of the gear transmission. In this section, the nonlinear vibration system (Figure 10) of a spur gear pair is set as the object of this study, and the proposed reliability method is applied to the reliability analysis of the gear vibration.Figure 10
Nonlinear vibration system of a spur gear pair.


### 5.1. Response of the Gear

The differential equation of the gear vibration system is expressed as (the detailed procedure of derivation is seen in [27])(24)100010−111d2ypdt2d2ygdt2d2ydt2+2ξp0ξmp0ξg−ξmg00ξmdypdtdygdtdydt+k110k130k22−k2300k33fpypfgygfmy=00F−d2edt2,where the parameters of equation (24) are expressed as follows: yt=y¯t¯/b, ypt=y¯pt¯/b, ygt=y¯gt¯/b, ξm=cm/2meωn, ξp=cp/2mpωn, ξg=cg/2mgωn, ξmp=cm/2mpωn, ξmg=cm/2mgωn, cm=2ξmmeωn, et=e¯t¯/b, k11=kp/mpωn2, k13=kmt¯/mpωn2, k22=kg/mgωn2, k23=kmt¯/mgωn2, k33=kmt¯/meωn2, and F=F¯/bmeωn2, where mp and mg are the masses of the gear and pinion, respectively, and me is the equivalent mass of the gear pair; kp, kg, and km are the stiffnesses corresponding to mp, mg, and me; F¯ is the gear meshing force; t=t¯ωn is dimensionless time, in which t is the time; and ωn=k0/me is the inherent frequency, where k0 is the average value of the time-varying meshing stiffness between the gear and pinion. The static transmission error et is converted to a sine function with the Fourier series as et=e¯t¯/b=e¯0/b+e¯r/bsinωmt+φ. In the same way, k13t, k23t, and k33t is, respectively, expressed as k13t=kmt¯/mpωn2=mekmt¯/mpk0=me/mp+me/mpεsinωmt+φk, k23t=kmt¯/mgωn2=mekmt¯/mgk0=me/mg+me/mgεsinωmt+φk, k33t=kmt¯/meωn2=kmt¯/k0=1+εsinωmt+φk.The fourth-order Runge–Kutta numerical integration method edited by Matlab sofware is used to solve the vibration differential equation, and then the response of gear pair is obtained. It is noted that the process of computing the equation (24) is time-consuming and it requires about 7.260 seconds to calculate the implicit performance function. Furthermore, the response of this gear system is highly nonlinear.


### 5.2. Limit State Function of Gear Vibration

The meshing stiffnessε, meshing damping ξm, tooth side clearance bm, and rotational speed ωm of the gear vibration system are considered as random factors. According to [26], ε, ξm, bm, and ωm obey the Gauss distribution, and the means and standard deviations of all random variables are listed in Table 4.Table 4
Mean and standard deviation of all variables.
Variable
Distribution type
Mean
Standard deviation
ε
Gauss
0.2
0.02
ξ
m
Gauss
0.1
0.02
b
m
Gauss
0.07
0.01
ω
m
Gauss
0.3
0.05If the fluctuation of the dynamic transmission errory(t) exceeds the threshold value, the gear vibration system is in the failure state. Based on this failure mode, the performance function of the gear vibration system is expressed as(25)gx=y−Δy,where Δy is the fluctuation value of the dynamic transmission error, i.e., Δy=ymax−ymin, in which ymax and ymin, respectively, represent the maximum and minimum values of the dynamic transmission error during the engagement of the gear pair, namely, ymax=maxyt and ymin=minyt and y is the threshold value, which is selected as 0.055 mm in this paper. When y>Δy, gx>0, so the gear is in the safe state. When y≤Δy,gx≤0, and then the gear is in the failure state.It is known from Section5.1 that y(t) is an implicit response and needs to be obtained by solving equation (25) with the fourth-order Runge–Kutta numerical integration method, so the performance function gx is also an implicit response as well, and each calculation of the gear performance function requires about 7.260 seconds (i3 cpu). Furthermore, the response of the gear performance function is highly nonlinear. For example, Figure 10 is the bifurcation diagram of the gear vibration system’s response corresponding to the tooth side clearance bm. It is shown from Figure 11 that as the tooth side clearance bm increases; firstly, the system’s motion changes from single period (bm ∈ [0.01, 0.032]) to five times period (bm ∈ [0.32, 0.061]), and then it becomes the quasiperiodic motion (bm ∈ [0.061, 0.068]); finally, it is to be the chaotic motion state. Namely, the motion behavior of this gear system is highly nonlinear. And, the performance function response of the gear vibration is very sensitive to the slight changes of the system’s parameters.Figure 11
Bifurcation diagram.Therefore, this is the type problem of mechanical reliability analysis which is time-consuming and nonlinear. Traditional mechanical reliability methods (such as MCS, FORM, and RSM) will become nonapplied because of the problems of large computation and low precision.


### 5.3. Results and Discussion

To validate the efficiency and correctness of the proposed method, the gear vibration reliability analysis is also performed by the traditional mechanical reliability methods (such as crude MCS, FORM [1], and RSM + MCS [5]) and the active learning reliability methods (such as HPF + minx∈SUx≥2, HPF + PC≥0.99, and HPF + improved condition (Pthreshold = 0.99)), where the result of the crude MCS is regarded as the true one; FORM is the first order second-moment reliability method; RSM + MCS means that the polynomial response surface based on Latin DOE is established, and then the constructed response surface is used to calculate the failure probability using MCS.Table5 shows the results of different methods. It can be found from Table 5 that FORM, RSM + MCS, HPF +minx∈SUx≥2, HPF+PC≥0.99, and HPF + improved condition cost the cuputime about 473.1 s, 4150 s, 2412 s, 3455 s, and 1827 s, respectively. The computational efficiency of FORM, RSM + MCS, and the active learning reliability methods (HPF+minx∈SUx≥2, HPF+PC≥0.99, and HPF + improved condition) can be acceptable in the engineering. But the results of RSM + MCS and FORM are totally different from those of the crude MCS, and their relative errors are 0.9693 and 0.8831, respectively, due to the highly nonlinear performance function. Table 5 also shows that the results of the active learning reliability methods are all consistent with that of crude MCS, and the number of calls to the real performance function require 303(15 + 288), 416(15 + 401), and 240(15 + 225), respectively. And, they cost only 2412 s, 3455 s, and 1827 s cuputime. In comparison with 7.264 × 105 s of the crude MCS, the proposed active learning reliability method can improve the efficiency and accuracy dramatically. Comparing with the two stopping criteria minx∈SUx≥2 and PC≥0.99, the improved stopping criterion requires the least number of calls to the real performance function among three criteria, and its accuracy is enough in the engineering (ε=0.0060). The improved stopping condition could get a relatively accurate result faster than others. And, it is more suitable for reliability analysis with a time-consuming performance function.Table 5
Comparison of different methods.
Method
n
call
p
^
f
ε
CPU time
MCS
1 × 105
1.406 × 10−2 (δ = 0.0264)
—
7.264 × 105 s
FORM [1]
45
4.339 × 10−4 (β = 3.33)
0.9693
4.731 × 102 s
RSM + MCS [5]
500
1.650 × 10−3
0.8832
4.150 × 103 s
HPF+minx∈SUx≥2
303
1.404 × 10−2
0.0014
2.412 × 103 s
HPF+PC≥0.99
416
1.407 × 10−2
0.0007
3.455 × 103 s
HPF + improved condition (the proposed method)
240
1.409 × 10−2
0.0060
1.827 × 103 s


## 5.1. Response of the Gear

The differential equation of the gear vibration system is expressed as (the detailed procedure of derivation is seen in [27])(24)100010−111d2ypdt2d2ygdt2d2ydt2+2ξp0ξmp0ξg−ξmg00ξmdypdtdygdtdydt+k110k130k22−k2300k33fpypfgygfmy=00F−d2edt2,where the parameters of equation (24) are expressed as follows: yt=y¯t¯/b, ypt=y¯pt¯/b, ygt=y¯gt¯/b, ξm=cm/2meωn, ξp=cp/2mpωn, ξg=cg/2mgωn, ξmp=cm/2mpωn, ξmg=cm/2mgωn, cm=2ξmmeωn, et=e¯t¯/b, k11=kp/mpωn2, k13=kmt¯/mpωn2, k22=kg/mgωn2, k23=kmt¯/mgωn2, k33=kmt¯/meωn2, and F=F¯/bmeωn2, where mp and mg are the masses of the gear and pinion, respectively, and me is the equivalent mass of the gear pair; kp, kg, and km are the stiffnesses corresponding to mp, mg, and me; F¯ is the gear meshing force; t=t¯ωn is dimensionless time, in which t is the time; and ωn=k0/me is the inherent frequency, where k0 is the average value of the time-varying meshing stiffness between the gear and pinion. The static transmission error et is converted to a sine function with the Fourier series as et=e¯t¯/b=e¯0/b+e¯r/bsinωmt+φ. In the same way, k13t, k23t, and k33t is, respectively, expressed as k13t=kmt¯/mpωn2=mekmt¯/mpk0=me/mp+me/mpεsinωmt+φk, k23t=kmt¯/mgωn2=mekmt¯/mgk0=me/mg+me/mgεsinωmt+φk, k33t=kmt¯/meωn2=kmt¯/k0=1+εsinωmt+φk.The fourth-order Runge–Kutta numerical integration method edited by Matlab sofware is used to solve the vibration differential equation, and then the response of gear pair is obtained. It is noted that the process of computing the equation (24) is time-consuming and it requires about 7.260 seconds to calculate the implicit performance function. Furthermore, the response of this gear system is highly nonlinear.


## 5.2. Limit State Function of Gear Vibration

The meshing stiffnessε, meshing damping ξm, tooth side clearance bm, and rotational speed ωm of the gear vibration system are considered as random factors. According to [26], ε, ξm, bm, and ωm obey the Gauss distribution, and the means and standard deviations of all random variables are listed in Table 4.Table 4
Mean and standard deviation of all variables.
Variable
Distribution type
Mean
Standard deviation
ε
Gauss
0.2
0.02
ξ
m
Gauss
0.1
0.02
b
m
Gauss
0.07
0.01
ω
m
Gauss
0.3
0.05If the fluctuation of the dynamic transmission errory(t) exceeds the threshold value, the gear vibration system is in the failure state. Based on this failure mode, the performance function of the gear vibration system is expressed as(25)gx=y−Δy,where Δy is the fluctuation value of the dynamic transmission error, i.e., Δy=ymax−ymin, in which ymax and ymin, respectively, represent the maximum and minimum values of the dynamic transmission error during the engagement of the gear pair, namely, ymax=maxyt and ymin=minyt and y is the threshold value, which is selected as 0.055 mm in this paper. When y>Δy, gx>0, so the gear is in the safe state. When y≤Δy,gx≤0, and then the gear is in the failure state.It is known from Section5.1 that y(t) is an implicit response and needs to be obtained by solving equation (25) with the fourth-order Runge–Kutta numerical integration method, so the performance function gx is also an implicit response as well, and each calculation of the gear performance function requires about 7.260 seconds (i3 cpu). Furthermore, the response of the gear performance function is highly nonlinear. For example, Figure 10 is the bifurcation diagram of the gear vibration system’s response corresponding to the tooth side clearance bm. It is shown from Figure 11 that as the tooth side clearance bm increases; firstly, the system’s motion changes from single period (bm ∈ [0.01, 0.032]) to five times period (bm ∈ [0.32, 0.061]), and then it becomes the quasiperiodic motion (bm ∈ [0.061, 0.068]); finally, it is to be the chaotic motion state. Namely, the motion behavior of this gear system is highly nonlinear. And, the performance function response of the gear vibration is very sensitive to the slight changes of the system’s parameters.Figure 11
Bifurcation diagram.Therefore, this is the type problem of mechanical reliability analysis which is time-consuming and nonlinear. Traditional mechanical reliability methods (such as MCS, FORM, and RSM) will become nonapplied because of the problems of large computation and low precision.


## 5.3. Results and Discussion

To validate the efficiency and correctness of the proposed method, the gear vibration reliability analysis is also performed by the traditional mechanical reliability methods (such as crude MCS, FORM [1], and RSM + MCS [5]) and the active learning reliability methods (such as HPF + minx∈SUx≥2, HPF + PC≥0.99, and HPF + improved condition (Pthreshold = 0.99)), where the result of the crude MCS is regarded as the true one; FORM is the first order second-moment reliability method; RSM + MCS means that the polynomial response surface based on Latin DOE is established, and then the constructed response surface is used to calculate the failure probability using MCS.Table5 shows the results of different methods. It can be found from Table 5 that FORM, RSM + MCS, HPF +minx∈SUx≥2, HPF+PC≥0.99, and HPF + improved condition cost the cuputime about 473.1 s, 4150 s, 2412 s, 3455 s, and 1827 s, respectively. The computational efficiency of FORM, RSM + MCS, and the active learning reliability methods (HPF+minx∈SUx≥2, HPF+PC≥0.99, and HPF + improved condition) can be acceptable in the engineering. But the results of RSM + MCS and FORM are totally different from those of the crude MCS, and their relative errors are 0.9693 and 0.8831, respectively, due to the highly nonlinear performance function. Table 5 also shows that the results of the active learning reliability methods are all consistent with that of crude MCS, and the number of calls to the real performance function require 303(15 + 288), 416(15 + 401), and 240(15 + 225), respectively. And, they cost only 2412 s, 3455 s, and 1827 s cuputime. In comparison with 7.264 × 105 s of the crude MCS, the proposed active learning reliability method can improve the efficiency and accuracy dramatically. Comparing with the two stopping criteria minx∈SUx≥2 and PC≥0.99, the improved stopping criterion requires the least number of calls to the real performance function among three criteria, and its accuracy is enough in the engineering (ε=0.0060). The improved stopping condition could get a relatively accurate result faster than others. And, it is more suitable for reliability analysis with a time-consuming performance function.Table 5
Comparison of different methods.
Method
n
call
p
^
f
ε
CPU time
MCS
1 × 105
1.406 × 10−2 (δ = 0.0264)
—
7.264 × 105 s
FORM [1]
45
4.339 × 10−4 (β = 3.33)
0.9693
4.731 × 102 s
RSM + MCS [5]
500
1.650 × 10−3
0.8832
4.150 × 103 s
HPF+minx∈SUx≥2
303
1.404 × 10−2
0.0014
2.412 × 103 s
HPF+PC≥0.99
416
1.407 × 10−2
0.0007
3.455 × 103 s
HPF + improved condition (the proposed method)
240
1.409 × 10−2
0.0060
1.827 × 103 s


## 6. Conclusions

The aim of this study is to solve the problem of large computation and low precision for the mechanical reliability analysis which is time-consuming and nonlinear response. A Kriging-based active learning reliability method is presented in this paper. The efficiency of the Kriging-based active learning reliability method mainly depends on the proposed hybrid point-selected learning function (HPF) which is inspired by EFF, ERF, and U. The analytic expression of HPF is deduced by averaging the hybrid indicator throughout the real space. The uncertainty of HPF could measure the features like the sampling density, the probability to be wrongly predicted, and the local and global uncertainty close to the limit state. The new point-selected criterion can help select the next point effectively, and only a small number of true function evaluations are required to build an accurate Kriging performance function. The efficiency and correctness of HPF are validated by two explicit examples. Finally, the proposed method is applied to reliability analysis of the gear vibration (involving time-consuming and highly nonlinear response). By comparing with traditional mechanical reliability methods, the results show that the proposed method can solve the problems of large computation and low precision.For many engineering problems, the response of the system is often calculated by the numerical method (time-consuming) and its behavior is usually highly nonlinear. The proposed active reliability method is independent of the gear reliability analysis and treats the system response as a black box. Thus, the proposed method can be also used to solve other engineering problems to perform reliability analysis.



---

*Source: 7672623-2019-08-18.xml*