Risk Evaluation Method Based on Fault Propagation and Diffusion

Abstract

The high reliability demand of the machining center emphasizes the accuracy of the fault risk evaluation. In the traditional fault risk evaluation research of the machining center, the influence of fault mode is mostly based on subjective recommendation or does not consider the propagation and diffusion of fault, which makes the risk evaluation results different from the real situation. Therefore, this paper presents a framework to evaluate the fault risk for machining center components. A certain type of machining center is considered as a case study. The fault mode frequency ratio of components is calculated by fault mode analysis. The fault rate calculation is conducted based on the Johnson method. Considering that different fault modes have different influences on fault propagation breadth and depth, the hypergraph theory is used to build a hypernetwork model. The propagation and diffusion influence degree are defined to describe the propagation and diffusion process of faults. Then, the comprehensive influence degree of fault mode is calculated. The risk evaluation is realized by considering the component fault rate, fault mode frequency ratio, and the comprehensive influence degree of fault mode. The method proposed in this paper can provide a reference for the formulation of risk strategies for the machining center.

1. Introduction

As the “core equipment” of the manufacturing industry, the machining center has the characteristics of high precision, high efficiency, and high flexibility [1]. Affected by the internal structure and the external working conditions, its use process is more and more complicated. Minor damage to the components will cause a series of cascade reactions in the whole machine [2]. If the component fails and cannot be handled in time, it will affect the whole production schedule. In serious cases, it may even threaten the personal safety of on-site operators. This puts forward higher requirements for the reliability and safety of the machining center. Therefore, risk evaluation plays a significant role in ensuring the safe and stable operation of the machining center.

Risk is the possibility of encountering hazards in the future [3,4,5]. Risk evaluation is not only to test, measure, and estimate risks but also to reprocess the results of risk identification. At present, the commonly used risk evaluation methods for complex systems are based on expert knowledge [6], fault mode and effect analysis (FMEA) [7], fault mode effect and criticality analysis (FMECA) [8], fault tree [9], probabilistic risk analysis [10], Bayesian network [11], and other methods [12].

Combining the functional fault identification and propagation of the system with fuzzy logic, Wang et al. proposed a functional fault risk analysis method for natural gas fracturing system based on expert knowledge [13]. Wang et al. improved the FMEA method and constructed a hybrid FMEA risk analysis framework by applying the cloud model, Choquet integral, and gain–loss dominance scoring method [14]. This method considered the risk of groups and individuals to prioritize fault risks and carried out machine tool risk analysis accordingly. On the basis of the traditional FMEA method, Li et al. transformed the method into a cloud model and then applied the best–worst method (BWM) to obtain the cloud weight [15]. Combined with the cloud distance to evaluate the expert’s weight on the fault mode risk factors, the comprehensive cloud expression of each fault mode was obtained to evaluate the risk of the machine tool. Wang et al. applied multi-criteria decision-making and probabilistic hesitant fuzzy linguistic term sets to improve traditional FMEA [16]. This method increased the cognitive uncertainty and group risk and then achieved the risk evaluation of fault modes via expert groups. Yucesan et al. applied the fuzzy best–worst method to weigh the risk parameters in FMEA and integrated the Fuzzy Bayesian Network (FBN) with FMEA to clarify the correlation between equipment fault events [17]. Then, they carried out the risk analysis of industrial equipment. Shen et al. proposed an improved FMEA method based on multi-attribute group decision making for the problems of unreasonable allocation of expert weights [18]. Rouabhia-Essalhi, Boukrouh, and Ghemari applied the FMECA method in the risk analysis of industrial equipment [19]. By taking measures to reduce the risk priority number to an acceptable risk value, the potential fault of equipment was reduced, thereby reducing the downtime of equipment. Aiming at the lack of evaluation data and poor interpretability in the risk evaluation process of complex systems, Zhu et al. constructed a risk evaluation model based on belief rule base (BRB) and fault tree analysis (FTA) [20]. To identify system risks, Kwag, Gupta, and Dinh constructed a performance-based probabilistic risk evaluation method, which uses Bayesian statistics to clarify the correlation between non-Boolean relationships and different levels of system events [21]. Considering that many factors will affect the reliability of high-speed trains, Zeng et al. constructed a risk evaluation model based on a Bayesian network array to reduce the impact of faults on the operation reliability of trains [22]. Based on fault tree analysis, integrated fuzzy set theory, Bayesian network, and Markov chain, Mohammadi et al. conducted the top event probability and reliability evaluation of cranes to determine the appropriate risk evaluation methods [23].

As described in the above literature, the research on risk evaluation of complex systems has made important progress, but it also faces certain challenges. In the process of fault risk evaluation of complex systems, the fault propagation breadth and depth of different modes are different, and there are certain functional dependencies between components. It is necessary to further explore the influence of fault modes. In addition, most of the risk evaluation models established in the past are static models, which cannot realize real-time dynamic evaluation of fault risk. Therefore, this paper proposes a novel framework for risk evaluation. A case study on the CNC machine center is conducted to elaborate the framework. The framework integrates the fault mode frequency ratio, the fault rate function, and the comprehensive influence degree of fault mode to dynamically evaluate the component risk.

The main contributions of this paper are as follows: (1) A framework to evaluate the risk of the CNC machining center component is presented. (2) A component fault mode hypernetwork model is constructed based on the hypergraph theory, considering the heterogeneity of the network between components and the network between fault modes. (3) A calculation method of comprehensive influence of fault mode is developed, which takes into account the propagation and diffusion of faults. (4) A case study using the real data of a CNC machining center is conducted, and we finally realized the evaluation of the risk of a certain type of CNC machining center at different times.

