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

Abstract

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.

Keywords

Failure probability of gears reliability analysis of gears Copula function multiple failure modes

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.¹ 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. That means certain degrees of correlations exist among failure modes.^2–7 In the past, many scholars studied the reliability as well as the failure mode problems of mechanic parts.^8–15 Ditlevsen⁸ 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¹⁶ 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.¹⁷ 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.¹⁸ proposed a non-empirical 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.¹⁹ proposed a combined machine learning strategy to evaluate the fatigue of turbine blade disks. Luo et al.²⁰ 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.

Because the Copula function has a concise form and is suitable for being used as a model to describe the non-linear relationship between reliability variables, it has received more attention than the linear correlation coefficient method recently.^21–23 Tang et al.²⁴ used Copula functions to analyze the reliability of a system’s structure. Mi et al.²⁵ 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.²¹ used 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.²⁶ 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,²⁷ and the exponential reduction inertia weight PSO (e²PSO) algorithm based MCopula function.²⁸

Theoretical bases of IMCMFCM methodology

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.²⁹

Supposed 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}, . . ., u_{N})$ that connects the respective marginal distribution functions $F_{1} (x_{1}), F_{2} (x_{2}), . . ., F_{N} (x_{N})$ with the joint distribution function $F (x_{1}, x_{2}, . . ., x_{N})$ , such that for any $x_{1}, x_{2}, . . ., x_{N} \in R$ , the following formula holds:

F (x_{1}, x_{2}, . . ., x_{N}) = C (F_{1} (x_{1}), F_{2} (x_{2}), . . ., F_{N} (x_{N}))

(1)

If $F_{1} (x_{1}), F_{2} (x_{2}), . . ., F_{N} (x_{N})$ are all continuous functions, then C may be uniquely determined. Further, if $F_{1} (x_{1}), F_{2} (x_{2}), . . ., F_{N} (x_{N})$ are all univariate distribution functions, then $F (x_{1}, x_{2}, . . ., x_{N})$ in equation (1) is the joint distribution function of the marginal distribution function $F_{1} (x_{1}), F_{2} (x_{2}), . . ., F_{N} (x_{N})$ . The probability form of the Copula function is expressed as:

\begin{matrix} P (X_{1} \leq x_{1}, X_{2} \leq x_{2}, . . ., X_{n} \leq x_{n}) \\ = F (X_{1} \leq x_{1}, X_{2} \leq x_{2}, . . ., X_{n} \leq x_{n}) \\ = C (u_{1}, u_{2}, . . ., u_{N}) \end{matrix}

(2)

where $u_{i} = F_{i} (x_{i}), (i = 1, 2, \dots, N)$ . Eventually, the failure probability of a system with k failure modes can be expressed as³⁰:

\begin{matrix} P = P (⋃_{j = 1}^{k} f_{j} (x) \leq 0) = \sum_{j = 1}^{k} P_{j} - \sum_{1 \leq j < h \leq k} C (P_{j}, P_{h}) \\ + \sum_{1 \leq j < h < t \leq k} C (P_{j}, P_{h}, P_{t}) - . . . - (- 1)^{k - 1} C (P_{1}, P_{2}, . . ., P_{k}) \end{matrix}

(3)

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.³¹ Different 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.³² Usually, 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:

C_{M} (u_{1}, u_{2}, . . ., u_{N}) = \sum_{k = 1}^{s} λ_{k} C_{k} (u_{1}, u_{2}, . . ., u_{N}; θ_{k})

(4)

