# Research on Reliability Modelling for Heavy CNC Machine Tools under Uncertain Variables Based on Universal Generating Function

Authors: Honghua Sun; Lipeng Xu; Chunwei Li; Hongxia Chen; Guangjun Jiang; Xiaobo Nie
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Due to the lack of historical data, the inability to carry out a large number of reliability tests on the machine tool, the lack of reliability information collected, and the change of the working environment of the machine tool, there are a lot of uncertainties in the machine tool. In the existing research, only unilateral uncertainty factors are usually considered, while in practical problems, random uncertainty and cognitive uncertainty exist at the same time. Therefore, in this study, the random uncertainty and cognitive uncertainty in the system are characterized by random variables, fuzzy variables, probability box variables, and interval variables. For the structural function with more variables, the dimension is reduced first, and then the reliability model of the universal generating function (UGF) is established. Taking the heavy CNC machine tool as an example, the structural reliability analysis model of the fatigue strength of the milling shaft based on the UGF is constructed. The local sensitivity analysis and global sensitivity analysis of the variables affecting the fatigue strength of the milling shaft are carried out, and the factors that have the greatest impact on the fatigue strength of the milling shaft are obtained. The case study shows the effectiveness of the method proposed in this study.


---

## Body

## 1. Introduction

CNC machine tools are widely used in modern manufacturing industry, known as the “machine tool” of modern manufacturing industry, and its performance level also directly affects the rapid development of manufacturing industry in the whole country [1]. Heavy CNC machine tool has the characteristics of complex system structure, difficult process, long R&D cycle, high manufacturing cost, large driving load, and complex working condition stress. Due to the high R&D cost, long development cycle, insufficient experimental conditions, small batch customization, and other factors of Heavy CNC machine tools, it is impossible to obtain sufficient test data and field use data, which makes the reliability data include random uncertainty and cognitive uncertainty. How to deal with uncertainty is the key and difficulty of reliability evaluation. Uncertainty includes random uncertainty and cognitive uncertainty. Random uncertainty comes from the randomness and volatility of the system itself, which can be quantified by probability theory; cognitive uncertainty is due to cognitive bias, incomplete information, and other reasons. The subjective information given by experts or engineers often has cognitive uncertainty, which can be gradually reduced with the deepening of understanding and the increase of information [2]. At present, the main theories dealing with cognitive uncertainty include fuzzy mathematics, possibility theory, interval analysis, evidence theory, and Bayesian method. These theories still have certain limitations in dealing with cognitive uncertainty, and these theories need to make various assumptions about the probability distribution of component parameters or the independence of components [3], the reliability evaluation method based on various assumptions has a certain degree of unreliability. In addition, none of the above theories can effectively quantify random uncertainty and cognitive uncertainty at the same time. How to establish a unified model to quantify random uncertainty and cognitive uncertainty is the main problem of complex system reliability evaluation.Xiahou et al. [4] considered the uncertainty of component degradation parameters in the system, quantitatively characterized the component degradation parameters by evidence theory, and established the system reliability evaluation model based on evidence theory. Yang [5] combined evidence theory with Bayesian network to expand the traditional fault tree with cognitive uncertainty. Mi [6] combined evidence theory with common cause failure and proposed a reliability analysis method of complex systems under the coexistence of cognitive uncertainty and common cause failure. Lü [7] proposed a random fuzzy uncertainty reliability analysis method. The fuzzy variable is transformed into an interval variable by using the cut set method, and the Taylor series is expanded to calculate the upper and lower bounds of reliability index and failure probability. Mi [8] uses the probability boxes to quantify uncertain parameters. Zhang [9] quantified the uncertainty of parameters through fuzzy numbers and probability boxes and proposed a method based on entropy equivalent transformation and saddle point approximation for reliability modeling and analysis. Wang [10] transformed fuzzy variables into interval variables and established a reliability model containing random variables and fuzzy variables. Tang [11] proposed a mixed reliability discretization analysis method of uncertain structures based on evidence theory to solve the coexistence of random, fuzzy, and other uncertainties in structural reliability analysis. It can be seen that many methods have been applied to uncertainty quantification. These research results provide important technical means for problems such as the lack of complex reliability data.As a bridge between continuous mathematics and discrete mathematics, UGF is a common method in reliability analysis. This method was first proposed by Ushakov [12] and then extended by Levitin [13]. Mi [14] combined evidence theory with UGF and proposed a reliability model of a complex system based on trust UGF. Li [15] combined interval analysis with UGF and proposed an interval UGF model for system reliability analysis. Xiao [16] established a UGF model under a variety of mixed variables. Ding [17] proposed fuzzy UGF. Shen [18] used Markov model and UGF to analyze the reliability of bearings and solve the problem of state space explosion. Ma [19] combined the common cause failure problem with the fuzzy UGF on the basis of the fuzzy UGF. Yuan [20] proposed a reliability evaluation model of a complex system based on trust UGF and quantified the cognitive uncertainty in the system by using evidence theory.At present, the related research on uncertainty quantification has achieved some results. However, there is no way to deal with both cognitive and random uncertainty. There are few UGF models under imprecise probability [18–20]. Based on the existing, this project studies the UGF model of uncertain variables, studies the random and cognitive uncertainty quantification based on the UGF, and establishes the reliability evaluation model. The feasibility and reliability of the two methods are fully verified by Monte Carlo analysis under the condition of small uncertainty.


## 2. Mathematical Description of UGF

The UGF represents discrete variables in the form of polynomials, defines the calculation rules and combination operators of discrete variables, and obtains the general form of system construction through recursive calculation.Let all the possibilities of the discrete variableX be x1,x2,…,xk, and the corresponding probabilities are p1,p2,…pk respectively. Then the mapping of xi⟶pi is called the probability quality function and satisfies ∑i=1kpi=1. The UGF of a discrete variable X is defined as(1)mt=Eetx=∑i=1ketxipi.LetZ represent et, then the UGF of X is expressed as(2)uz=∑i=1kpi⋅zxi.For two universal generating functions, the algorithm of [21] u1z=∑i=1kpi⋅zxi., u2z=∑i=1kpi⋅zxi:(3)u1z+u2z=∑i=1k1∑j=1k2pipjzxi+xj,u1z−u2z=∑i=1k1∑j=1k2pipjzxi−xj,u1z×u2z=∑i=1k1∑j=1k2pipjzxi×xj,u1z÷u2z=∑i=1k1∑j=1k2pipjzxi÷xj.According to the algorithm of the UGF, considering “n” independent discrete variables X=X1,X2,…,Xn, the UGF ugz of the function M=gX can be expressed as(4)ugz=uX1z⊗uX2z⊗,…,⊗uXnz=∑j1=1k1∑j2=1k2⋯∑jnkn∏i=1npijizgXij1,…,Xijn,where, ⊗ is a compound operator, gXij1,…,Xijn is calculated according to the four algorithms of M=gX. It can be seen from equation (4) that the number of items in the equation is ∏i=1nki. For simple calculation, equation (4) can be rewritten as(5)uz=∑i=1kpizμi,where K=∏i=1nki.


### 2.1. UGF for Different Types of Variables

#### 2.1.1. UGF Representation of Random Variables