The rest of this paper is organized as follows: Section 2 presents the framework for risk evaluation. Section 3 elaborates on the fault mode analysis. Section 4 calculates the fault rate of the component. Section 5 provides the method of comprehensive influence of fault mode. In Section 6, the risk evaluation of the component is conducted. Conclusions are given in Section 7.

2. A Framework for Risk Evaluation

The machining center considered in this paper is a typical manufacturing equipment. Its healthy operation is crucial for the entire production system. We collected fault data from customers who purchased this type of machining center. The time span of the fault data collection was one year. Figure 1 presents a framework for risk evaluation, and the main technologies included in this framework are as follows:

(1)

Fault mode analysis. Fault mode analysis is the first step in developing the risk evaluation of machining center components. The fault frequency and fault mode of the components are recorded during the operation of the machining center. Then, we calculate the frequency ratio of component fault modes. This step acts as a basis for risk evaluation of the machining center.

(2)

Fault rate calculation. Fault rate acts as a core role in component risk evaluation. We use the Johnson method to correct the fault data rank. The two-parameter Weibull model fits in as a hypothetical distribution of the fault data. The least-squares method can obtain accurate parameter estimates. Eventually, the fault rate function of each component is determined.

(3)

Comprehensive influence degree of fault mode. The comprehensive fault mode influence degree is to reflect the influence degree of each fault mode more accurately. First, we construct a hypernetwork model of the component fault mode. Second, the propagation and diffusion influence of fault modes is calculated. Finally, the comprehensive influence degree of the fault mode can be obtained.

(4)

Risk evaluation. In this paper, the risk evaluation considers the component fault rate, the frequency ratio of fault mode, and the comprehensive influence degree of fault modes. The accuracy and reliability of the risk evaluation results can provide a reference for the formulation of the risk strategy of the machining center components.

3. Fault Mode Analysis

In this section, we conduct the fault mode analysis of the components. According to the functional structure mapping relationship of the machining center, we divide it into components, and clarify the composition and relationship of each component. We compiled the operation record table and the fault record table of the machining center. The record table mainly includes the start time of the record, fault occurrence time, fault location, fault causes, and solutions. Then, the fault information is collected, summarized, and sorted out. Therefore, we carried out fault mode analysis, explored the fault mechanism, and obtained the fault mode frequency ratio of components.

This paper takes a certain series of machining centers as the research object. It is divided into 11 components: tool magazine (M), spindle system (S), cooling system (W), hydraulic system (D), electrical system (V), worktable (T), pneumatic system (G), XYZ axis feed system (J), lubrication system (L), auxiliary system (K), and numerical control system (NC). Via fault analysis, we draw the fault distribution histogram of the machining center components, as shown in Figure 2.

Via the fault mode analysis, the frequency histogram of the fault mode is obtained, as shown in Figure 3, and the code is detailed in

Table 1. Summary table of fault mode of the machining center.

Number	Mode Code	Fault Mode
1	101	Parts are damaged
2	103	Damage to liquid, gas, oil parts
3	301	Liquid, gas, and oil leakage
4	201	Fasteners are loose
5	203	The pretight parts are loose
6	404	Improper pressure adjustment
7	105	Protective plate and cover are damaged
8	108	The fuse is damaged
9	409	Control line connection is wrong
10	601	Temperature rises too high
⋮	⋮	⋮
29	506	Moving parts

.

As can be found in Figure 3, the most common fault mode of this series of machining centers is 101, 103, and 301. Therefore, it is of great significance to ensure the reliable quality of components for the reliable operation of the machining center.

4. Fault Rate Calculation

The Johnson method is applied to correct the fault time rank of machining center components [24]. It is assumed that the fault data of machining center components obey the Weibull distribution [25]. The reliability modeling process is shown in Figure 4. As shown in Figure 4, the Johnson method is applied to correct the order of the collected fault data, and then a two-parameter Weibull model is constructed. Using the least squares method for parameter estimation.

The relationship between the fault rate function and the reliability function is shown below.

$$\lambda \left(t\right)=\frac{f\left(t\right)}{R\left(t\right)}=\frac{\gamma }{{\theta }^{\gamma }}{t}^{\gamma -1},t\ge 0$$

(1)

where $\lambda \left(t\right)$ is fault rate function, $f\left(t\right)$ is probability density function, $R\left(t\right)$ is reliability function, $t$ is time variable, $\theta$ is scale parameter, and $\gamma$ is shape parameter.

According to the reliability modeling process in Figure 4, combined with the collected fault data of components, the component fault rate function can be obtained, as shown in

Table 2. Fault rate function of components.