where $C_{k} (u_{1}, u_{2}, . . ., u_{N}; θ_{k})$ is the k^th preselected single Copula function, $θ_{k}$ is the linking coefficient, and $λ_{k}$ is the weight coefficient that satisfies $0 \leq λ_{k} \leq 1, \sum_{k = 1}^{s} λ_{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:

C_{M} (u, v) = \sum_{k = 1}^{3} λ_{k} C_{k} (u, v; θ_{k})

(5)

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,34 Secondly, 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 $λ_{k}$ and $θ$ 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 $C (u, v; θ)$ and its density function is $c (u, ν, θ)$ , where

c (u, ν, θ) = \frac{\partial C (u, ν; θ)}{\partial u \partial ν} .

(6)

Then, the joint distribution function of the random variables $X_{1}$ and $X_{2}$ can be expressed as:

F (x_{1}, x_{2}; θ) = C [F_{1} (x_{1}), F_{2} (x_{2}); θ] .

(7)

Further, the joint density function of $F (x_{1}, x_{2}; θ)$ can be obtained as:

f (x_{1}, x_{2}; θ) = \frac{\partial F}{\partial x_{1} \partial x_{2}} = c [F_{1} (x_{1}), F_{2} (x_{2}); θ] f_{1} (x_{1}) f_{2} (x_{2}) .

(8)

The likelihood function of the samples $(X_{1 j}, X_{2 j})$ is :

L (θ) = Π_{j = 1}^{n} f (x_{1 j}, x_{2 j}; θ) = Π_{j = 1}^{n} c [F_{1} (x_{1 j}), F_{2} (x_{2 j}); θ] f_{1} (x_{1 j}) f_{2} (x_{2 j})

(9)

The logarithmic likelihood function are:

\begin{matrix} \ln L (θ) = \sum_{j = 1}^{n} \ln C [F_{1} (x_{1 j}), F_{2} (x_{2 j}); θ] + \sum_{j = 1}^{n} \ln f_{1} (x_{1 j}) \\ + \sum_{j = 1}^{n} \ln f_{2} (x_{2 j}) \end{matrix}

(10)

The optimal value of $θ$ happens at where

θ = \arg max \ln L (θ) .

(11)

The optimal weight coefficient in the MCopula function can be obtained by minimize the Euclidean distance between the empirical Copula function ( $C (\cdot)$ ) and the MCopula function ( $C_{M} (\cdot)$ ) which is defined by

O_{LS} = \sqrt{\frac{1}{n} \sum_{i = 1}^{n} {(C (F_{1} (x_{1 j}), F_{2} (x_{2 j})) - C_{M} (F_{1} (x_{1 j}), F_{2} (x_{2 j})))}^{2}},

(12)

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).