Let the probability density function of a random variableX be fx, where the distribution of samples is known and bounded, and it is recorded as a0,an. Divide a0,an , into n subintervals, as shown in Figure 1. The corresponding mean and probability of each subinterval are, respectively:(6)μi=a0+a12,…,an−1+an2,Pi=∫ai−1aifxdx,i=1,2,…,n.Figure 1
Probability distribution function of random variables.The universal generating function of random variables from equation (2) can be expressed as(7)uz=∑i=1npizμi.


#### 2.1.2. UGF Representation of Probability Box Variables

Let the probability density function of the probability box variableX⌢ be fx⌢ and the distribution interval be a¯,a¯. To determine the mean value and probability of each interval of X, the interval a¯,a¯ can be divided into l cells on average, which can be expressed as(8)a1¯,a1¯,a2¯,a2¯,…,al¯,al¯.The mean value of each interval is(9)μi=μ1,μ2,…,μl=a1¯+a1¯2,a2¯+a2¯2,…,al¯+al¯2.The probability value corresponding to the mean value of each interval is(10)pi=p1,p2,…,pl.And, the probability value meets the following conditions:(11)p1+p2+⋯+pi=1,piL≤pi≤piUi=1,2,…,l.piU=min∫ai¯ai¯fx⌢dx⌢i=1,2,…,l, piU=max∫ai¯ai¯fx⌢dx⌢i=1,2,…,l thus, the UGF of the probability box variable X⌢ is expressed as follows:(12)μz=∑i=1lpizμi.


#### 2.1.3. UGF Representation of Interval Variables

For a variableY, if the upper and lower bounds of the variable are known, it is recorded as Y=b¯,b¯ and the interval Y is divided into r cells, which are recorded as b¯i,b¯i,i=1,2,…,r at this time, the mean value and probability of each interval of interval variable Y can be expressed as(13)bi=b1,b2,…,br=b1¯+b1¯2,b2¯+b2¯2,…,br¯+br¯2,pi=p1,p2,…,pr,Therefore, the UGF of interval variables is expressed as(14)μz=∑i=1rpizbi.


#### 2.1.4. UGF Representation of Fuzzy Variables

Fuzzy information is also a common type of information in reliability engineering. Fuzzy information is described by the membership function. Triangular fuzzy number is a common membership function, which is represented byX˜ represents the fuzzy variable, and the triangular fuzzy membership function is shown in Figure 2:Figure 2
Membership function of triangular fuzzy number.(15)μx˜=x˜−ab−a,a≤x˜<b,1,x˜=b,x˜−cb−c,b<x˜≤c,0,others.For the fuzziness problems contained in variables, the existing theoretical methods often use the form of cut sets to transform the membership function into the form of intervals. In this paper, the method of entropy equivalence will be used to transform the fuzziness problems in variables into random problems.Shannon defines the probability entropy of random variables under uncertainty as:(16)PEx=−∫−∞+∞fxlnfxdx.For random variables with normal distribution, if the meanμx and variance σx is known, equation (16) can be transformed into(17)PEx=12+ln2πσx.Fuzzy variableX˜ the probability entropy of a is defined as(18)PEx˜=−∫acμ′x˜lnμ′x˜dx˜,where, μ′x˜=μx˜/∫acμx˜dx˜.To obtain the probability density function of equivalent normal random variables, it is necessary to determine the fuzzy variableX˜. According to reference [22], the mean and standard deviation of a are equivalent to the normal random variable x and the fuzzy variable X˜ have the same probability entropy, PEx=PEx˜Obtained by combining equation (17):(19)σx˜=12πePEx˜−0.5.Becauseμ′x˜ satisfies the property of probability density, it can be called a probability density function:(20)μx˜=∫acx˜μ′x˜dx=∫acx˜μx˜∫acμx˜dxdx˜=1∫acμx˜dx˜∫acx˜μx˜dx˜.The probability density function of the equivalent normal random variable of the fuzzy variable is expressed as(21)fx=12πσx˜e−x−μx˜/2σx˜2.The fuzzy random variable is divided intoK intervals, and the UGF is expressed as(22)μz=∑i=1kpizμi.


### 2.2. Structural Reliability Modeling with Mixed Uncertainty

Variables describing mixed uncertainty usually include: random variableX, fuzzy variable X˜, interval variable Y, probability box variable X⌢, if (X,X˜, Y,X⌢) is a variable that affects the state and nature of the structure, then the functional function of the structure can be expressed as(23)Z=gX,X˜,Y,X⌢.According to the stress strength interference model, whenZ<0 is, the structure fails; when Z>0 the structure is in a safe state; when Z=0, the structure reaches the limit state. In Figure 3, three limit states of the structural function in the case of two-dimensional variables.Figure 3
Structural limit state of two-dimensional variables.In general, the failure rate of the structure is obtained by multiple integration of the joint probability density function of each variable in the region ofZ<0. When the structure contains different types of variables and the number of variables is large, it is very difficult to solve it by multiple integrals. The UGF simplifies the calculation to a certain extent. The structural failure probability under various uncertainties based on the UGF can be expressed as(24)F=PfgX,X˜,Y,X⌢<0.When the number of variables affecting reliability is small, the calculation using equation (24) is relatively simple, but when the number of variables affecting reliability is large, it will lead to the problem of large amount of calculation. Therefore, the approximate dimension reduction method is used to reduce the dimension of the multivariable computing function, obtain the approximate function of the multivariable computing function, and then build a UGF model to solve the failure efficiency.


#### 2.2.1. Dimensionality Reduction of Multivariable Function

Rahman and Xu [23] proposed an approximate dimensionality reduction method for multivariable functions. The main idea of the multivariable dimensionality reduction method is to decompose the original multivariable function into the sum of multiple single variables. Multivariable dimensionality reduction methods are mainly divided into two categories: Multivariable dimensionality reduction method based on mean point expansion and multivariable dimensionality reduction method based on the maximum possible failure point. This paper mainly adopts the multivariable dimensionality reduction method based on mean point expansion. Let the mean value of random variables in n intervals μ=μx1,μx2,…,μxn, and then the dimension reduction approximation model of the function gX of the mean value is as follows:(25)gX≅g^X=∑i=1ngμx1,μx2,…,Xi,μxn−n−1gμx1,μx2,…,μxn.For variables with different types, the mean values of (X,X˜, Y, X⌢) are, respectively, μx=μx1,μx2,…,μxn, μx˜=μx˜1,μx˜2,…,μx˜l, μY=μY1,μY2,…,μYr, μx⌢=μx⌢1,μx⌢2,…,μx⌢k it can be seen from equation (25) that the dimension reduction approximation of the functional function under a variety of mixed variables can be expressed as(26)gX,X˜,Y,X⌢≅g^X,X˜,Y,X⌢=g1X+g2X˜+g3Y+g4X⌢.If the fuzzy variable is transformed into a normal random variable, based on the equivalent method of constant entropy, equation (26) can be rewritten as(27)gX,X˜,Y,X⌢≅g^X,X˜,Y,X⌢=g⇀X+g3Y+g4X⌢,where g⇀X=g1X+g2X˜−n+l+r+k−1gμX,μX˜,μY,μX⌢.


#### 2.2.2. Estimate the Probability Density Function ofZ=g⇀X Based on the Maximum Entropy Estimation Method