Code	Fault Rate
M	$\lambda \left(t\right)=\left(1.312/{2458.200}^{1.312}\right)t^{\left(1.312-1\right)}$
J	$\lambda \left(t\right)=\left(1.295/2606.440^{1.295}\right)t^{\left(1.295-1\right)}$
S	$\lambda \left(t\right)=\left(0.641/{7832.664}^{0.641}\right)t^{\left(0.641-1\right)}$
V	$\lambda \left(t\right)=\left(0.771/{7935.696}^{0.771}\right)t^{\left(0.771-1\right)}$
K	$\lambda \left(t\right)=\left(0.931/{5960.776}^{0.931}\right)t^{\left(0.931-1\right)}$
G	$\lambda \left(t\right)=\left(0.524/{\mathrm{27,503.400}}^{0.524}\right)t^{\left(0.524-1\right)}$
W	$\lambda \left(t\right)=\left(0.910/{4145.720}^{0.910}\right)t^{\left(0.910-1\right)}$
NC	$\lambda \left(t\right)=\left(0.542/{7912.399}^{0.542}\right)t^{\left(0.542-1\right)}$
L	$\lambda \left(t\right)=\left(0.791/{8161.512}^{0.791}\right)t^{\left(0.791-1\right)}$
D	$\lambda \left(t\right)=\left(0.656/{8905.062}^{0.656}\right)t^{\left(0.656-1\right)}$
T	$\lambda \left(t\right)=\left(0.773/{6869.766}^{0.773}\right)t^{\left(0.773-1\right)}$

. Since the fault rate of components is a function of time, different times will correspond to different values. This paper arbitrarily takes 4000 h as an example to illustrate the method. The fault rate of the component at 4000 h can be obtained based on Table 2, as shown in

Table 3. Fault rate of components at 4000 h.

Code	Fault Rate	Code	Fault Rate
V	0.000114	G	0.000048
W	0.000220	J	0.000564
K	0.000161	M	0.000621
D	0.000097	S	0.000104
L	0.000112	T	0.000127
NC	0.000094

.

5. Comprehensive Influence Degree of Fault Mode

In traditional risk evaluation, the influence degree of fault mode is mostly based on the reference value given by recommended standards or quantitative estimation [26]. In this paper, a hypernetwork model is used to calculate the propagation and diffusion value of each fault mode.

5.1. Construction of Component Fault Mode Hypernetwork

The concept of the hypergraph was proposed by C. Berge in 1970. Hypergraphs can be used to analyze multi-layer heterogeneous networks. A network that can be described by a hypergraph is a hypernetwork [27,28]. A component of the machining center has multiple fault modes, and a fault mode may occur in multiple components. The components and fault modes are associated in the form of a set. The component fault mode network is a heterogeneous network composed of two types of nodes with different properties. Therefore, the machining center component fault mode hypernetwork model can be established to describe the interactions and effects between component fault modes.

During the operation of the machining center, different components may have different fault modes, and the same mode may occur multiple times in the same component. Therefore, this paper uses the number of modes occurring in the components as weights to build a weighted component fault mode hypergraph with components as hyperedges and fault modes as nodes.

Defining the component fault mode hypergraph as $H=\left(V,E,C\right)$, then

1.

$V=\left\{v_1,v_2,\cdots ,v_g\right\}$, $g$ is the number of fault modes and $v_i\in V\left(1\le i\le g\right)$ is the fault mode with code $i$, called the node of the hypergraph.

2.

$E=\left\{e_1,e_2,\cdots ,e_h\right\}$ is a subset of $V$, $e_j\in E\left(1\le j\le h\right)$ represents the component of the code $j$, which is called hyperedge, and has the following properties:

(1)

$e_j\subset V\left(e_j\ne \varnothing ,1\le j\le h\right)$, which means that the component can have one or more modes.

(2)

${\displaystyle \underset{j=1}{\overset{h}{\cup }}{e}_j=}V$, the set of elements of all hyperedges in the hypergraph is the set of modes.

(3)

$v_i\in e_j$ means that fault mode $v_i$ is a unit of $e_j$, indicating that the component has the $v_i$th fault mode.

(4)

If $v_i,v_j\in e_j$, node $v_i$ and $v_k$ are adjacent in component $e_j$, denoted as $v_i\stackrel{e_j}{\leftrightarrow }{v}_k$. If there is no component making node $v_i$ adjacent to $v_k$, denoted as $v_i\stackrel{\varnothing }{\leftrightarrow }{v}_k$.

(5)

If $e_i\cap e_j\ne \varnothing$, the hyperedge $e_i$ is adjacent to $e_j$, denoted as $e_i\stackrel{}{\leftrightarrow }{e}_j$, and the adjacent hyperedges have a common node. For $\forall v_i$, $\exists e_j$ make $v_i\in e_j$, there is no isolated fault mode that does not belong to any component.

3.

$C=\left\{c_{\mathit{ij}}|v_i\in e_j,1\le i\le g,1\le j\le h\right\}$ represents the set of times the fault mode occurs in the component. $c_{\mathit{ij}}$ is the weight of the node $v_i$ belonging to the hyperedge $e_j$, if $v_i\in e_j$, the value of $c_{\mathit{ij}}$ is not 0, otherwise it is 0.

4.

The fault mode node $v_i\in V$ is described by the fault mode code $i$ and the fault mode occurrence probability ${\mu }_{i\in j}\left(t\right)$. The larger ${\mu }_{i\in j}\left(t\right)$ is, the greater the probability of the fault mode $v_i$ occurring in the component $e_j$ is.

5.

Component hyperedge $e_j\in E$, described by code $j$ and reliability $R_{e_j}\left(t\right)$. The larger $R_{e_j}\left(t\right)$ is, the higher the reliability of component $e_j$ is, and the smaller the probability of fault is.

6.

The component fault mode hypergraph $H$ can be expressed in two forms: set or matrix. Since the matrix can express the relationship between elements more clearly, this paper uses the form of matrix to describe the relationship between components and fault modes. The correlation matrix corresponding to the hypergraph $H$ is expressed by $\mathit{H}$.

