A novel reliability analysis methodology based on IPSO-MCopula model for gears with multiple failure modes

In order to accurately and quickly predict the failure probability of gears with multiple failure modes, a novel reliability analysis methodology based on the mixed Copula (MCopula) function model is proposed to deal with the complex correlation among different failure functions. Firstly, we construct a novel MCopula model based on three famous Copula functions: Gumbel Copula, Clayton Copula, and Frank Copula functions. Secondly, we use and improve the particle swarm optimization (PSO) method to optimally calculate the weight coefficients in the proposed MCopula model. Thirdly, the maximum likelihood estimation (MLE) method is adopted to estimate related parameters in the proposed MCopula model. Finally, we verify the proposed reliability analysis methodology with a standard life-prediction case of a strut system and a practical life-problem case of a gear pair system. Comparison results of both cases show that, by using the proposed methodology, the failure probability of a gear pair system with multiple failure correlations can be quickly calculated through a small number of samples and can be estimated as accurately as that by the Monte Carlo scheme. Consequently, our proposed novel methodology successfully analyzes the reliability problems for a gear pair system with multiple failure modes. The proposed methodology can be further applied to solve the reliability problem for other machine parts.


Introduction
As an important transmission component, the metal gear is widely and vastly used in automation machines, which is one of the most common failure parts encountered in a mechanical system. 1 Therefore, it is of great significance to study the issue of analyzing and further improving the reliability of gears so as to help extending the service life of whole machine.Traditional reliability theories of mechanical systems or parts were based on the so called ''failure modes.''And, it was usually assumed that different failure modes of a machine part are independent with each other.However, in practical applications, it was found that different failure modes of machine parts may be affected by the same factors, such as the load, temperature, and material.][10][11][12][13][14][15] Ditlevsen 8 first proposed the narrow bound method for system reliability calculations using some widely accepted failure mode correlations.However, this method becomes very complex as the failure mode increases, and the calculation results are limited for factors within a small interval.Cui and Blockley 16 proposed a reliability analysis method based on the interval probability theory, which had a narrower range than the previous narrow bound method.On the basis of Ditlevsen's theory, Yu et al. 17 proposed a reliability analysis method which uses the linear coefficients scheme to establish the correlation between the primary and the secondary failure modes.And, they verified their proposed method by taking a transmission shaft as an example.Zhou et al. 18 proposed a nonempirical model based on the Nataf transformation to quantitatively analyze the common-cause failure probability of load-dependent systems, which transformed one complex correlation problem into two independent problems.Zhu et al. 19 proposed a combined machine learning strategy to evaluate the fatigue of turbine blade disks.Luo et al. 20 proposed a method combining enhanced uniform importance sampling coupled and support vector regression (EUIS-SVR) for structural reliability analysis with low failure probability.In the above investigations, most of the proposed failure models adopted the assumption that a linear correlation exits among reliability variables.However, these simple models will cause large calculation error of failure probability in structural reliability analysis since the actual relationship among reliability variables are far more complex for this kind of problem.
2][23] Tang et al. 24 used Copula functions to analyze the reliability of a system's structure.Mi et al. 25 proposed a reliability evaluation model which specifically considered the effect of strength degradation, and established the correlation among different failure modes via Copula functions.Sun et al. 21used the Copula function method to capture the dependent relations among residual life distributions, and performed some tests on a microwave machine to verify their proposed method.Jiang et al. 26 analyzed the reliability of a machine structure based on the vine-Copula function.Their proposed method has been proven to be effective in dealing with complex multi-dimensional correlation problems.Although the single Copula function method introduced above has better accuracy and is more efficient than the traditional linear coefficient method, it is still necessary to find more comprehensive reliability analysis methods when facing highly non-linear problems of finding the failure mode correlations for the complex machine parts such as gears.
To better solve the reliability analysis problem of transmission gears, we propose a novel methodology of the integral mixed Copula (MCopula) function with multiple-failure-correlation modes (IMCMFCM).In this proposed IMCMFCM methodology, we majorly adopt the maximum likelihood estimation (MLE) method and together with the improved particle swarm optimization (IPSO) method to estimate the related parameters appeared in the MCopula function.Besides, to identify the performance of our proposed IMCMFCM methodology, some theoretical as well as practical reliability analysis cases are carried out and the results are compared with those obtained by the standard PSO (SPSO) algorithm based MCopula function, the linearly decreasing inertia weight PSO (LDIPSO) algorithm based MCopula function, 27 and the exponential reduction inertia weight PSO (e 2 PSO) algorithm based MCopula function. 28