Z=g⇀X, X is a random variable, so Z is also a random variable. The maximum entropy estimation method does not need any artificial assumption when approaching the unknown probability density function. It is an effective and high-precision method. In 1957, jayncs proposed the calculation method of maximum entropy. The model of estimating Z=g⇀X with maximum entropy is as follows:(28)maxE=max∫ab−fzlnpzdz.The following constraints are met:(29)∫abfzdz=1,(30)mi=∫abzifzdz,i=1,2,…,s.The value ofmi can be solved by simulation and its model is(31)mi=1N∑j=1Ng⇀Xji.fz is the probability density function and mi is the i-th order origin moment of fz. According to the Lagrange method, the probability density function fz of the maximum entropy distribution is expressed as(32)fz=expλ0+∑i=15λizi,λ0=ln∫abexp∑i=1sλizidz.λ0, λi, i=1,2,…,s is a Lagrange multiplier which can be obtained by the nonlinear programming method. The probability density function expression of multiple variables with the maximum entropy distribution satisfying the constraint conditions in equation (31):(33)fz=expλ0+λ1z+λ2z2+λ3z3+⋯+λszs.Then it is transformed into the form of a universal generating function according to equation (2).Since the functiong3Y is an interval variable, the minimum and maximum values are determined through the linear programming optimization model, that is,(34)g3Lg3U=minmaxgY,s.t.Y1L≤Y1≤Y1U,⋮⋮⋮YrL≤Yr≤YrU.According to equations (26) and (34), the maximum and minimum values of the functional functions gX,X˜,Y,X⌢ in the interval variables are, respectively.(35)gX,X˜,Y,X⌢=gmin=g⇀X+g4X⌢+g3L,gX,X˜,Y,X⌢=gmax=g⇀X+g4X⌢+g3U.The minimum and maximum failure probability values can be obtained by using the calculation method of UGF as shown in Figure4.Figure 4
Flow chart of failure probability procedure.


## 2.1. UGF for Different Types of Variables

### 2.1.1. UGF Representation of Random Variables

Let the probability density function of a random variableX be fx, where the distribution of samples is known and bounded, and it is recorded as a0,an. Divide a0,an , into n subintervals, as shown in Figure 1. The corresponding mean and probability of each subinterval are, respectively:(6)μi=a0+a12,…,an−1+an2,Pi=∫ai−1aifxdx,i=1,2,…,n.Figure 1
Probability distribution function of random variables.The universal generating function of random variables from equation (2) can be expressed as(7)uz=∑i=1npizμi.


### 2.1.2. UGF Representation of Probability Box Variables

Let the probability density function of the probability box variableX⌢ be fx⌢ and the distribution interval be a¯,a¯. To determine the mean value and probability of each interval of X, the interval a¯,a¯ can be divided into l cells on average, which can be expressed as(8)a1¯,a1¯,a2¯,a2¯,…,al¯,al¯.The mean value of each interval is(9)μi=μ1,μ2,…,μl=a1¯+a1¯2,a2¯+a2¯2,…,al¯+al¯2.The probability value corresponding to the mean value of each interval is(10)pi=p1,p2,…,pl.And, the probability value meets the following conditions:(11)p1+p2+⋯+pi=1,piL≤pi≤piUi=1,2,…,l.piU=min∫ai¯ai¯fx⌢dx⌢i=1,2,…,l, piU=max∫ai¯ai¯fx⌢dx⌢i=1,2,…,l thus, the UGF of the probability box variable X⌢ is expressed as follows:(12)μz=∑i=1lpizμi.


### 2.1.3. UGF Representation of Interval Variables

For a variableY, if the upper and lower bounds of the variable are known, it is recorded as Y=b¯,b¯ and the interval Y is divided into r cells, which are recorded as b¯i,b¯i,i=1,2,…,r at this time, the mean value and probability of each interval of interval variable Y can be expressed as(13)bi=b1,b2,…,br=b1¯+b1¯2,b2¯+b2¯2,…,br¯+br¯2,pi=p1,p2,…,pr,Therefore, the UGF of interval variables is expressed as(14)μz=∑i=1rpizbi.


### 2.1.4. UGF Representation of Fuzzy Variables

Fuzzy information is also a common type of information in reliability engineering. Fuzzy information is described by the membership function. Triangular fuzzy number is a common membership function, which is represented byX˜ represents the fuzzy variable, and the triangular fuzzy membership function is shown in Figure 2:Figure 2
Membership function of triangular fuzzy number.(15)μx˜=x˜−ab−a,a≤x˜<b,1,x˜=b,x˜−cb−c,b<x˜≤c,0,others.For the fuzziness problems contained in variables, the existing theoretical methods often use the form of cut sets to transform the membership function into the form of intervals. In this paper, the method of entropy equivalence will be used to transform the fuzziness problems in variables into random problems.Shannon defines the probability entropy of random variables under uncertainty as:(16)PEx=−∫−∞+∞fxlnfxdx.For random variables with normal distribution, if the meanμx and variance σx is known, equation (16) can be transformed into(17)PEx=12+ln2πσx.Fuzzy variableX˜ the probability entropy of a is defined as(18)PEx˜=−∫acμ′x˜lnμ′x˜dx˜,where, μ′x˜=μx˜/∫acμx˜dx˜.To obtain the probability density function of equivalent normal random variables, it is necessary to determine the fuzzy variableX˜. According to reference [22], the mean and standard deviation of a are equivalent to the normal random variable x and the fuzzy variable X˜ have the same probability entropy, PEx=PEx˜Obtained by combining equation (17):(19)σx˜=12πePEx˜−0.5.Becauseμ′x˜ satisfies the property of probability density, it can be called a probability density function:(20)μx˜=∫acx˜μ′x˜dx=∫acx˜μx˜∫acμx˜dxdx˜=1∫acμx˜dx˜∫acx˜μx˜dx˜.The probability density function of the equivalent normal random variable of the fuzzy variable is expressed as(21)fx=12πσx˜e−x−μx˜/2σx˜2.The fuzzy random variable is divided intoK intervals, and the UGF is expressed as(22)μz=∑i=1kpizμi.


## 2.1.1. UGF Representation of Random Variables

Let the probability density function of a random variableX be fx, where the distribution of samples is known and bounded, and it is recorded as a0,an. Divide a0,an , into n subintervals, as shown in Figure 1. The corresponding mean and probability of each subinterval are, respectively:(6)μi=a0+a12,…,an−1+an2,Pi=∫ai−1aifxdx,i=1,2,…,n.Figure 1
Probability distribution function of random variables.The universal generating function of random variables from equation (2) can be expressed as(7)uz=∑i=1npizμi.


## 2.1.2. UGF Representation of Probability Box Variables

Let the probability density function of the probability box variableX⌢ be fx⌢ and the distribution interval be a¯,a¯. To determine the mean value and probability of each interval of X, the interval a¯,a¯ can be divided into l cells on average, which can be expressed as(8)a1¯,a1¯,a2¯,a2¯,…,al¯,al¯.The mean value of each interval is(9)μi=μ1,μ2,…,μl=a1¯+a1¯2,a2¯+a2¯2,…,al¯+al¯2.The probability value corresponding to the mean value of each interval is(10)pi=p1,p2,…,pl.And, the probability value meets the following conditions:(11)p1+p2+⋯+pi=1,piL≤pi≤piUi=1,2,…,l.piU=min∫ai¯ai¯fx⌢dx⌢i=1,2,…,l, piU=max∫ai¯ai¯fx⌢dx⌢i=1,2,…,l thus, the UGF of the probability box variable X⌢ is expressed as follows:(12)μz=∑i=1lpizμi.