The $h$ row in the matrix $\mathit{H}$ corresponds to the component hyperedge $e_1,e_2,\cdots ,e_h$, and the $g$ column corresponds to the fault mode node $v_1,v_2,\cdots ,v_g$. Let $\mathit{H}\left[i,j\right]=c_{\mathit{ij}}$ denote the element in the $i$th row and $j$th column of the matrix.

The edge load attribute will be an important structural feature of key fault propagation path identification in subsequent systems. From the process of building a component fault mode hypernetwork model, combined with the collected fault information, a machining center weighted component fault mode hypergraph is obtained, as shown in Figure 5.

According to Figure 5, the correlation matrix $\mathit{H}$ can be obtained as follows:

$$\mathit{H}=\stackrel{\begin{array}{ccccccccccccc} & & M& J& S& V& K& G& W& \mathit{NC}& L& D& T\end{array}}{\begin{array}{c}101\\ 301\\ 103\\ 102\\ 605\\ ⋮\\ 108\\ 201\\ 203\\ 404\\ 506\end{array}\left[\begin{array}{ccccccccccc}6& 2& 2& 0& 3& 0& 4& 0& 1& 0& 3\\ 0& 0& 0& 0& 0& 2& 2& 0& 2& 8& 1\\ 0& 1& 0& 0& 0& 1& 1& 0& 4& 5& 1\\ 0& 0& 0& 7& 0& 0& 0& 2& 0& 0& 0\\ 2& 0& 2& 0& 0& 0& 0& 0& 0& 0& 1\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]}$$

(2)

5.2. Propagation and Diffusion Influence Degree of Fault Mode

When one fault mode occurs in a component, it will cause the occurrence of other fault modes of the component, and it will diffuse to other components via the hyperedge connection in the hypernetwork, causing faults. The influence of a certain fault mode on the entire network has always been a concern of researchers. The calculation of the fault influence degree can quantitatively guide the risk evaluation. Therefore, on the basis of establishing the component fault mode hypernetwork, the fault propagation and diffusion influence degree are studied. It should be noted that the following definitions should be made when calculating the fault propagation and diffusion influence degree of components:

1.

The propagation influence of fault mode. If fault mode $v_i$ occurs and there is $v_i\stackrel{e_j}{\leftrightarrow }{v}_k$, the occurrence of fault mode $v_i$ will cause the occurrence of $v_k$, and this process is called the propagation influence of fault mode.

2.

Fault propagation capability. The probability of occurrence of a fault mode in a component that causes other fault modes in the component to occur. The fault propagation capability is related to the occurrence probability of the mode $v_i$.

3.

The influence degree of fault propagation. The propagation influence degree represents the probability that the occurrence of a certain fault mode will cause the occurrence of a specific fault mode. It is related to the propagation capability of the faults.

4.

The diffusion influence of fault mode. If $v_i\stackrel{e_x}{\leftrightarrow }{v}_j$, $v_j\stackrel{e_y}{\leftrightarrow }{v}_k$, and $v_i\stackrel{\varnothing }{\leftrightarrow }{v}_k$, the phenomenon that the occurrence of fault mode $v_i$ eventually leads to the occurrence of mode $v_k$ is called fault mode diffusion influence.

5.

The influence degree of fault diffusion. The diffusion influence degree indicates the probability that the occurrence of a certain fault mode will cause the occurrence of other fault modes.

6.

Fault diffusion trajectory. The diffusion trajectory from $v_1$ to $v_{t+1}$ is composed of a set of fault mode node sequences $\left(v_1,v_2,\cdots ,v_t,v_{t+1}\right)$. Adjacent nodes belong to the same component and have an association relationship because they belong to the same component. Non-adjacent nodes belong to different components. The length of the trajectory from $v_i$ to $v_{i+1}$ is counted as one step, the length of the trajectory from $v_i$ to $v_{i+k}$ is counted as $k$ steps, and $v_{i+k}$ is called the $k$-step reachable node of $v_i$.

7.

The influence degree of fault propagation diffusion. The sum of the propagation and diffusion influences caused by a certain fault mode in the hypernetwork.

(1)

Fault mode propagation influence capability.

The fault propagation influence capability of a fault mode is related to the probability of occurrence of the fault mode and the number of times it occurs in a component. According to the established hypergraph, when a certain fault mode of a component occurs many times during the operation of the machining center, it indicates that the mode of the component has a strong fault propagation capability. The function $\phi \left(v_i,e_j\right)$ represents the fault propagation capability of fault mode $v_i$ in component $e_j$, is shown as follows:

$$\phi \left(v_i,e_j\right)={\mu }_{i\in j}\left(t\right)\times \frac{H\left[i,j\right]}{\sqrt{{\displaystyle \sum _{k=1}^n{\left(H\left[k,j\right]\right)}^2}}},\left(k\ne i\right)$$

(3)

where $0\le \phi \left(v_i,e_j\right)\le 1$. $\mathit{H}\left[i,j\right]=c_{\mathit{ij}}$ indicates the number of times of fault mode $v_i$ appearing in component $e_j$. ${\mu }_i\left(t\right)$ represents the occurrence probability of fault mode $v_i$, and the calculation equation is shown in Equation (4):

$${\mu }_{i\in j}\left(t\right)=\frac{{\alpha }_{\mathit{ji}}\left(t\right)\cdot {\lambda }_j\left(t\right)}{{\lambda }_Z\left(t\right)}$$