\begin{matrix} \min O_{LS} \\ s . t . {\begin{matrix} \sum_{k = 1}^{3} λ_{k} = 1; \\ 0 \leq λ_{k} \leq 1; k = 1, 2, 3 \end{matrix} \end{matrix}

(13)

where $λ_{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:

\begin{matrix} F_{min}, λ, = F + η \sum_{i = 1}^{n} H (q_{i} (λ)), \\ q_{i} (λ) = max (g_{i} (λ), 0), \\ η = {\begin{matrix} 1; q_{i} (λ) \leq 0.001 \\ 2; \begin{matrix} otherwise \end{matrix} \end{matrix}, \\ H (q_{i} (λ)) = {\begin{matrix} 0; q_{i} (λ) \leq 0.00001 \\ 10; 0.00001 \leq q_{i} (λ) \leq 0.001 \\ 50; \begin{matrix} otherwise \end{matrix} \end{matrix}, \end{matrix}

(14)

where $F$ is the objective function to be optimized, $η$ is the penalty coefficient, $q_{i} (\cdot)$ is the relative constrained penalty function, and $H (\cdot)$ 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 $θ$ and $λ$ 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:

{\begin{matrix} AIC = n \log (\frac{1}{n - 1} \sum_{i = 1}^{n} {(F_{ei} - F_{ci})}^{2}) + 2 K \\ d^{2} = \sum_{i = 1}^{n} {(C (u_{i}, v_{i}) - C_{M} (u_{i}, v_{i}))}^{2} \end{matrix}

(15)

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 (\cdot)$ is the empirical copula function, $C_{M} (\cdot)$ 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,³⁵ 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 $V_{i} = (v_{i 1}, v_{i 2}, . . ., v_{in})$ and $Q_{i} = (q_{i 1}, q_{i 2}, . . ., q_{in})$ , respectively. The particle updates the velocity and position according to the following formula:

\begin{matrix} v_{id} (t + 1) = w v_{id} (t) + c_{1} r_{1} [p_{id} - q_{id} (t)] + c_{2} r_{2} [p_{gd} - q_{id} (t)] \\ q_{id} (t + 1) = q_{id} (t) + v_{id} (t + 1), i = 1, 2, . . ., n \end{matrix}

(16)

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:

w = w_{min} + (w_{max} - w_{min}) \cdot e^{- {[t / (\frac{T}{4})]}^{2}}

(17)

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.

{\begin{matrix} c_{1} = c_{1, max} + \frac{c_{1, min} - c_{1, max}}{T} t \\ c_{2} = c_{2, min} + \frac{c_{2, max} - c_{2, i min}}{T} t \end{matrix}

(18)

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.

Step2 : MCopula model. We mix the Gumbel, Clayton, 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 $θ$ in the IPSO-MCopula model is solved by the MLE scheme. Meanwhile, the IPSO algorithm is used to solve the weight coefficient $λ$ 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 Figure 1

Figure 1.

Procedure of IMCMFCM methodology.

Results and discussion

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²⁷ and e²PSO²⁸ 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.³⁶ In 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).

Table 1.

Standard test functions.³⁶

Model	Mathematical expression	Range of search
Sphere	$f_{1} = \sum_{i = 1}^{n} {x_{i}}^{2}$	${(- 100, 100)}^{n}$
Rosenbrock	$f_{2} = \sum_{i = 1}^{n} [100 {(x_{i} + 1 - {x_{i}}^{2})}^{2} + {(x_{i} - 1)}^{2}]$	${(- 30, 30)}^{n}$
Griewank	$f_{3} = \frac{1}{4000} \sum_{i = 1}^{n} {x_{i}}^{2} - Π_{i = 1}^{n} \cos (\frac{x_{i}}{\sqrt{i}}) + 1$	${(- 600, 600)}^{n}$
Ackley	$f_{4} = - 20 \exp (- 0.2 \sqrt{\frac{1}{n} \sum_{i = 1}^{n} {x_{i}}^{2}} - \exp (\frac{1}{n} \sum_{i = 1}^{n} \cos (2 π x_{i}))) + 20 + e$	${(- 32, 32)}^{n}$

Figure 2.

Comparison of the optimization effect using different calculation algorithms: (a) sphere model, (b) rosenbrock model, (c) griewank model, and (d) ackley model.

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.

Table 2.

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

Model	Indicator	SPSO	LDIPSO	e ²PSO	IPSO
Sphere	Mean	0.528	0.071	0.271	0
Sphere	Standard deviation	0.2	0.126	0.228	0
Rosenbrock	Mean	82.032	36.258	58.59	28.721
Rosenbrock	Standard deviation	19.8	13.872	21.77	0.068
Griewank	Mean	0.025	0.008	0.018	0
Griewank	Standard deviation	0.011	0.007	0.010	0
Ackley	Mean	1.402	0.239	0.887	8.88 × 10⁻¹⁶
Ackley	Standard deviation	0.340	0.407	0.592	9.86 × 10⁻³²

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.

Table 3.

Average fitness by different models with different population sizes.

Model	$N = 40$	$N = 60$	$N = 80$	$N = 100$	$N = 200$
Sphere	0	0	0	0	0
Rosenbrock	28.721	28.66	28.622	28.482	27.861
Griewank	0	0	0	0	0
Ackley	8.88 × 10⁻¹⁶	8.88 × 10⁻¹⁶	8.88 × 10⁻¹⁶	8.88 × 10⁻¹⁶	8.88 × 10⁻¹⁶

N: population size.

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³ (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:

g_{1} = σ_{s} - \frac{P}{π dc},

(19)

and

g_{2} = \frac{π^{2} E [{(d + 2 c)}^{4} - d^{4}]}{64 l^{2} dt} .

(20)

Figure 3.

Hollow metal strut subjected to forces at two ends.

In equations (19) and (20), $σ_{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.

Table 4.

Parameters in equations (19) and (20).³

Variable	Mean	Standard deviation	Distribution type
$P (kN)$	120	12	Normal
$σ_{s} (GPa)$	490	4	Normal
$d (mm)$	42.5	0.212	Normal
$c (mm)$	2.4	0.012	Normal
$E (MPa)$	203	4	Normal
$l (mm)$	1090	5.45	Normal

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 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⁸, and the result obtained is used as the truth value to test other algorithms.

Table 5.

Failure probability and algorithm parameters by IPSO-MCopula model for a strut system.

$p_{f_{1}}, p_{f_{2}}$	$θ_{1}, θ_{2}, θ_{3}$	$λ_{1}, λ_{2}, λ_{3}$
0.01766, 0.02694	1.27, 1.73, 5.10	0.33, 0.21, 0.45

Table 6.

Failure probability and goodness of fit by different Copula function models for a strut system.

parameters	Method
	SPSO-MCopula	LDIPSO-MCopula	e ²PSO-MCopula	IPSO-MCopula	MC
$p_{f}$	0.038	0.038	0.039	0.039	0.039
$d^{2}$	2.6	2.6	0.0	0.0	*
AIC	0.256	0.272	0.256	0.255	*
Error (%)	−4171.7	−4117.4	−4169.4	−4175.0	*
T (s)	1.6855	1.5599	1.5905	1.3698	11.109

E: error.

Comparison basis

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 \cdot 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_{β}$ 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_{ε}$ is the overlap coefficient, b is the tooth width, u is the transmission ratio, $K_{H α}$ is the load-distribution coefficient for different teeth, $K_{υ}$ is the dynamic-load coefficient, $K_{A}$ is the coefficient of working condition, $K_{H β}$ is the directional load-distribution coefficient, $σ_{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_{υ}$ 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.

Figure 4.

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

Table 7.

Definite variables of the gear pair geometry.

Variable	$Z_{H}$	$Z_{ε}$	$Z_{β}$	$m_{n}$	$d_{1}$	$b$	$Z_{X}$	$Y_{ε}$	u	$Y_{β}$	$Y_{X}$
Values	2.5	0.875	1	8 mm	152 mm	45 mm	1	0.675	2	1	1

Table 8.

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

Variable	$σ_{H lim}$	$Z_{N}$	$Z_{R}$	$Z_{υ}$	$Z_{W}$	$Z_{L}$	$Z_{E}$	$F_{t}$	$K_{A}$	$K_{υ}$	$K_{H α}$	$K_{H β}$
Mean	700	1.5	1.02	1.05	1	0.95	189.8	13 316	1.25	1.08	1	1.22
Standard deviation	62	0.043	0.021	0.021	0.03	0.019	4.745	210.65	0.063	0.12	0.03	0.061

Table 9.

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

Variable	${σ_{F}}_{lim}$	$Y_{ST}$	$Y_{NT}$	$Y_{δ relT}$	$Y_{RelT}$	$F_{t}$	$Y_{Fa}$	$Y_{sa}$	$K_{A}$	$K_{υ}$	$K_{F α}$	$K_{F β}$
Mean	300	1.5	0.85	0.99	1.02	13 316	2.8	1.51	1.25	1.08	1	1.18
Standard deviation	60	0.05	0.032	0.028	0.03	210.65	0.089	0.062	0.063	0.12	0.03	0.059

Based on the concept proposed by Liu and Chen,³³ 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:

\begin{matrix} Z_{1} = g_{1} (X) = σ_{H lim C} - σ_{H} = σ_{H lim} Z_{N} Z_{R} Z_{υ} Z_{W} Z_{L} Z_{X} \\ - Z_{H} Z_{E} Z_{ε} Z_{β} \sqrt{\frac{F_{t}}{d_{1} b} \frac{u \pm 1}{u} K_{A} K_{υ} K_{H α} K_{H β}} \end{matrix}

(21)

(2) Marginal distribution function of the bending fatigue fracture of tooth root:

\begin{matrix} Z_{2} = g_{2} (X) = σ_{F lim C} - σ_{F} \\ = {σ_{F}}_{lim} Y_{ST} Y_{NT} Y_{δ relT} Y_{RelT} Y_{X} - \frac{F_{t}}{b m_{n}} Y_{Fa} Y_{sa} Y_{ε} Y_{β} K_{A} K_{υ} K_{F α} K_{F β} \end{matrix}

(22)

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⁷, 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.

Table 10.

Failure probability and algorithm parameters by IPSO-MCopula model for a gear pair system.

Sample sizes	$p_{f_{1}}, p_{f_{2}}$	$θ_{1}, θ_{2}, θ_{3}$	$λ_{1}, λ_{2}, λ_{3}$
500	0.10698, 0.00402	0.01476, 1.02231, 0.20463	0.42475, 0.52663, 0.06324
1000	0.10698, 0.00402	0.14476, 1.08730, 0.75403	0.20975, 0.70868, 0.08157
1500	0.10698, 0.00402	0.15185, 1.05790, 0.72730	0.32794, 0.30030, 0.37176
2000	0.10698, 0.00402	0.20793, 1.10352, 1.06825	0.17396, 0.53585, 0.29019
2500	0.10698, 0.00402	0.13176, 1.06691, 0.71566	0.46091, 0.51047, 0.02862

Table 11.

Failure probability and goodness of fit by different Copula function models for a gear pair system.

Sample sizes	Parameters	Method
Sample sizes	Parameters	SPSO-MCopula	LDIPSO-MCopula	e²PSO-MCopula	IPSO-MCopula	MC
500	$p_{f}$	0.1105	0.1105	0.1105	0.1105	0.1098
	$d^{2}$	0.2188	0.2204	0.2164	0.1805	—
	AIC	−1996.96	−1993.75	−2001.77	−2080.56	—
	Error (%)	0.6346	0.6336	0.6308	0.6124	—
	T (s)	0.148	0.157	0.156	0.144	2.761
1000	$p_{f}$	0.1104	0.1104	0.1104	0.1102	0.1098
	$d^{2}$	0.2637	0.2456	0.2689	0.2379	—
	AIC	−4145.46	−4207.20	−4128.47	−4234.62	—
	Error (%)	0.5253	0.4746	0.5277	0.3666	—
	T (s)	0.223	0.292	0.283	0.221	2.761
1500	$p_{f}$	0.1103	0.1102	0.1104	0.1102	0.1098
	$d^{2}$	0.2268	0.2246	0.2262	0.2219	—
	AIC	−6684.71	−6697.40	−6688.01	−6713.52	—
	Error (%)	0.3947	0.3523	0.5015	0.3018	—
	T (s)	0.397	0.397	0.394	0.325	65.995
2000	$p_{f}$	0.1102	0.1102	0.1103	0.1101	0.1098
	$d^{2}$	0.2263	0.2250	0.2161	0.2101	—
	AIC	−9171.11	−9181.20	−9251.10	−9292.55	—
	Error (%)	0.3290	0.3222	0.4036	0.2945	—
	T (s)	0.480	0.486	0.494	0.450	2.761
2500	$p_{f}$	0.1102	0.1102	0.1102	0.1101	0.1098
	$d^{2}$	0.2975	0.3009	0.2984	0.2945	—
	AIC	−11,114.91	−11,090.46	−11,108.23	−11,137.05	—
	Error (%)	0.2895	0.2914	0.2863	0.2567	—
	T (s)	0.726	0.793	0.788	0.612	2.761

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.

Figure 5.

Comparisons between the theoretical and the empirical marginal distributions: (a) contact-fatigue failure mode of tooth surface and (b) bending failure mode of tooth root.

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.

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 $θ$ and the IPSO algorithm to estimate the weight coefficient of $λ$ 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:

For the standard test function, the proposed IPSO has faster convergence speed and smaller average fitness values compared to SPSO, LDIWPSO and e²PSO.

The proposed MCopula model includes copulas with different tail relationships, which can be well used to describe the correlation between multiple failure modes.

The proposed IPSO-MCopula model, with only using a small amount of sample data, has higher computational efficiency and reliability-prediction ability than that obtained by the traditional MC method with using a large number of samples.

The proposed IPSO-McCopula model has achieved 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.

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.

Footnotes

Handling Editor: Chenhui Liang

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.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Natural Science Foundation of Fujian Province (No. 2022J011181), the Special Project of Central Government Guiding Local Science and Technology Development (No. 2021L3029), the Leading Technology Project of Sanming City (No. 2023-G-2 and No. 2022-G-13), the Major science and technology projects of Fujian Province (No. 2022HZ026025), the Department of Education Science and Technology Program of Fujian Province (No. JAT200641), and the National Natural Foundation Cultivation Project of Sanming University (No. PYT2006). Also, this work was supported by the Project of the Department of Science and Technology of Fujian Province, China (2021G02013, 2020H0049, 2021H0060), and in part by the Sanming University of Fujian Province, China (19YG05, 19YG04).

ORCID iD

Wang Kun-Chieh

References

Yan

Shen

Sun

, et al. Knowledge transfer for rotary machine fault diagnosis. IEEE Sens J 2019; 20: 8374–8393.

Zhang

Zhao

A stress-strength time-varying correlation interference model for structural reliability analysis using copulas. IEEE Trans Reliab 2017; 66: 351–365.

Wang

Zhao

Xia

, et al. Analyzing reliability and sensitivity to reliability of rescue attachment with multiple failure modes. J Braz Soc Mech Sci Eng 2021; 43: 1–16.

Zhang

Wilson

System reliability and component importance under dependence: a copula approach. Technometrics 2017; 59: 215–224.

Bichon

McFarland

Mahadevan

Efficient surrogate models for reliability analysis of systems with multiple failure modes. Reliab Eng Syst Saf 2011; 96: 1386–1395.

Xie

Qian

NX.

Reliability model for complex mechanical component design. Appl Mech Mater 2013; 365: 28–31.

XZ.

Assessing the correlated performance functions of an engineering system via probabilistic analysis. Struct Saf 2015; 52: 10–19.

Ditlevsen

Narrow reliability bounds for structural systems. J Struct Mech 1979; 7: 453–472.

Valdebenito

Schuëller

GI.

A survey on approaches for reliability-based optimization. Struct Multidiscipl Optim 2010; 42: 645–663.

10.

Ditlevsen

Bjerager

Methods of structural systems reliability. Struct Saf 1986; 3: 195–229.

11.

Ditlevsen

Model uncertainty in structural reliability. Struct Saf 1982; 1: 73–86.

12.

Darestani

Padgett

Shafieezadeh

Parametrized wind-surge-wave fragility functions for wood utility poles. J Struct Eng 2022; 148: 04022057.

13.

Song

Bai

GC.

Deep learning regression-based stratified probabilistic combined cycle fatigue damage evaluation for turbine bladed disks. Int J Fatigue 2022; 159: 106812.

14.

Song

Choy

, et al. Multivariate ensembles-based hierarchical linkage strategy for system reliability evaluation of aeroengine cooling blades. Aerosp Sci Technol 2023; 138: 108325.

15.

Song

Bai

GC.

Recent advances in reliability analysis of aeroengine rotor system: a review. Int J Struct Integr 2022; 13: 1–29.

16.

Cui

Blockley

DI.

On the bounds for structural system reliability. Struct Saf 1991; 9: 247–259.

17.

Sun

Zhang

Reliability computing new theory of mechanical elements with multi failure mode. Chin J Mech Eng 2003; 39: 134–138.

18.

Zhou

Xie

Han

, et al. Method for syetem risk prediction with load dependency based on Nataf transformation. J Mech Eng 2009; 45: 137–141.

19.

Zhu

Niu

Keshtegar

, et al. Machine learning-based probabilistic fatigue assessment of turbine bladed disks under multisource uncertainties. Int J Struct Integr 2023; 14: 1000–1024.

20.

Luo

Zhu

Keshtegar

, et al. An enhanced uniform simulation approach coupled with SVR for efficient structural reliability analysis. Reliab Eng Syst Saf 2023; 237: 109377.

21.

Sun

Wang

, et al. Remaining useful life prediction for a machine with multiple dependent features based on Bayesian dynamic linear model and copulas. IEEE Access 2017; 5: 16277–16287.

22.

Bin

Shahzad

Bing

, et al. Probabilistic computational model for correlated wind farms using copula theory. IEEE Access 2018; 6: 14179–14187.

23.

Jin

Mesiar

Cognitive integrals with its generalized and adapted forms. IEEE Trans Fuzzy Syst 2018; 26: 1960–1969.

24.

Tang

Zhou

, et al. Impact of copulas for modeling bivariate distributions on system reliability. Struct Saf 2013; 44: 80–90.

25.

Zhang

, et al. A copula-function-based structural system reliability analysis method. Acta Mech 2022; 233: 1371–1391.

26.

Jiang

Zhang

Han

, et al. A vine-copula-based reliability analysis method for structures with multidimensional correlation. J Mech Des 2015; 137: 061405.

27.

Shi

Eberhart

. Parameter selection in particle swarm optimization. In: Proceedings of the 7th international conference on evolutionary programming, San Diego, CA, 1998, VII, pp. 591–600. Springer Berlin Heidelberg.

28.

Chen

Huang

Jia

, et al. Natural exponential inertia weight strategy in particle swarm optimization. In: 6th World Congress on Intelligent Control and Automation, Dalian, 2006, pp. 3672–3675. New York: IEEE.

29.

Nelsen

RB.

An introduction to copulas. 2nd ed. New York, NY: Springer series in statistics, 2006.

30.

Sun

Luo

Zhang

, et al. Reliability analysis of stochastic structure with multi-failure modes based on mixed copula. Eng Fail Anal 2019; 105: 930–944.

31.

Jiang

Zhang

Wang

, et al. Structural reliability analysis using a copula-function-based evidence theory model. Comput Struct 2014; 143: 19–31.

32.

Gunawan

Khaled

Kohn

Mixed marginal copula modeling. J Bus Econ Stat 2020; 38: 137–147.

33.

Liu

Chen

Dynamic reliability evaluation of high-speed train gearbox based on copula function. IEEE Access 2022; 10, 51792–51803.

34.

Wang

Pham

Modeling the dependent competing risks with multiple degradation processes and random shock using time-varying copulas. IEEE Trans Reliab 2011; 61: 13–22.

35.

Kennedy

Eberhart

Particle swarm optimization. In Proceedings of ICNN’95-international conference on neural networks, Perth, WA, Australia, 27 November–1 December 1995, vol. 4, pp. 1942–1948. New York: IEEE.

36.

Yaghoobi

Esmaeili

An improved artificial bee colony algorithm for global numerical optimization. Int J Bio-Inspir Com 2017; 9: 251–258.