## 2.1.3. UGF Representation of Interval Variables

For a variableY, if the upper and lower bounds of the variable are known, it is recorded as Y=b¯,b¯ and the interval Y is divided into r cells, which are recorded as b¯i,b¯i,i=1,2,…,r at this time, the mean value and probability of each interval of interval variable Y can be expressed as(13)bi=b1,b2,…,br=b1¯+b1¯2,b2¯+b2¯2,…,br¯+br¯2,pi=p1,p2,…,pr,Therefore, the UGF of interval variables is expressed as(14)μz=∑i=1rpizbi.


## 2.1.4. UGF Representation of Fuzzy Variables

Fuzzy information is also a common type of information in reliability engineering. Fuzzy information is described by the membership function. Triangular fuzzy number is a common membership function, which is represented byX˜ represents the fuzzy variable, and the triangular fuzzy membership function is shown in Figure 2:Figure 2
Membership function of triangular fuzzy number.(15)μx˜=x˜−ab−a,a≤x˜<b,1,x˜=b,x˜−cb−c,b<x˜≤c,0,others.For the fuzziness problems contained in variables, the existing theoretical methods often use the form of cut sets to transform the membership function into the form of intervals. In this paper, the method of entropy equivalence will be used to transform the fuzziness problems in variables into random problems.Shannon defines the probability entropy of random variables under uncertainty as:(16)PEx=−∫−∞+∞fxlnfxdx.For random variables with normal distribution, if the meanμx and variance σx is known, equation (16) can be transformed into(17)PEx=12+ln2πσx.Fuzzy variableX˜ the probability entropy of a is defined as(18)PEx˜=−∫acμ′x˜lnμ′x˜dx˜,where, μ′x˜=μx˜/∫acμx˜dx˜.To obtain the probability density function of equivalent normal random variables, it is necessary to determine the fuzzy variableX˜. According to reference [22], the mean and standard deviation of a are equivalent to the normal random variable x and the fuzzy variable X˜ have the same probability entropy, PEx=PEx˜Obtained by combining equation (17):(19)σx˜=12πePEx˜−0.5.Becauseμ′x˜ satisfies the property of probability density, it can be called a probability density function:(20)μx˜=∫acx˜μ′x˜dx=∫acx˜μx˜∫acμx˜dxdx˜=1∫acμx˜dx˜∫acx˜μx˜dx˜.The probability density function of the equivalent normal random variable of the fuzzy variable is expressed as(21)fx=12πσx˜e−x−μx˜/2σx˜2.The fuzzy random variable is divided intoK intervals, and the UGF is expressed as(22)μz=∑i=1kpizμi.


## 2.2. Structural Reliability Modeling with Mixed Uncertainty

Variables describing mixed uncertainty usually include: random variableX, fuzzy variable X˜, interval variable Y, probability box variable X⌢, if (X,X˜, Y,X⌢) is a variable that affects the state and nature of the structure, then the functional function of the structure can be expressed as(23)Z=gX,X˜,Y,X⌢.According to the stress strength interference model, whenZ<0 is, the structure fails; when Z>0 the structure is in a safe state; when Z=0, the structure reaches the limit state. In Figure 3, three limit states of the structural function in the case of two-dimensional variables.Figure 3
Structural limit state of two-dimensional variables.In general, the failure rate of the structure is obtained by multiple integration of the joint probability density function of each variable in the region ofZ<0. When the structure contains different types of variables and the number of variables is large, it is very difficult to solve it by multiple integrals. The UGF simplifies the calculation to a certain extent. The structural failure probability under various uncertainties based on the UGF can be expressed as(24)F=PfgX,X˜,Y,X⌢<0.When the number of variables affecting reliability is small, the calculation using equation (24) is relatively simple, but when the number of variables affecting reliability is large, it will lead to the problem of large amount of calculation. Therefore, the approximate dimension reduction method is used to reduce the dimension of the multivariable computing function, obtain the approximate function of the multivariable computing function, and then build a UGF model to solve the failure efficiency.


### 2.2.1. Dimensionality Reduction of Multivariable Function

Rahman and Xu [23] proposed an approximate dimensionality reduction method for multivariable functions. The main idea of the multivariable dimensionality reduction method is to decompose the original multivariable function into the sum of multiple single variables. Multivariable dimensionality reduction methods are mainly divided into two categories: Multivariable dimensionality reduction method based on mean point expansion and multivariable dimensionality reduction method based on the maximum possible failure point. This paper mainly adopts the multivariable dimensionality reduction method based on mean point expansion. Let the mean value of random variables in n intervals μ=μx1,μx2,…,μxn, and then the dimension reduction approximation model of the function gX of the mean value is as follows:(25)gX≅g^X=∑i=1ngμx1,μx2,…,Xi,μxn−n−1gμx1,μx2,…,μxn.For variables with different types, the mean values of (X,X˜, Y, X⌢) are, respectively, μx=μx1,μx2,…,μxn, μx˜=μx˜1,μx˜2,…,μx˜l, μY=μY1,μY2,…,μYr, μx⌢=μx⌢1,μx⌢2,…,μx⌢k it can be seen from equation (25) that the dimension reduction approximation of the functional function under a variety of mixed variables can be expressed as(26)gX,X˜,Y,X⌢≅g^X,X˜,Y,X⌢=g1X+g2X˜+g3Y+g4X⌢.If the fuzzy variable is transformed into a normal random variable, based on the equivalent method of constant entropy, equation (26) can be rewritten as(27)gX,X˜,Y,X⌢≅g^X,X˜,Y,X⌢=g⇀X+g3Y+g4X⌢,where g⇀X=g1X+g2X˜−n+l+r+k−1gμX,μX˜,μY,μX⌢.


### 2.2.2. Estimate the Probability Density Function ofZ=g⇀X Based on the Maximum Entropy Estimation Method

Z=g⇀X, X is a random variable, so Z is also a random variable. The maximum entropy estimation method does not need any artificial assumption when approaching the unknown probability density function. It is an effective and high-precision method. In 1957, jayncs proposed the calculation method of maximum entropy. The model of estimating Z=g⇀X with maximum entropy is as follows:(28)maxE=max∫ab−fzlnpzdz.The following constraints are met:(29)∫abfzdz=1,(30)mi=∫abzifzdz,i=1,2,…,s.The value ofmi can be solved by simulation and its model is(31)mi=1N∑j=1Ng⇀Xji.fz is the probability density function and mi is the i-th order origin moment of fz. According to the Lagrange method, the probability density function fz of the maximum entropy distribution is expressed as(32)fz=expλ0+∑i=15λizi,λ0=ln∫abexp∑i=1sλizidz.λ0, λi, i=1,2,…,s is a Lagrange multiplier which can be obtained by the nonlinear programming method. The probability density function expression of multiple variables with the maximum entropy distribution satisfying the constraint conditions in equation (31):(33)fz=expλ0+λ1z+λ2z2+λ3z3+⋯+λszs.Then it is transformed into the form of a universal generating function according to equation (2).Since the functiong3Y is an interval variable, the minimum and maximum values are determined through the linear programming optimization model, that is,(34)g3Lg3U=minmaxgY,s.t.Y1L≤Y1≤Y1U,⋮⋮⋮YrL≤Yr≤YrU.According to equations (26) and (34), the maximum and minimum values of the functional functions gX,X˜,Y,X⌢ in the interval variables are, respectively.(35)gX,X˜,Y,X⌢=gmin=g⇀X+g4X⌢+g3L,gX,X˜,Y,X⌢=gmax=g⇀X+g4X⌢+g3U.The minimum and maximum failure probability values can be obtained by using the calculation method of UGF as shown in Figure4.Figure 4
Flow chart of failure probability procedure.