(4)

where ${\alpha }_{\mathit{ji}}\left(t\right)$ is the frequency of fault mode $v_i$ in component $e_j$ at time $t$. ${\lambda }_j\left(t\right)$ represents the fault rate of component $e_j$ at time $t$, and ${\lambda }_z\left(t\right)$ represents the fault rate of the whole machine at time $t$.

According to Equation (4), combined with Table 1 and Table 3, the fault mode frequency ratio and fault mode occurrence probability of components at 4000 h are obtained and summarized, as shown in

Table 4. Fault mode frequency ratio and fault mode occurrence probability of components at 4000 h.

Code	Mode	${\mathit{\alpha }}_{\mathit{j}\mathit{i}}$	${\mathit{\mu }}_{\mathit{i}\in \mathit{j}}$	Code	Mode	${\mathit{\alpha }}_{\mathit{j}\mathit{i}}$	${\mathit{\mu }}_{\mathit{i}\in \mathit{j}}$	Code	Mode	${\mathit{\alpha }}_{\mathit{j}\mathit{i}}$	${\mathit{\mu }}_{\mathit{i}\in \mathit{j}}$
M	101	0.3158	0.1403	T	0101	0.3333	0.0303	W	101	0.2857	0.0450
	503	0.1579	0.0701		0103	0.1111	0.0101		302	0.2143	0.0337
	202	0.1053	0.0468		0301	0.1111	0.0101		301	0.1429	0.0225
	605	0.1053	0.0468		0503	0.1111	0.0101		805	0.1429	0.0225
	201	0.0526	0.0234		0603	0.1111	0.0101		103	0.0714	0.0112
	205	0.0526	0.0234		0605	0.1111	0.0101		601	0.0714	0.0112
	501	0.0526	0.0234		0607	0.1111	0.0101		610	0.0714	0.0112
	502	0.0526	0.0234	L	103	0.5714	0.0458	J	101	0.2857	0.1153
	506	0.0526	0.0234		301	0.2857	0.0229		103	0.1429	0.0577
	616	0.0526	0.0234		101	0.1429	0.0114		502	0.1429	0.0577
S	410	0.2143	0.0159	V	102	0.7778	0.0634		601	0.1429	0.0577
	101	0.1429	0.0106		108	0.1111	0.0093		603	0.1429	0.0577
	202	0.1429	0.0106		409	0.1111	0.0093		803	0.1429	0.0577
	605	0.1429	0.0106	K	0101	0.5000	0.0576	G	403	0.4286	0.0146
	203	0.0714	0.0053		0501	0.1667	0.0192		301	0.2857	0.0097
	503	0.0714	0.0053		0205	0.1667	0.0192		103	0.1429	0.0049
	609	0.0714	0.0053		0105	0.1667	0.0192		404	0.1429	0.0049
	603	0.0714	0.0053	D	0301	0.6154	0.0427	NC	102	0.6670	0.0447
	606	0.0714	0.0053		0103	0.3846	0.0267		501	0.3330	0.0223

.

Using the correlation properties of the correlation matrix, the propagation capability matrix $\mathit{B}$ of the fault mode in each component can be calculated as $\mathit{B}\left[i,j\right]=\phi \left(v_i,e_j\right)$. The stronger the propagation capability of a certain fault mode, the greater the probability of causing other fault modes, and the larger the corresponding element value in matrix $\mathit{B}$.

$$B=\stackrel{\begin{array}{cc} & \begin{array}{ccccccc} \hspace{1em}\hspace{1em}\hspace{1em}M \hspace{1em}& \hspace{1em}\hspace{1em}J \hspace{1em}& \hspace{1em}\hspace{1em}S \hspace{1em}& \hspace{1em}\hspace{1em}\cdots \hspace{1em}& \hspace{1em}L \hspace{1em}& \hspace{1em}\hspace{1em}D \hspace{1em}& \hspace{1em}\hspace{1em}T \hspace{1em}\end{array}\end{array}}{\begin{array}{cc}\begin{array}{c}101\\ 301\\ 103\\ ⋮\\ 203\\ 404\\ 506\end{array}& \left[\begin{array}{ccccccc}0.1096& 0.0769& 0.0042& \cdots & 0.0025& 0& 0.0235\\ 0& 0& 0& \cdots & 0.0100& 0.0362& 0.0026\\ 0& 0.0192& 0& \cdots & 0.0400& 0.0142& 0.0026\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0.0010& \cdots & 0& 0& 0\\ 0& 0& 0& \cdots & 0& 0& 0\\ 0.0030& 0& 0& \cdots & 0& 0& 0\end{array}\right]\end{array}}$$

(5)

(2)

Calculation of fault mode propagation influence degree

According to the fault mode propagation capability matrix, the propagation influence matrix between the fault modes can be further calculated. Corresponding to the fault mode-component hypergraph, which can represent the fault propagation situation between the fault modes. The fault mode propagation influence matrix represented by the matrix $\mathit{A}$ is shown in Equation (6):

$$\mathit{A}=\mathit{B}{\mathit{B}}^T-\mathit{I}$$

(6)

where ${\mathit{B}}^T$ is the transpose of matrix $\mathit{B}$ and $\mathit{I}$ is the diagonal matrix of $\mathit{B}{\mathit{B}}^T$. Since the propagation of faults occurs between different modes, it is not necessary to consider the diagonal elements in $\mathit{A}$, and set it to 0. The elements in the matrix can be represented as