Copula function
The Copula function is a joint distribution function including a marginal distribution function, which may be used to provide a non-linear correlation between variables.The fundamental theory about the Copula function is given as follows. 29upposed that an N-dimensional random variable X 1 , X 2 , :::, X N ð Þis known, then there exist a Copula function C u 1 , u 2 , :: , such that for any x 1 , x 2 , :::, x N 2 R, the following formula holds: The probability form of the Copula function is expressed as: Eventually, the failure probability of a system with k failure modes can be expressed as 30 : 1łj\h\t łk C P j ,P h ,P t À Á À ::: À ( À 1) kÀ1 C P 1 ,P 2 ,:: where f j x ð Þ is the j th failure function and P j is the failure probability of the j th failure mode.

MCopula function
Among the known Copula functions, the Archimedean Copula function was especially popular because of its excellent reliability prediction ability for machine parts.And, the Gumbel Copula, Clayton Copula, and Frank Copula functions were three commonly used Archimedean Copula functions. 31Different Copula functions can describe different tail features of random variables.The Gumbel Copula function is very sensitive to the upper tail feature of random variables.On the contrary, the Clayton Copula function is sensitive to the lower tail feature of random variables.The Frank Copula function is suitable for describing the symmetric tail feature of random variables. 32Usually, a single Copula function is valid for describing only partial characteristics of random variables.However, it can neither be used to describe the whole characteristics of random variables with simple correlations nor to or describe the partial characteristic of random variable with complex correlations.To overcome this weakness, we now consider a linear combination of different Copula functions which is a new kind of MCopula function, denoted as C M , and is defined as follows: where C k u 1 , u 2 , :::, u N ; u k ð Þis the k th preselected single Copula function,u k is the linking coefficient, and l k is the weight coefficient that satisfies 0łl k ł1, It is noted that different Copula functions may be used to correlate random variables of reliability that have different characteristics.Hence, based on the requirement of high prediction performance in the reliability analysis, we carefully choose three commonly used Copula functions: the binary Gumbel Copula function (C 1 ), the Clayton Copula function (C 2 ), and the Frank Copula function (C 3 ) to construct the novel MCopula function (C M ), as follows:

Parameters estimation
To optimally determine the parameters appeared in the proposed MCopula function, we firstly use the Monte Carlo (MC) scheme to randomly generate the sample data and the MLE scheme to estimate the dependent parameters, where both schemes were quite successfully applied in reliability analysis problems. 3,33,34Secondly, the weight coefficient is determined by the IPSO algorithm proposed by this study.Details are described as follows.
Sample data generation via the MC scheme All random variables (X i ) in the marginal distribution of each failure mode (F i (x i )) in the concerned structure system are sampled by the MC scheme according to their distribution characteristics.
Determination of l k and u via the maximum likelihood scheme Supposed that the marginal distribution functions of the random variables X 1 and X 2 are F 1 x 1 ð Þ and F 2 x 2 ð Þ, respectively.And, we assume that the Copula function is Then, the joint distribution function of the random variables X 1 and X 2 can be expressed as: Further, the joint density function of F(x 1 , x 2 ; u) can be obtained as: The likelihood function of the samples X 1j , X 2j À Á is : The logarithmic likelihood function are: The optimal value of u happens at where The optimal weight coefficient in the MCopula function can be obtained by minimize the Euclidean distance between the empirical Copula function (C(Á)) and the MCopula function (C M (Á)) which is defined by where n is the number of samples.The constructed optimization problem is shown in equation ( 13) and can be solved by the least-square method (LS).minO LS s:t: where l k is the weight coefficient of the MCopula function.
The above optimization equation is a constrained nonlinear equation.By introducing the penalty function, we may transform this constrained problem into an unconstrained one.Eventually, we have the following transformed objective function as: H(q i (l)) = 0; q i (l) ł 0:00001 10; 0:00001 ł q i (l) ł 0:001 50; otherwise where F is the objective function to be optimized, h is the penalty coefficient, q i (Á) is the relative constrained penalty function, and H(Á) is a piecewise penalty function.
Besides, to improve the solution accuracy of the MCopula function, we adopt the known MLE method and IPSO algorithm to find the optimal u and l which play the deterministic roles in calculation accuracy.