## 2.2.1. Dimensionality Reduction of Multivariable Function

Rahman and Xu [23] proposed an approximate dimensionality reduction method for multivariable functions. The main idea of the multivariable dimensionality reduction method is to decompose the original multivariable function into the sum of multiple single variables. Multivariable dimensionality reduction methods are mainly divided into two categories: Multivariable dimensionality reduction method based on mean point expansion and multivariable dimensionality reduction method based on the maximum possible failure point. This paper mainly adopts the multivariable dimensionality reduction method based on mean point expansion. Let the mean value of random variables in n intervals μ=μx1,μx2,…,μxn, and then the dimension reduction approximation model of the function gX of the mean value is as follows:(25)gX≅g^X=∑i=1ngμx1,μx2,…,Xi,μxn−n−1gμx1,μx2,…,μxn.For variables with different types, the mean values of (X,X˜, Y, X⌢) are, respectively, μx=μx1,μx2,…,μxn, μx˜=μx˜1,μx˜2,…,μx˜l, μY=μY1,μY2,…,μYr, μx⌢=μx⌢1,μx⌢2,…,μx⌢k it can be seen from equation (25) that the dimension reduction approximation of the functional function under a variety of mixed variables can be expressed as(26)gX,X˜,Y,X⌢≅g^X,X˜,Y,X⌢=g1X+g2X˜+g3Y+g4X⌢.If the fuzzy variable is transformed into a normal random variable, based on the equivalent method of constant entropy, equation (26) can be rewritten as(27)gX,X˜,Y,X⌢≅g^X,X˜,Y,X⌢=g⇀X+g3Y+g4X⌢,where g⇀X=g1X+g2X˜−n+l+r+k−1gμX,μX˜,μY,μX⌢.


## 2.2.2. Estimate the Probability Density Function ofZ=g⇀X Based on the Maximum Entropy Estimation Method

Z=g⇀X, X is a random variable, so Z is also a random variable. The maximum entropy estimation method does not need any artificial assumption when approaching the unknown probability density function. It is an effective and high-precision method. In 1957, jayncs proposed the calculation method of maximum entropy. The model of estimating Z=g⇀X with maximum entropy is as follows:(28)maxE=max∫ab−fzlnpzdz.The following constraints are met:(29)∫abfzdz=1,(30)mi=∫abzifzdz,i=1,2,…,s.The value ofmi can be solved by simulation and its model is(31)mi=1N∑j=1Ng⇀Xji.fz is the probability density function and mi is the i-th order origin moment of fz. According to the Lagrange method, the probability density function fz of the maximum entropy distribution is expressed as(32)fz=expλ0+∑i=15λizi,λ0=ln∫abexp∑i=1sλizidz.λ0, λi, i=1,2,…,s is a Lagrange multiplier which can be obtained by the nonlinear programming method. The probability density function expression of multiple variables with the maximum entropy distribution satisfying the constraint conditions in equation (31):(33)fz=expλ0+λ1z+λ2z2+λ3z3+⋯+λszs.Then it is transformed into the form of a universal generating function according to equation (2).Since the functiong3Y is an interval variable, the minimum and maximum values are determined through the linear programming optimization model, that is,(34)g3Lg3U=minmaxgY,s.t.Y1L≤Y1≤Y1U,⋮⋮⋮YrL≤Yr≤YrU.According to equations (26) and (34), the maximum and minimum values of the functional functions gX,X˜,Y,X⌢ in the interval variables are, respectively.(35)gX,X˜,Y,X⌢=gmin=g⇀X+g4X⌢+g3L,gX,X˜,Y,X⌢=gmax=g⇀X+g4X⌢+g3U.The minimum and maximum failure probability values can be obtained by using the calculation method of UGF as shown in Figure4.Figure 4
Flow chart of failure probability procedure.


## 3. Fatigue Strength Reliability Analysis of Milling Shaft Based on UGF

### 3.1. Reliability Analysis of Fatigue Strength of Sub Milling Shaft

For the heavy CNC boring and milling machine, the milling shaft is the core component of the machine tool. When the machine tool is working, the motor drives the milling shaft to rotate through the gear transmission mechanism, and the milling shaft bears the meshing force from the external gear in the process of rotation; at the same time, the milling shaft bears the cutting force, which is transmitted to the milling shaft by the milling cutter. There is little difference between the maximum outer diameter and the minimum inner diameter of the milling shaft, so the milling shaft can be simplified as a hollow shaft with an inner diameterd2 of 130 mm and an outer diameter d1 of 221 mm. Under the combination of bending and torsion, the stress diagram of the milling shaft is shown in Figure 5:Figure 5
Stress analysis of milling shaft.Fc is the external load (kN); M1 and M2 are the bending moment on the dangerously section (kN⋅m); Fr is the radial force at the meshing part of the gear, and Ft is the circumferential force at the meshing part of the gear; a, b, and c are the dimensions of the milling shaft (m).The functional function of the structure under the fatigue strength failure mode of the milling shaft is as follows:(36)F=Prσ−1K−MH2+MV2+αMT2W<0,where MT—torque of milling shaft (kN · m); MH—bending moment on the horizontal plane of milling shaft (kN · m); MV—bending moment on the vertical plane of milling shaft (kNm); W—bending section coefficient m3; K—the correction coefficient is related to the surface quality and size of the milling shaft; σ−1—fatigue strength of standard smooth specimen MPa;The statistical table of fatigue strength related parameters of machine tool milling shaft is shown in Table1:Table 1
Information on fatigue strength variables.
VariableParameter 1Parameter 2Distribution typeσ−1100025Normal distributionMT2.050.1Normal distributionMH0.45,0.550.1Probability box distributionMV0.4,0.420.1Probability box distributionW0.00210.0023Interval distributionα0.50.6Interval distributionK1.31.5Interval distributionFor random variables, parameter 1 is the mean of normal distribution and parameter 2 is the variance of normal distribution; for the probability box variable, parameter 1 is the mean interval of normal distribution, and parameter 2 is the variance of normal distribution; for interval variables, parameter 1 is the interval lower limit, marked asa, and parameter 2 is the interval upper limit, marked as b.There are many dependent variables. First, approximate dimensionality reduction ofgσ−1,MT,MH,MV,W,α,K at the mean value:(37)gσ−1,MT,MH,MV,W,α,K=σ−11.4−0.42+0.55×MT21/2+1649.35−MH2+1.441/20.0022−MV2+1.521/20.0022+1000K−1.3W−0.42+2.05α1/20.0022.For the random variableσ−1, MT, let z=σ−1/1.4−0.42+0.55×MT21/2+1649.35, the first four order origin moments of Z are obtained by (30), where the first four order origin moments are shown in Table 2:Table 2
First four order origin moments of random variableZ.
First order momentSecond order momentThird order momentFourth order moment1.7724×1033.1421×1065.5717×1099.8825×1012The probability density function ofZ obtained from (33) is(38)fz=exp−1001.257597+0.749403×z+1.142073×10−7×z2+3.000513×10−10×z3−3.377152×10−11×z4.DivideZ into 5 cells, and the mean value and probability of each interval are(39)z¯=1715,1745,1775,1805,1835,pz=0.060868,0.260884,0.404255,0.221053,0.041135..For interval variablesK, W, α, let z1=1000/K−1.3/W−0.42+2.05α1/2/0.0022 the maximum and minimum value of z1 can be obtained through optimization.(40)minmaxz1=1000K−1.3W−0.42+2.05α1/20.0022,1.3≤K≤1.5,0.0021≤W≤0.0023,0.5≤α≤0.6.According to the3σ principle, the distribution interval μ¯−3σ,μ¯+3σ of the probability box variable MH is divided into 7 cells, and the mean value and probability of each cell of MH are obtained:(41)M¯H=0.2,0.3,0.4,0.5,0.6,0.7,0.8,pMH=pMH1,pMH2,…,pMH7,pMH1+pMH2+⋯+pMH7=1pMHj∈0.0013,0.0214,0.0214,0.1359,0.1359,0.3413,0.3413,0.34130.1359,0.3413,0.0214,0.1359,0.0013,0.0214..Similarly,MV is divided into 5 cells, and the mean value and probability of each cell of MV are obtained:(42)M¯V=0.162,0.286,0.410,0.534,0.658,pMV=pMV1,pMV2,…,pMV5,pMV1+pMV2+⋯+pMV5=1,pMVj∈0.0243,0.0379,0.2108,0.2623,0.4627,0.4647,0.2108,0.2623,0.0243,0.0379.According to equation (24), the calculation program is compiled with MATLAB software to obtain the minimum and maximum values of the fatigue strength and failure probability of the milling shaft. The calculation results with Monte Carlo method are shown in Table 3:Table 3
Calculation results of failure probability calculated by different methods.
Failure probabilityPaper methodMonte Carlo methodPfL0.3586×10−40.43×10−4PfU0.04560.03272The fatigue strength of the milling shaft is calculated by using the multivariable dimension reduction method and the UGF. Compared with the results calculated by the Monte Carlo method, the error of the two methods is small, and the method proposed in this paper does not need a large number of data points and has high efficiency.