$$A=\{\begin{array}{l}0, \hspace{1em} i=j\\ {\displaystyle \sum _{k=1}^n\left(B\left[i,k\right]\cdot B\left[k,j\right]\right),i\ne j}\end{array}$$

(7)

This matrix characterizes the fault propagation between fault modes within the same component. At this time, the fault is only propagated inside the component. The off-diagonal elements represent the influence of fault propagation between modes $v_i$ and $v_j$, that is, the probability of causing another fault mode when one of the fault modes occurs. If the value of the off-diagonal element is 0, it means that there is no fault propagation between the fault modes.

We bring the matrix $\mathit{B}$ into Equation (7) to obtain the matrix $\mathit{A}$, as shown below. The one-step fault propagation trajectory of each fault mode can be obtained from the propagation influence matrix $\mathit{A}$. Taking mode 102 as an example, its one-step propagation trajectory is $v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}$, $v_{102}\stackrel{e_V}{\to }{v}_{108}$, and $v_{102}\stackrel{e_V}{\to }{v}_{409}$.

$$A=\begin{array}{c} \hspace{1em}101 \hspace{1em}301 \hspace{1em}103 \hspace{1em}\cdots \hspace{1em}203 \hspace{1em}404 \hspace{1em}506\\ \begin{array}{cc}\begin{array}{c}101\\ 301\\ 103\\ ⋮\\ 203\\ 404\\ 506\end{array}& \left[\begin{array}{ccccccc}0& 0.0311& 0.1700& \cdots & 0.0004& 0& 0.0334\\ 0.0311& 0& 0& \cdots & 0& 0.0006& 0\\ 0.1700& 0.0192& 0& \cdots & 0& 0.0002& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0.0004& 0& 0& \cdots & 0& 0& 0\\ 0& 0.0006& 0.0002& \cdots & 0& 0& 0\\ 0.0334& 0& 0& \cdots & 0& 0& 0\end{array}\right]\end{array}\end{array}\times {10}^{-2}$$

(8)

(3)

Calculation of fault mode diffusion influence degree

The fault mode propagation influence matrix $\mathit{A}$ indicates the situation of fault propagation via “shared” hyperedges. In fact, the fault not only propagates once but also further diffuses. As shown in Figure 6, the relationship between nodes and edges is $v_1,v_2\in e_1$, $v_2,v_3\in e_2$. $v_1$ can propagate faults to $v_2$, $v_2$ can propagate faults to $v_3$, and $v_1$ cannot directly propagate faults to $v_3$, but $v_2$ also belongs to the hyperedges $e_1$ and $e_2$. $v_1$ fault can still transmit the fault to $v_3$ via the “springboard” effect of $v_2$. This situation is known as fault diffusion.

It should be noted that the fault probability of the k-step reachable node is a conditional probability based on the fault of the (k − 1) step reachable node.

$$\mathit{P}={\displaystyle \sum _{r=1}^{n,r\ne i,j}\left({\mathit{A}}^{\left(k-1\right)}\left[i,r\right]\cdot \mathit{A}\left[r,j\right]\right)}$$

(9)

In the process of fault propagation, the reliability of the component will also affect the diffusion probability of the fault mode. The higher the reliability of the component, the stronger the ability to resist the fault. Therefore, combining the conditional probability with the fault probability, the diffusion influence matrix is obtained as follows:

$${\mathit{A}}^{\left(k\right)}\left[i,j\right]=\{\begin{array}{l}{\mathit{A}}^{\left(k\right)}\left[i,j\right]=0, i=j\\ {\displaystyle \sum _{r=1}^{n,r\ne i,j}\left({\mathit{A}}^{\left(k-1\right)}\left[i,r\right]\cdot \mathit{A}\left[r,j\right]\right), i\ne j,{\mathit{A}}^{\left(k-1\right)}\left[i,j\right]\ne 0}\\ {\displaystyle \sum _{r=1}^{n,r\ne i,j}\left\{{\mathit{A}}^{\left(k-1\right)}\left[i,r\right]\cdot \mathit{A}\left[r,j\right]\cdot \left(1-R_{e_j}\left(t\right)\right)\right\}, i\ne j,{\mathit{A}}^{\left(k-1\right)}\left[i,j\right]=0}\end{array}$$

(10)

According to Equations (1) and (10), combined with Table 4, the two-step fault diffusion influence matrix ${\mathit{A}}^{\left(2\right)}$ can be obtained, as shown below.

$$A^{(2)}=\begin{array}{c} \hspace{1em}101 \hspace{1em}301 \hspace{1em}103 \hspace{1em}\cdots \hspace{1em}203 \hspace{1em}404 \hspace{1em}506\\ \begin{array}{cc}\begin{array}{c}101\\ 301\\ 103\\ ⋮\\ 203\\ 404\\ 506\end{array}& \left[\begin{array}{ccccccc}0& 0.1736& 0.2688& \cdots & 0.0017& 0.0005& 0.0397\\ 0.1736& 0& 0.0542& \cdots & 0.0001& 0.0002& 0.0104\\ 0.2688& 0.0542& 0& \cdots & 0.0008& 0.0006& 0.0567\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0.0017& 0.0001& 0.0007& \cdots & 0& 0& 0.0002\\ 0.0004& 0.0002& 0.0006& \cdots & 0& 0& 0\\ 0.0397& 0.0090& 0.0503& \cdots & 0.0002& 0& 0\end{array}\right]\end{array}\end{array}\times {10}^{-5}$$

(11)

We still take mode 102 as an example, whose two-step diffusion trajectory is shown in

Table 5. Two-step diffusion trajectory of mode 102.

Diffusion Trajectory	Diffusion Trajectory
$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_k}{\to }{v}_{101}$	$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{\to }{v}_{605}$
$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_K}{\to }{v}_{105}$	$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{\to }{v}_{503}$
$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_K}{\to }{v}_{205}$	$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{⇢}{v}_{201}$
$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{\to }{v}_{101}$	$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{⇢}{v}_{502}$
$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{\to }{v}_{202}$	$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{⇢}{v}_{506}$
$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{\to }{v}_{205}$	$v_{102}\stackrel{e_{\mathit{NC}}}{\to }{v}_{501}\stackrel{e_M}{⇢}{v}_{616}$