Test of goodness of fit t
In order to estimate the accuracy of the proposed bivariate Copula function, we adopt a commonly-used index of goodness of fit.Among published methods of K-S, AIC, and the minimum distance for calculate the goodness of fit, 4,26,32 we select the popular AIC criterion and the minimum distance schemes as the testing methods, as follows: where F ei is the empirical joint probability function of n-dimensional random variables, F ci is the joint probability function of the MCopula function, K is the number of parameters in MCopula function, C Á ð Þ is the empirical copula function, C M Á ð Þ is the MCopula function.To note that the smaller the values of AIC and d 2 , the higher the goodness of fit of the MCopula function.

IPSO
In the traditional PSO algorithm for optimization calculations, 35 all particles with adaptive velocities and positions follow the current optimal particle to search the target in space.The velocity and position of the i th particle in the n-dimensional searching space are Vi = v i1 , v i2 , :::, v in ð Þ and Qi = q i1 , q i2 , :::, q in ð Þ , respectively.The particle updates the velocity and position according to the following formula: where w is the inertia weight, c 1 and c 2 are the learning factors, r 1 and r 2 are the random numbers between 0 and 1, p id is the best position of i th particle in the d th iteration, and p gd is the optimal solution of the whole population in the d th iteration.This traditional PSO algorithm has the defects of low searching accuracy and slow convergence.In order to overcome the shortcomings of the traditional PSO algorithm, here we make two modifications of traditional PSO as follows and we call this as an improved particle swarm optimization algorithm.
Improvement of inertia weight in calculation.It is known that larger inertia weight can raise the global searching ability of the PSO algorithm, while a smaller inertia weight can enhance the local searching ability of the PSO algorithm.An improper inertia weight may easily cause the PSO calculation falling into premature convergence and occurring oscillation near the point where the global optimal solution happens.To solve this problem, we introduce an exponential decreasing inertia weight coefficient as following: where t is the current iteration number and T is the total iteration numbers.
Improvement of learning factor in calculation.A proper learning factor can enhance the self-learning speed in calculations and help for quickly converging to the optimal solution.Therefore, we define a new asynchronous learning factor as follows.
where c i, min and c i, max means the minimum and the maximum of asynchronous learning factors c 1 and c 2 , respectively.
Through the combination of the asynchronous changes in the learning factor and the weight coefficient, the IPSO algorithm has a stronger global searching ability than before in the initial stage.And, it also has the ability of fast convergence to the global optimal solution in the later stage.

IMCMFCM methodology
Our proposed IMCMFCM methodology contains five major steps, as illustrated in the followings.
Step1: Sample data generation.Firstly, the MC scheme is used to generate the sample data.Then, based on these data, we define a cumulative distribution function by the kernel density estimation scheme.