### 3.2. Sensitivity Analysis of Milling Shaft

If the reliability is sensitive to the changes of some factors, attention should be paid to the design, manufacturing, and processing of the structure to strictly ensure its accuracy. On the contrary, if the uncertainty of variables has no obvious impact on the structural reliability, it can be treated as a fixed value in the analysis, to reduce the difficulty of analysis and improve the calculation efficiency. Reliability sensitivity analysis can provide useful guidance for reliability design, construction, and maintenance of structures. The local sensitivity method and global sensitivity analysis method are used to analyze the fatigue strength of the milling shaft, and then the effects of several important variables that have the greatest impact on the fatigue strength of the milling shaft on the functional function are determined.


#### 3.2.1. Local Sensitivity Analysis Based on Monte Carlo

Local sensitivity is the value of the partial derivative of the input variable to the failure probability at a given parameter value. When the Monte Carlo method is used for analysis, the failure probability is expressed as(43)Pf=1N∑i=1NIg<0.I⋅ is a decision function, and its value can be expressed as follows:(44)g<0,I⋅=1,g≥0,I⋅=0.According to the definition of sensitivity, the expression is(45)∂Pf∂θx=∫…∫∂fx∂θxdx.θx Is the characteristic parameter of the random variable distribution, and fx is the probability density function of the basic random variable.By transforming the integral in equation (45) into the expression of mathematical expectation, the Monte Carlo method can be used to calculate the structural failure sensitivity as follows:(46)∂Pf∂θx=∫…∫∂fx∂θx1fxfxdx=∫…∫I⋅∂fx∂θx1fxfxdx=EI⋅fx∂fx∂θx=1N∑i=1NI⋅fx∂fx∂σx.The sensitivity analysis of the fatigue strength limit state equation of the milling shaft is carried out in combination with the Monte Carlo method. Takingσ−1,MT,MV,MH as random variables and K,W,α as interval variables, the local sensitivity analysis results are shown in Table 4:Table 4
Local sensitivity analysis results of fatigue strength of milling shaft.
ParameterLocal sensitivityσ−13.32×10−8−1.42×10−8MT−1.61×10−51.59×10−5MH−1.05×10−51.00×10−6MV−4.72×10−67.77×10−6W−0.0050.005α1×10−5−1×10−5K−3.33×10−63.33×10−6If the sensitivity is positive, increasingθx value will reduce the failure probability, negative sensitivity and increase θx will increase the failure probability.


#### 3.2.2. Global Sensitivity Analysis Based on Variance

The local sensitivity analysis method is used to test the influence of the change of a single parameter at the nominal value point on the model results, while the global sensitivity analysis extends the change range of the parameters to the whole definition domain. In this paper, the variance sensitivity method is mainly used to analyze the global sensitivity.Variance sensitivity analysis uses variance to describe the uncertainty of input variables. It can allocate the variance of output variables to each input variable and the interaction between different input variables, and decompose the response functionY=gX into(47)Y=gX=g0+∑i=1ngiXi+∑1≤i<j≤ngi,jXi,Xj+L+g1,2,…,nX1,…,Xn.In this formula,g0=EY, gi=EY|Xi−g0, gi,j=EY|Xi,Xj−g0−gi−gj. The variance of the output variable y can be decomposed into(48)VY=∑i=1nVi+∑1≤i<j≤nVi,j+⋯+V1,2,…,n.In this formula,Vi=Vgi=VEY|Xi, Vij=Vgij=VEY|Xi,Xj−Vi−Vj etc.Sobol defines the sensitivity index as(49)Si=ViVY=VEY|XiVY.Si indicate the influence degree of input variable Xi on the output variable Y.The reliability sensitivity analysis of Heavy CNC is carried out in combination with the variance sensitivity analysis method, and the fatigue strength limit state equation of milling shaft is as follows:(50)gσ−1,MT,MV,MH,K,W,α=σ−1K−MH2+MV2+αMT2W.The fatigue strength variance of the milling shaft is calculated by using equations (47) and (48), as shown in Table 5.Table 5
Variance and partial variance.
VpfVσ−1VMTVMHVMVVWVαVK4438.82342.95453.41355.99209.37249.62525.642297.91Apply equation (48) to obtain the main effect sensitivity of each variable, as shown in Table 6.Table 6
Global sensitivity analysis results based on main effect.
Basic variableσ−1MTMHMVWαKMain effect sensitivitySi0.0770.1020.0800.0470.0560.1180.517It can be seen from the table that the uncertainty of the correction coefficientK has the greatest impact on the fatigue strength of the milling shaft, and the others are α, MT, MH, σ−1,W, and MV in turn.


## 3.1. Reliability Analysis of Fatigue Strength of Sub Milling Shaft