. It is important to note that the fault will not propagate continuously. When the influence value is less than a certain value, the fault will not propagate. According to literature research [29], this paper assumes that the fault will not diffuse when the influence degree of fault diffusion is less than 10⁻⁸.

In theory, the two-step diffusion trajectory of mode 102 can realize fault diffusion, but via calculation, it is found that the two-step diffusion influence value of mode 102 to some nodes is less than 10⁻⁸, as shown in the dotted line node in Table 5. At this time, the fault will not be transmitted. Therefore, the two-step diffusion trajectory of mode 102 has a total of 8 paths reachable.

Meanwhile, according to Table 5, it can be seen that there may be multiple diffusion paths between the two modes. For example, mode 102 can pass $e_{\mathit{NC}}$ and then propagate the fault to the mode 101 through $e_K$, or it can pass $e_{\mathit{NC}}$ and then propagate to mode 101 through $e_M$.

Similarly, the three-step fault diffusion influence matrix ${\mathit{A}}^{\left(3\right)}$ of each fault mode can be calculated, as shown below.

$$A^{(3)}=\begin{array}{c} \hspace{1em}101 \hspace{1em}301 \hspace{1em}103 \hspace{1em}\cdots \hspace{1em}203 \hspace{1em}404 \hspace{1em}506\\ \begin{array}{cc}\begin{array}{c}101\\ 301\\ 103\\ ⋮\\ 203\\ 404\\ 506\end{array}& \left[\begin{array}{ccccccc}0& 0.0263& 0.0537& \cdots & 0.0002& 0.0002& 0.0026\\ 0.1095& 0& 0.0422& \cdots & 0.0001& 0& 0.0071\\ 0.4542& 0.0111& 0& \cdots & 0.0004& 0& 0.0158\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0.0014& 0.0001& 0.0004& \cdots & 0& 0& 0.0001\\ 0.0002& 0.0001& 0.0001& \cdots & 0& 0& 0\\ 0.0939& 0.0064& 0.0156& \cdots & 0.0001& 0& 0\end{array}\right]\end{array}\end{array}\times {10}^{-7}$$

(12)

From the matrix ${\mathit{A}}^{\left(3\right)}$, it can be found that when the fault diffuses in the third step, there have been many modes with the fault diffusion influence degree less than 10⁻⁸. After the calculation of the fault diffusion in the fourth step, the fault diffusion influence degree of all modes is less than the threshold of 10⁻⁸. Therefore, the fault of the machining center occurs in three steps at most at 4000 h.

As the number of diffusion steps of fault mode increases, the diffusion influence of fault is constantly decreasing, but the influence of the comprehensive fault mode is on the rise. Therefore, although the influence of diffusion is becoming smaller and smaller, its influence on the whole fault propagation cannot be ignored.

5.3. Comprehensive Influence Degree of Fault Mode

It is assumed that after the occurrence of fault mode $v_i$, k-step propagation and diffusion are conducted, and its influence on other nodes in the hypernetwork is expressed by the comprehensive influence degree of the fault mode:

$$U\left(i,k\right)=u\left(i,1\right)+u\left(i,2\right)+\cdots +u\left(i,k\right)$$

(13)

where, $1\le i\le g,1\le k\le h$. When k = 1, it means that the fault has one-step propagation, and only needs to calculate $U\left(i,1\right)$. At this time, the direct propagation influence of the fault can be obtained according to the influence degree matrix of the direct propagation of the fault. The calculation equation is as follows:

$$U\left(i,k\right)=u\left(i,1\right)={\displaystyle \sum _{j=1}^g\mathit{A}\left[i,j\right]} \left(i\ne j\right)$$

(14)

When $k>1$, it means that the fault has diffused further. It should be noted that when calculating $U\left(i,k\right)$, it is not necessary to re-calculate the influence of the node that has failed in the step k.

According to the multi-step propagation diffusion influence matrix, combined with the hypergraph model, the propagation diffusion influence value of each step of the fault mode can be obtained. Furthermore, the comprehensive influence value of the fault mode can be calculated by Equation (13), which is summarized in

Table 6. Comprehensive influence degree value at 4000 h.

Mode	Comprehensive Influence Degree	Mode	Comprehensive Influence Degree
101	0.01763	601	0.00307
301	0.00156	607	0.00010
103	0.00425	610	0.00012
102	0.00056	616	0.00054
605	0.00224	803	0.00301
503	0.00437	409	0.00008
202	0.00221	606	0.00003
501	0.00129	609	0.00003
410	0.00017	105	0.00034
603	0.00307	108	0.00008
302	0.00086	201	0.00054
403	0.00009	203	0.00003
205	0.00089	404	0.00002
805	0.00045	506	0.00054
502	0.00352