Comparison of the IPSO algorithm with other algorithms
To confirm the performance of our proposed IPSO algorithm in searching global optimal values, we test the calculation performance of the IPSO, SPSO (inertia weight w = 0:5, learning factor c 1 = c 2 = 2), LDIWPSO 27 and e 2 PSO 28 algorithms by using the standard test functions as shown in Table 1 where f 1 (x) and f 2 (x) are the unimodal functions, and f 3 (x) and f 4 (x) are the multimodal functions. 36In calculations with the software package of MATLAB R2020b, we set the dimension of the test function as 30, the size of initial population particles as 40, and the maximum number of iterations as 1000.The convergence conditions of using Sphere, Rosenbrock, Griewank, and Ackley models are respectively shown in Figure 2(a) to (d).
It is seen form Figure 2 that the calculation by the IPSO algorithm converges quicker and the best fitness is smaller than those by other three algorithms.The smaller the best fitness value, the better the optimization effect of the algorithm.Besides, in order to gain further insight into the performance of our proposed IPSO algorithm, we use the mean and standard deviation as the performance indicators, and adopt a maximum of 1000 iterations in each calculation.Considering the randomness of the algorithms, each algorithm was repeated 20 times, and the average value and the standard deviation of the best fitness was obtained, as shown in Table 2.The test results show that, for all function models, the standard deviation obtained by using the IPSO algorithm are much smaller and more accurate than those obtained by other algorithms.We can conclude that the IPSO algorithm performs well not only for multimodal functions, but also for unimodal functions.
Next, we further test the proposed IPSO algorithm by examining the best fitness as the population size changes.In these tests, we still adopt the same four kinds of model functions in IPSO algorithm as previously done, which are Sphere, Rosenbrock, Griewank, and Ackley models.We set a limitation of 1000 maximum iterations in calculations.Considering the randomness of the algorithms, each algorithm was repeated 20 times, the average value of the best fitness obtained by IPSO with using four different models are shown in Table 3.It is found that, for different adopted model functions, using the IPSO algorithm in calculations may always make the fitness correctly converge to its ideal value regardless of the population size.However, as the population size increases, the convergence speed of average fitness reduces.Obviously, a proper choice of population size is quite important in best fitness calculations.As a consequence, we have carried out the test and verified that the IPSO algorithm has good optimization performance.

Comparison of IPSO-MCopula model with other models for a strut system
We now verify the performance of the proposed IPSO-MCopula model through the example of a hollow metal strut subjected to forces at two ends 3 (as shown in Figure 3) where the failure functions of the stability failure mode (g 1 ) and the strength failure mode (g 2 ) are respectively expressed as: and In equations ( 19) and (20), s s is the yielding strength of material, P is the axial load, d is the inner diameter of the strut, c is the wall thickness of the strut, E is the elastic modulus of material, and l is the strut length.Each of the above parameters is assumed to have a normal distribution and specified as listed in Table 4.
In our proposed IMCMFCM methodology, firstly we adopt the MC scheme to generate 2000 sample data, and the IPSO method with a population size of 80 and  36

Mathematical expression Range of search
Sphere a maximum number of iterations of 1000 is used to optimize the algorithm parameters in calculation.Then, we employ the improved fourth-order moment method to calculate the failure probability p f i for each failure mode.Through the manipulation of Step 2 to 5, mentioned in Section 3, the failure probability of the strut system is calculated and the obtained results are shown in Tables 5 and 6.In addition, the sample of  MC method is 10 8 , and the result obtained is used as the truth value to test other algorithms.
For the calculation of failure probability (p f ), our proposed IPSO-MCopula function model exhibits the highest calculation accuracy (the smallest error near 0.0%) comparing to other single or mixed Copula function model.For the computation time (T), it is found that the calculation speed for all Copula function models are faster than MC method, except IPSO-MCopula takes the least computation time.Further, it is seen that the AIC and d 2 obtained by using the IPSO-MCopula function model are much smaller than those obtained by the other six Copula function models, which indicates that our proposed IPSO-MCopula function model has the best goodness of fit among all Copula function models.In short, our proposed IMCMFCM methodology including the novel IPSO-MCopula function model can not only fully describe the correlation among multiple failure modes, but also has high accuracy and fast speed in calculations.