For the heavy CNC boring and milling machine, the milling shaft is the core component of the machine tool. When the machine tool is working, the motor drives the milling shaft to rotate through the gear transmission mechanism, and the milling shaft bears the meshing force from the external gear in the process of rotation; at the same time, the milling shaft bears the cutting force, which is transmitted to the milling shaft by the milling cutter. There is little difference between the maximum outer diameter and the minimum inner diameter of the milling shaft, so the milling shaft can be simplified as a hollow shaft with an inner diameterd2 of 130 mm and an outer diameter d1 of 221 mm. Under the combination of bending and torsion, the stress diagram of the milling shaft is shown in Figure 5:Figure 5
Stress analysis of milling shaft.Fc is the external load (kN); M1 and M2 are the bending moment on the dangerously section (kN⋅m); Fr is the radial force at the meshing part of the gear, and Ft is the circumferential force at the meshing part of the gear; a, b, and c are the dimensions of the milling shaft (m).The functional function of the structure under the fatigue strength failure mode of the milling shaft is as follows:(36)F=Prσ−1K−MH2+MV2+αMT2W<0,where MT—torque of milling shaft (kN · m); MH—bending moment on the horizontal plane of milling shaft (kN · m); MV—bending moment on the vertical plane of milling shaft (kNm); W—bending section coefficient m3; K—the correction coefficient is related to the surface quality and size of the milling shaft; σ−1—fatigue strength of standard smooth specimen MPa;The statistical table of fatigue strength related parameters of machine tool milling shaft is shown in Table1:Table 1
Information on fatigue strength variables.
VariableParameter 1Parameter 2Distribution typeσ−1100025Normal distributionMT2.050.1Normal distributionMH0.45,0.550.1Probability box distributionMV0.4,0.420.1Probability box distributionW0.00210.0023Interval distributionα0.50.6Interval distributionK1.31.5Interval distributionFor random variables, parameter 1 is the mean of normal distribution and parameter 2 is the variance of normal distribution; for the probability box variable, parameter 1 is the mean interval of normal distribution, and parameter 2 is the variance of normal distribution; for interval variables, parameter 1 is the interval lower limit, marked asa, and parameter 2 is the interval upper limit, marked as b.There are many dependent variables. First, approximate dimensionality reduction ofgσ−1,MT,MH,MV,W,α,K at the mean value:(37)gσ−1,MT,MH,MV,W,α,K=σ−11.4−0.42+0.55×MT21/2+1649.35−MH2+1.441/20.0022−MV2+1.521/20.0022+1000K−1.3W−0.42+2.05α1/20.0022.For the random variableσ−1, MT, let z=σ−1/1.4−0.42+0.55×MT21/2+1649.35, the first four order origin moments of Z are obtained by (30), where the first four order origin moments are shown in Table 2:Table 2
First four order origin moments of random variableZ.
First order momentSecond order momentThird order momentFourth order moment1.7724×1033.1421×1065.5717×1099.8825×1012The probability density function ofZ obtained from (33) is(38)fz=exp−1001.257597+0.749403×z+1.142073×10−7×z2+3.000513×10−10×z3−3.377152×10−11×z4.DivideZ into 5 cells, and the mean value and probability of each interval are(39)z¯=1715,1745,1775,1805,1835,pz=0.060868,0.260884,0.404255,0.221053,0.041135..For interval variablesK, W, α, let z1=1000/K−1.3/W−0.42+2.05α1/2/0.0022 the maximum and minimum value of z1 can be obtained through optimization.(40)minmaxz1=1000K−1.3W−0.42+2.05α1/20.0022,1.3≤K≤1.5,0.0021≤W≤0.0023,0.5≤α≤0.6.According to the3σ principle, the distribution interval μ¯−3σ,μ¯+3σ of the probability box variable MH is divided into 7 cells, and the mean value and probability of each cell of MH are obtained:(41)M¯H=0.2,0.3,0.4,0.5,0.6,0.7,0.8,pMH=pMH1,pMH2,…,pMH7,pMH1+pMH2+⋯+pMH7=1pMHj∈0.0013,0.0214,0.0214,0.1359,0.1359,0.3413,0.3413,0.34130.1359,0.3413,0.0214,0.1359,0.0013,0.0214..Similarly,MV is divided into 5 cells, and the mean value and probability of each cell of MV are obtained:(42)M¯V=0.162,0.286,0.410,0.534,0.658,pMV=pMV1,pMV2,…,pMV5,pMV1+pMV2+⋯+pMV5=1,pMVj∈0.0243,0.0379,0.2108,0.2623,0.4627,0.4647,0.2108,0.2623,0.0243,0.0379.According to equation (24), the calculation program is compiled with MATLAB software to obtain the minimum and maximum values of the fatigue strength and failure probability of the milling shaft. The calculation results with Monte Carlo method are shown in Table 3:Table 3
Calculation results of failure probability calculated by different methods.
Failure probabilityPaper methodMonte Carlo methodPfL0.3586×10−40.43×10−4PfU0.04560.03272The fatigue strength of the milling shaft is calculated by using the multivariable dimension reduction method and the UGF. Compared with the results calculated by the Monte Carlo method, the error of the two methods is small, and the method proposed in this paper does not need a large number of data points and has high efficiency.


## 3.2. Sensitivity Analysis of Milling Shaft

If the reliability is sensitive to the changes of some factors, attention should be paid to the design, manufacturing, and processing of the structure to strictly ensure its accuracy. On the contrary, if the uncertainty of variables has no obvious impact on the structural reliability, it can be treated as a fixed value in the analysis, to reduce the difficulty of analysis and improve the calculation efficiency. Reliability sensitivity analysis can provide useful guidance for reliability design, construction, and maintenance of structures. The local sensitivity method and global sensitivity analysis method are used to analyze the fatigue strength of the milling shaft, and then the effects of several important variables that have the greatest impact on the fatigue strength of the milling shaft on the functional function are determined.


### 3.2.1. Local Sensitivity Analysis Based on Monte Carlo

Local sensitivity is the value of the partial derivative of the input variable to the failure probability at a given parameter value. When the Monte Carlo method is used for analysis, the failure probability is expressed as(43)Pf=1N∑i=1NIg<0.I⋅ is a decision function, and its value can be expressed as follows:(44)g<0,I⋅=1,g≥0,I⋅=0.According to the definition of sensitivity, the expression is(45)∂Pf∂θx=∫…∫∂fx∂θxdx.θx Is the characteristic parameter of the random variable distribution, and fx is the probability density function of the basic random variable.By transforming the integral in equation (45) into the expression of mathematical expectation, the Monte Carlo method can be used to calculate the structural failure sensitivity as follows:(46)∂Pf∂θx=∫…∫∂fx∂θx1fxfxdx=∫…∫I⋅∂fx∂θx1fxfxdx=EI⋅fx∂fx∂θx=1N∑i=1NI⋅fx∂fx∂σx.The sensitivity analysis of the fatigue strength limit state equation of the milling shaft is carried out in combination with the Monte Carlo method. Takingσ−1,MT,MV,MH as random variables and K,W,α as interval variables, the local sensitivity analysis results are shown in Table 4:Table 4
Local sensitivity analysis results of fatigue strength of milling shaft.
ParameterLocal sensitivityσ−13.32×10−8−1.42×10−8MT−1.61×10−51.59×10−5MH−1.05×10−51.00×10−6MV−4.72×10−67.77×10−6W−0.0050.005α1×10−5−1×10−5K−3.33×10−63.33×10−6If the sensitivity is positive, increasingθx value will reduce the failure probability, negative sensitivity and increase θx will increase the failure probability.