.

As can be found from Table 6, the modes with high comprehensive influence degree at 4000 h are 101, 503, 103, 502, 601, 603, and 803. However, the traditional methods for determining the influence of fault modes are mostly based on the subjective evaluation of humans, lacking objective quantitative indicators.

6. Risk Evaluation

The main factors affecting the fault risk of components are the fault rate of the component, the frequency ratio of fault modes, and the influence degree of fault modes. According to Equation (15), the risk of components can be calculated as

$$R_j\left(t\right)={\displaystyle \sum _{i=1}^g{\alpha }_{\mathit{ji}}\left(t\right){\lambda }_j\left(t\right)}{U}_j\left(i,k\right)$$

(15)

where $R_j\left(t\right)$ is the value of risk for the component $j$. $t$ is a time variable. $g$ is the number of types of fault modes. $a_{\mathit{ji}}\left(t\right)$ is the fault mode frequency ratio of component $j$ in mode $i$. ${\lambda }_j\left(t\right)$ represents the fault rate of component $j$, and $U\left(i,k\right)$ is the $i$th fault mode comprehensive influence degree.

According to Equation (15), combined with the analysis of the calculation results in Table 5 and Table 6, the risk values of the components at 4000 h are calculated, as shown in

Table 7. Risk value of machining center component (10⁻⁴).

Code	${\mathit{R}}_{\mathit{j}}$	Code	${\mathit{R}}_{\mathit{j}}$
V	0.00052	G	0.00053
W	0.01329	J	0.04204
K	0.01487	M	0.04416
D	0.00252	S	0.00388
L	0.00604	T	0.00966
NC	0.00075

.

It can be seen from Table 7 that when the machining center runs for 4000 h, the risk values of components are in descending order: $M>J>K>W>T>L>S>D>NC>$$G>V$. The component with the highest risk of fault is M, then J, and the least is V. Therefore, attention should be paid to the maintenance and health management of the tool magazine system and XYZ axis feed system.

The risk evaluation method for components of machining centers proposed in this paper is a time-varying model, which can obtain the risk value corresponding to components at any time, and fully consider the influence of propagation and diffusion of fault modes on components. In order to show its superiority, 300 h and 1000 h are selected as examples for application promotion. Considering the influence of fault propagation and diffusion, the fault risk value of the component is calculated. The comparison chart of the sorting trend of component risk value at different moments is obtained, as shown in Figure 7.

As can be found in Figure 7, the sorting results of the fault risk of the components at different times are different. The sorting result of the fault risk value at 300 h is $NC>D>L>K>W>T>J>M>V>S>G$, 1000 h is $J>M>W>K>T>L>D>S>$ $NC>V>G$ and 4000 h is $M>J>K>W>T>L>S>D>NC>G>V$. From this, it can be seen that in the early stages of machining center operation, the components with a higher risk of fault are the numerical control system and hydraulic system. As time goes by, there is a greater risk of faults in the tool magazine and feed system. This is because, with the extension of the operation time of the machining center, the comprehensive influence degree of the tool magazine and the feed system is increasing. This is in line with the engineering.

7. Conclusions

In this paper, we presented a framework for the evaluation of the fault risk. Fault mode analysis is conducted to explore the fault mechanism and obtain the fault mode frequency ratio of components. The fault rate function is built based on Weibull distribution. A hypergraph model is established to calculate the comprehensive influence degree of fault mode. The risk evaluation is carried out based on the fault mode frequency ratio, fault rate, and the comprehensive influence of fault. A real fault data set of the CNC machining center has been used to present the procedure of risk evaluation. The main concluding remarks are given as follows:

(1)

The most common fault mode of this series of machining centers is 101, 103, and 301. Therefore, it is of great significance to ensure the reliable quality of components for the reliable operation of the machining center.

(2)

In the process of fault propagation, although the diffusion influence is becoming smaller and smaller, its cumulative influence on the whole fault propagation cannot be ignored. The modes with a high comprehensive influence degree at 4000 h are 101, 503, 103, 502, 601, 603, and 803.

(3)

Based on the proposed fault risk evaluation method, the risk value of the machining center at 4000 h is obtained, and the sorting result is $M>J>K>W>T>L>S>D>\mathit{NC}>G>V$. The component with the highest risk of fault is M, then J, and the least one is V. Therefore, attention should be paid to the maintenance and health management of the tool magazine system and XYZ axis feed system.

(4)

Via the calculation, we found that the sorting of the risk value of the machining center components is different at different times. At 300 h is $NC>D>L>K>W>T>J>M>V>S>G$, 1000 h is $J>M>W>K>T>L>D>S>NC>V>G$, and 4000 h is $M>J>K>W>T>L>S>D>NC>G>V$. The influence of fault is not static, but dynamic.

However, this paper still has some shortcomings and needs to be improved in the future. (1) The research on fault propagation and diffusion in this paper is based on statistical data, which can be integrated into real-time online monitoring data to obtain more accurate risk evaluation results. (2) This paper does not consider the impact of different maintenance strategies on the risk of component fault in machining centers. However, different levels of maintenance will inevitably have different impacts on component faults. Therefore, in future research, a more comprehensive risk evaluation framework should be established considering the impact of maintenance strategies.