Application of the IPSO-MCopula model to gear reliability analysis
To further verify our proposed IMCMFCM methodology, we now carry out a practical reliability analysis for a gear pair in a mechanical gearbox (shown in Figure 4).We firstly set some proper conditions for the gear pair in consideration.The number of teeth and modulus of the pinion are 19 and 8 mm, respectively, and its material is Cr.The large gear has 38 teeth and    is made from steel.Both gears have a tooth width of 45 i and a dimension accuracy of level 7.The transmitted torque is 1862 N Á m, and the applied circumferential force is 13,316 N. Geometric parameters of the gear pair can be divided to the definite-variable and the random-variable parts.The definite-variable part is listed in Table 7, and the random-variable part is listed in Tables 8 and 9, where Z H is the coefficient of nodal region, Z b is the coefficient of helix angle, Z E is the elasticity coefficient, F t is the tangential force applied on the end face, d 1 is the diameter of indexing circle of pinion, Z e is the overlap coefficient, b is the tooth width, u is the transmission ratio, K Ha is the loaddistribution coefficient for different teeth, K y is the dynamic-load coefficient, K A is the coefficient of working condition, K Hb is the directional loaddistribution coefficient, s H lim is the contact fatigue strength of gear, Z W is the working coefficient of material hardening, Z L is the working coefficient of lubricant condition, Z y is the gear-speed coefficient, Z R is the surface-roughness coefficient of gear, Z X is the dimension coefficient of gear, and Z N is the life coefficient of contact strength of gear.
Based on the concept proposed by Liu and Chen, 33 we establish two marginal distribution functions, one for describing the failure mode of the contact fatigue failure of tooth surface (g 1 (X )) and the other for that of the bending fatigue fracture of tooth root (g 2 (X )), as follows.
(1) Marginal distribution function of the contact fatigue failure of tooth surface: (2) Marginal distribution function of the bending fatigue fracture of tooth root: With using the proposed IMCMFCM methodology, we firstly generate a certain samples by the MC scheme,  and then theoretically build the marginal distributions from the generated data by the kernel density estimation method for the aforementioned two failure modes.Secondly, we calculate the related parameters in the proposed IPSO-MCopula function model by the IPSO and the MLE method with a population size of 80 and a maximum iteration number of 1000.In order to test the influence of different sample sizes on the calculation efficiency and accuracy of the algorithm, the calculation results for different sample sizes are listed in Tables 10 and 11.In addition, the sample of MC method is 10 7 , and the result obtained is used as the truth value to test other algorithms.In Table 10, the failure probability p f i for each failure mode in consideration is calculated by the improved fourth moment method, and the obtained results, including other index parameters.The experimental results show that with the increase of the number of samples, the calculation accuracy of various algorithms will increase, but the calculation efficiency will decrease.Especially, it is seen from Table 11 that the AIC and d 2 obtained by using our proposed IPSO-MCopula function model are smaller than those obtained by the other three Copula functions under the same sample size, which indicates that our proposed IPSO-MCopula model has the best goodness of fit among all.It is worth noting that regardless of the algorithm, the computational efficiency will decrease with the increase of this number.Therefore, when the computational accuracy meets the requirements, there is no need to use a large sample size.
Figure 5(a) (contact fatigue failure mode) and (b) (bending failure mode) show the marginal distribution obtained by IPSO-MCopula and its comparison with empirical values when the sample size is 2000.It is seen that the fitted conditions of the built marginal distribution functions are satisfactory for both failure modes.In conjunction with the contour plots of the probability density function (PDF) and the cumulative distribution function (CDF) for the mixed Copula functions, as shown in Figure 6, it can be found that the PDF and the CDF for the two failure modes of the gear pair are very sensitive to the change of weight coefficient, especially the upper-tail shape in PDF or CDF contour.Among the calculation results of failure probability (p f ) by different Copula function models, our proposed IPSO-MCopula function model exhibits the highest calculation accuracy (the smallest error of 0. 294%) comparing to other Copula functions when the sample size is 2000.And, it (p f =0.1101) is quite consistent with that by the MC scheme (p f =0.1098).In short, even for the practical reliability problems of a gear pair, our proposed IMCMFCM methodology including the novel IPSO-MCopula model performs well for multiple failure modes with a sample size of 2000, and has high accuracy with fast convergence in calculations.In addition, when analyzing smaller failure probability, the MC method requires a larger sample size, and the computational efficiency of the proposed method will be more advantageous.