### 3.2.2. Global Sensitivity Analysis Based on Variance

The local sensitivity analysis method is used to test the influence of the change of a single parameter at the nominal value point on the model results, while the global sensitivity analysis extends the change range of the parameters to the whole definition domain. In this paper, the variance sensitivity method is mainly used to analyze the global sensitivity.Variance sensitivity analysis uses variance to describe the uncertainty of input variables. It can allocate the variance of output variables to each input variable and the interaction between different input variables, and decompose the response functionY=gX into(47)Y=gX=g0+∑i=1ngiXi+∑1≤i<j≤ngi,jXi,Xj+L+g1,2,…,nX1,…,Xn.In this formula,g0=EY, gi=EY|Xi−g0, gi,j=EY|Xi,Xj−g0−gi−gj. The variance of the output variable y can be decomposed into(48)VY=∑i=1nVi+∑1≤i<j≤nVi,j+⋯+V1,2,…,n.In this formula,Vi=Vgi=VEY|Xi, Vij=Vgij=VEY|Xi,Xj−Vi−Vj etc.Sobol defines the sensitivity index as(49)Si=ViVY=VEY|XiVY.Si indicate the influence degree of input variable Xi on the output variable Y.The reliability sensitivity analysis of Heavy CNC is carried out in combination with the variance sensitivity analysis method, and the fatigue strength limit state equation of milling shaft is as follows:(50)gσ−1,MT,MV,MH,K,W,α=σ−1K−MH2+MV2+αMT2W.The fatigue strength variance of the milling shaft is calculated by using equations (47) and (48), as shown in Table 5.Table 5
Variance and partial variance.
VpfVσ−1VMTVMHVMVVWVαVK4438.82342.95453.41355.99209.37249.62525.642297.91Apply equation (48) to obtain the main effect sensitivity of each variable, as shown in Table 6.Table 6
Global sensitivity analysis results based on main effect.
Basic variableσ−1MTMHMVWαKMain effect sensitivitySi0.0770.1020.0800.0470.0560.1180.517It can be seen from the table that the uncertainty of the correction coefficientK has the greatest impact on the fatigue strength of the milling shaft, and the others are α, MT, MH, σ−1,W, and MV in turn.


## 3.2.1. Local Sensitivity Analysis Based on Monte Carlo

Local sensitivity is the value of the partial derivative of the input variable to the failure probability at a given parameter value. When the Monte Carlo method is used for analysis, the failure probability is expressed as(43)Pf=1N∑i=1NIg<0.I⋅ is a decision function, and its value can be expressed as follows:(44)g<0,I⋅=1,g≥0,I⋅=0.According to the definition of sensitivity, the expression is(45)∂Pf∂θx=∫…∫∂fx∂θxdx.θx Is the characteristic parameter of the random variable distribution, and fx is the probability density function of the basic random variable.By transforming the integral in equation (45) into the expression of mathematical expectation, the Monte Carlo method can be used to calculate the structural failure sensitivity as follows:(46)∂Pf∂θx=∫…∫∂fx∂θx1fxfxdx=∫…∫I⋅∂fx∂θx1fxfxdx=EI⋅fx∂fx∂θx=1N∑i=1NI⋅fx∂fx∂σx.The sensitivity analysis of the fatigue strength limit state equation of the milling shaft is carried out in combination with the Monte Carlo method. Takingσ−1,MT,MV,MH as random variables and K,W,α as interval variables, the local sensitivity analysis results are shown in Table 4:Table 4
Local sensitivity analysis results of fatigue strength of milling shaft.
ParameterLocal sensitivityσ−13.32×10−8−1.42×10−8MT−1.61×10−51.59×10−5MH−1.05×10−51.00×10−6MV−4.72×10−67.77×10−6W−0.0050.005α1×10−5−1×10−5K−3.33×10−63.33×10−6If the sensitivity is positive, increasingθx value will reduce the failure probability, negative sensitivity and increase θx will increase the failure probability.


## 3.2.2. Global Sensitivity Analysis Based on Variance

The local sensitivity analysis method is used to test the influence of the change of a single parameter at the nominal value point on the model results, while the global sensitivity analysis extends the change range of the parameters to the whole definition domain. In this paper, the variance sensitivity method is mainly used to analyze the global sensitivity.Variance sensitivity analysis uses variance to describe the uncertainty of input variables. It can allocate the variance of output variables to each input variable and the interaction between different input variables, and decompose the response functionY=gX into(47)Y=gX=g0+∑i=1ngiXi+∑1≤i<j≤ngi,jXi,Xj+L+g1,2,…,nX1,…,Xn.In this formula,g0=EY, gi=EY|Xi−g0, gi,j=EY|Xi,Xj−g0−gi−gj. The variance of the output variable y can be decomposed into(48)VY=∑i=1nVi+∑1≤i<j≤nVi,j+⋯+V1,2,…,n.In this formula,Vi=Vgi=VEY|Xi, Vij=Vgij=VEY|Xi,Xj−Vi−Vj etc.Sobol defines the sensitivity index as(49)Si=ViVY=VEY|XiVY.Si indicate the influence degree of input variable Xi on the output variable Y.The reliability sensitivity analysis of Heavy CNC is carried out in combination with the variance sensitivity analysis method, and the fatigue strength limit state equation of milling shaft is as follows:(50)gσ−1,MT,MV,MH,K,W,α=σ−1K−MH2+MV2+αMT2W.The fatigue strength variance of the milling shaft is calculated by using equations (47) and (48), as shown in Table 5.Table 5
Variance and partial variance.
VpfVσ−1VMTVMHVMVVWVαVK4438.82342.95453.41355.99209.37249.62525.642297.91Apply equation (48) to obtain the main effect sensitivity of each variable, as shown in Table 6.Table 6
Global sensitivity analysis results based on main effect.
Basic variableσ−1MTMHMVWαKMain effect sensitivitySi0.0770.1020.0800.0470.0560.1180.517It can be seen from the table that the uncertainty of the correction coefficientK has the greatest impact on the fatigue strength of the milling shaft, and the others are α, MT, MH, σ−1,W, and MV in turn.


## 4. Conclusions

A reliability model based on UGF is established. An example shows that the method is effective.The results of local sensitivity analysis are positive and negative. Positive sensitivity shows that the impact on reliability is positive, and negative sensitivity has a negative impact on reliability. VariableMTMVMH mean sensitivity is negative and variance sensitivity is positive, indicating that increasing the mean will increase the failure probability, and increasing the variance will reduce the failure probability. Increasing the mean value of variable σ−1 decreases the failure probability, and increasing the variance of variable σ−1 increases the failure probability. For interval variables, the value of K and α shall be increased and the value of W shall be reduced to reduce the failure probability.In the global sensitivity analysis results, the main effect sensitivity index value of the comprehensive correction coefficientK is the largest, followed by the stress correction coefficient α, the bending moment MH on the horizontal plane, the torque MT on the milling axis and the fatigue strength of the standard smooth specimen σ−1 are the second, and the bending moment MV and bending section coefficient W on the vertical plane are the smallest. Therefore, the uncertainty of the correction coefficient K has a great impact on the uncertainty of the fatigue strength of the milling shaft. Increasing the correction coefficient K can reduce the failure probability of the fatigue strength of the milling shaft the most.



---

*Source: 3756824-2022-03-19.xml*