Conclusion
In this study, a reliability analysis methodology of IMCMFCM for multiple failure modes of gears is proposed based on the novel IPSO-MCopula model which is an optimal combination of the well-known Clayton Copula, Gumbel Copula, and Frank Copula functions in order to more precisely establish the reliability correlation for different failure modes of gears.To optimally determine the calculation parameters in the IPSO-MCopula function model, which are closely related to the computational accuracy and speed, we adopt the MLE method to estimate the dependent parameter of u and the IPSO algorithm to estimate the weight coefficient of l based on the minimum sum of the squares of deviations.To verify the usability of our proposed IMCMFCM methodology, firstly we performed the reliability analysis for multiple failure modes of a standard strut.Secondly, we performed the reliability analysis for multiple failure modes of the gear pair in a practical gearbox system.Finally, We have summarized the following innovative findings: 1.For the standard test function, the proposed IPSO has faster convergence speed and smaller average fitness values compared to SPSO, LDIWPSO and e 2 PSO. 2. The proposed MCopula model includes copulas with different tail relationships, which can be At present, this paper only studies the reliability of the system when the failure mode is explicit, and does not study the failure mode is implicit.This situation can be dealt with by methods such as approximate model, which is also the focus of our future research.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
, and Frank Copula functions to from a new MCopula function.Furthermore, we develop an IPSO algorithm to enhance the calculation performance.The integration of the above two schemes are called the IPSO-MCopula model.Step3 : Calculation of parameters.The dependent parameter u in the IPSO-MCopula model is solved by the MLE scheme.Meanwhile, the IPSO algorithm is used to solve the weight coefficient l in the MCopula model.Step4 : Test of goodness of fit.We test the proposed IPSO-MCopula model with the AIC criterion and the minimum distance criterion.Step5 : Verifications.Use the IPSO-MCopula model to solve the correlations among random variables for different cases.The above manipulation procedure of IMCMFCM Methodology is also shown in Figure1

Figure 3 .
Figure 3. Hollow metal strut subjected to forces at two ends.

Figure 4 .
Figure 4. Configuration of a gear pair used in reliability analysis (cyan and blue parts).

Figure 6 .
Figure 6.Distributions and contour plot of CDF and PDF of the mixed Copula: (a) CDF of the mixed Copula, (b) contour plot of CDF of the mixed Copula, (c) PDF of the mixed Copula, and (d) contour plot of PDF of the mixed Copula.

Table 1 .
Standard test functions.

Table 2 .
Mean and standard deviation by standard test models with different algorithms.

Table 3 .
Average fitness by different models with different population sizes.
N: population size.

Table 6 .
Failure probability and goodness of fit by different Copula function models for a strut system.
E: error.* Comparison basis

Table 7 .
Definite variables of the gear pair geometry.

Table 8 .
Mean and standard deviation of random variables for the tooth contact fatigue failure mode.

Table 9 .
Mean and standard deviation of random variables for the tooth root bending failure mode.

Table 10 .
Failure probability and algorithm parameters by IPSO-MCopula model for a gear pair system.Sample sizes p f1 , p f2 u 1 , u 2 , u 3 l 1 , l 2 , l 3 4. The proposed IPSO-McCopula model hasachieved good results in the analysis of multiple failure modes of gear.As a consequence, the proposed reliability analysis methodology of IMCMFCM is convincible and satisfactory for theoretically and practically investigating the reliability prediction problem of gear systems.