Quantitative Fault Diagnostics of Hydraulic Cylinder Using Particle Filter

Abstract

Condition-based hydraulic cylinder maintenance necessitates quantitative fault diagnostics. However, existing methods are characterized by either qualitative or limited quantitative capabilities. In this paper, a quantitative fault diagnostic method using a particle filter for hydraulic cylinders is proposed. The problem of quantitative fault diagnostics is formally formulated in a stochastic framework to assess the health/fault state, and an architecture based on joint state-parameter estimation is proposed. Through the establishment and analysis of a nonlinear dynamic model of the hydraulic cylinder, the impact of time-varying parameters on the state variables is revealed. Three fault modes of the cylinder, including friction, internal leakage, and external leakage, are theoretically identified. The proposed method allows for a simultaneous quantitative diagnosis of these three fault modes. The performance of the proposed method is evaluated using meticulously designed experiments. The results demonstrate that the mean absolute percentage errors in the parameter estimations are below 9% (accuracy exceeding 91%), thus validating its feasibility and effectiveness.

1. Introduction

The hydraulic cylinder plays a vital role in hydraulic transmission and control systems. The health/fault state of the hydraulic cylinder directly affects the safety as well as the functioning of a hydraulic system. Common issues associated with hydraulic cylinders encompass excessive friction, internal leakage, and external leakage. These problems can be attributed to factors such as wear and tear on the reciprocating seal [1], temperature fluctuations [2], alterations in the properties of the working fluid [3], contaminants, surface texture [4], etc. According to statistics, hydraulic cylinder leakage faults make up over 50% of all crane faults [5]. To improve the availability, reliability, and safety of hydraulic systems, condition-based maintenance has emerged as a significant trend in the industry [6,7]. Fault diagnostics for hydraulic cylinders is a crucial aspect of this trend. Typically, sensors cannot directly measure internal leakages in hydraulic cylinders. Developing an effective condition monitoring method that can accurately predict conditions beyond the available measurements presents a significant challenge. Furthermore, for effective condition-based maintenance, there is typically a need for quantitative health/fault information from fault diagnostics. However, most existing methods can essentially only offer qualitative or limited quantitative diagnostics.

In the early stages, fault diagnostics for hydraulic cylinders primarily relied on knowledge-based artificial intelligence (AI) methods, such as fuzzy logic [8] and rule-based expert systems [9]. The primary drawback of this class of methods is the substantial amount of domain expertise required to formulate the rules. Furthermore, it was difficult for these rules to cover all potential scenarios comprehensively. As a result, knowledge-based AI methods waned in popularity. Over the past decade, hydraulic cylinder fault diagnostic techniques have gradually shifted towards two main categories: data-driven and model-based (also referred to as physics-based) methods.

Data-driven methods rely on labeled historical data that include training samples of both normal and faulty system behavior without depending on a mathematical model of the system. Once the data are collected, feature extraction techniques are applied to convert the raw data into a set of informative features. These features capture relevant information, patterns, or trends in the data that may indicate faults, making them crucial and drawing significant attention. Tan et al. [10] used experimental data to study fault diagnostics by extracting various features from water hydraulic cylinder vibration signals. Zhao et al. [11] explored energy fluctuations within hydraulic cylinders at different leakage levels using wavelet packet analysis to extract features. Tang et al. [12] established a mapping relationship between hydraulic cylinder leakage and energy entropy obtained from the wavelet packet transform. In addition to wavelet transform, other methods such as Fourier transform, short-time Fourier transform, Wigner-Ville distribution, Hilbert transform, and more are also employed for feature extraction [13]. However, it is worth noting that the wavelet transform is the most widely utilized method in the case of hydraulic cylinders. After feature extraction, data-driven methods often involve supervised machine learning (ML) classification algorithms. Typically, cylinder fault models and levels are categorized into discrete classes. For example, internal leakage is usually classified as none (representing a healthy state), small, medium, and large. Machine learning classification algorithms are then employed to map features to these predefined classes. Due to the increasing popularity of ML, these methods have gained significant traction in recent years. For instance, El-Betar et al. [14] used a feedforward neural network to detect internal leakage in a hydraulic cylinder. Jin et al. [5] proposed diagnosing hydraulic cylinder seal wear and internal leakage using wavelet transform and wavelet neural networks. Wang et al. [15] extracted features from pressure signals using wavelet decomposition and inputted these features into a BP neural network (BPNN) to identify internal leakage. Ma et al. [16] achieved a notably high level of accuracy in the detection of internal leakage by employing sensitive feature extraction based on wavelet packet transform in conjunction with a support vector machine. In addition to conventional ML, deep learning methods have garnered attention in the field of hydraulic cylinder fault diagnostics. For example, Wang et al. [17] proposed a combination of sliding-window spectrum feature extraction and a deep belief network (DBN) for fault diagnosis in hydraulic systems. Zhang and Chen [18] developed a fault diagnosis method based on an optimized DBN combined with complete ensemble empirical mode decomposition with adaptive noise techniques. Guo et al. [19] compared the performance of different regression algorithms in predicting internal leakage, including convolutional neural network (CNN), BPNN, radial basis function network, support vector regression, T-S neural network, and Elman neural network. The results showed that CNN outperformed other algorithms. Similarly, Na et al. [20] combined wavelet analysis-based feature extraction with CNN to detect leakage.

Data-driven methods can be effective when there is a lack of in-depth understanding of the system’s underlying physics. However, these methods are limited in their ability to cover only several discrete fault scenarios. When a fault is classified as “large”, it remains vague and unclear how large it is exactly for maintenance decision-makers. The information offered by data-driven methods may be too coarse for condition-based maintenance. Essentially, these methods are qualitative because they cannot provide a quantitative description of the fault’s severity. Moreover, their performance may suffer when encountering different operating conditions than those present in the training dataset. Since the training data cannot encompass all possible operating conditions, ML algorithms typically struggle to generalize to unforeseen data that deviates significantly from the training set. Lastly, obtaining high-quality training data in practice is quite challenging. For effective fault diagnostics, operational data that include faulty situations is required. This often necessitates a large amount of run-to-failure data, which is usually not feasible in reality. A common practice is to use simulation or experimental data with artificially injected faults. However, the extent to which these data accurately represent real situations is debatable. Even if run-to-failure data are available, labeling such a vast amount of data is also a formidable challenge.

In contrast to data-driven methods, model-based methods rely on accurate mathematical models of the system. Based on how the mathematical model is used, model-based methods can be further divided into analytical redundancy approaches and state estimator approaches. In the analytical redundancy approach, a mathematical model is meticulously developed to describe the expected normal behavior of the system. During operation, sensor data or measurements from the actual system are compared to the model’s output. Discrepancies between the measured values and those calculated by the model are referred to as residuals. Fault diagnostic algorithms analyze these residuals to identify deviations from the expected normal behavior. When a significant difference is detected, it signals the occurrence of some type of fault in the system. For instance, An and Sepehri [21] utilized an extended Kalman filter (KF) to estimate the states of the hydraulic cylinder and used the residuals to detect internal and external leakage. Sepasi and Sassani [22] enhanced this method by employing an unscented KF to better adapt to system models with higher non-linearity. Sun et al. [23] applied a series of extended KFs to diagnose multiple fault modes of an electro-hydraulic actuator. Analytical redundancy approaches are highly effective for fault detection. However, these methods are typically limited to detecting faults in the cylinder and may not be able to isolate specific faults or describe their severity. To address this limitation, state estimator-based approaches are developed. Garimella and Yao [24] developed an adaptive, robust observer capable of identifying several common faults that can occur in hydraulic cylinder drive units. These faults include inadequate supply pressure, decreased hydraulic compliance, and excessive hydraulic fluid leakage. Tan and Sepehri [25] proposed a versatile order-recursive estimation scheme to estimate the parameters of the cylinder model. These efforts can be considered early attempts at quantitative fault diagnostics for hydraulic cylinders. However, it is important to note that these methods are deterministic in nature. In reality, numerous unavoidable sources of uncertainty impact state estimation. These sources include uncertainties in the initial state, sensor noise, process noise, disturbance inputs, and more. To make informed maintenance decisions, decision-makers need to assess the quality of the estimation and then decide whether to accept or reject the estimation results. Therefore, for state estimation-based methods, providing probability distributions is preferred over a set of deterministic values. Furthermore, even though quantitative fault diagnostics has been partially utilized in existing studies, there is still a gap in formally formulating this type of problem.

To overcome these limitations, this paper aims to tackle the problem of quantitative fault diagnostics for hydraulic cylinders within a stochastic framework, thereby facilitating effective decision-making in the context of condition-based maintenance. The main new contributions of this study are as follows:

(1)

A physics-based quantitative fault diagnostics problem is formally formulated, and a joint state-parameter estimation-based architecture is proposed.

(2)

The fault modes and their impact on the hydraulic cylinder are analytically identified and revealed through the establishment of a nonlinear dynamic model.

(3)

A particle filter-based quantitative fault diagnostics method is proposed and validated, enabling accurate quantitative diagnosis of multiple faults for hydraulic cylinders.

The remainder of this paper is organized as follows. Section 2 presents the formulation and architecture of the physics-based quantitative fault diagnostics problem. A non-linear dynamic model of the hydraulic cylinder is developed, validated, and analyzed in Section 3. Section 4 is dedicated to addressing particle filter-based joint state-parameter estimation. In Section 5, the performance of the proposed method is evaluated and analyzed. Comparison with existing methods is discussed in Section 6. Concluding remarks are drawn in Section 7.

2. Quantitative Fault Diagnostics Approach

This section formally formulates the problem of the physics-based quantitative fault diagnostics approach and presents its architecture.

2.1. Problem Formulation

Assume the component or system can be described by a discrete-time model as follows:

$${\mathit{x}}_{k+1}=f\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{u}}_k,{\mathit{Q}}_k\right)$$

(1)

$${\mathit{z}}_k=h\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{u}}_k,{\mathit{R}}_k\right)$$

(2)

where $k\in \mathbb{N}$ is the discrete time variable, $\mathit{x}\in {ℝ}^{n_x}$ is the state vector, $\mathit{\theta }\in {ℝ}^{n_{\theta }}$ is the time-varying parameter vector, $\mathit{u}\in {ℝ}^{n_u}$ is the input vector, $\mathit{Q}\in {ℝ}^{n_Q}$ is the process noise, $\mathit{z}\in {ℝ}^{n_z}$ is the noisy measurement, $\mathit{R}\in {ℝ}^{n_R}$ is the measurement noise, $f\left(·\right)$ is the state transition function, and $h\left(·\right)$ is the measurement function. This representation represents a general nonlinear model without imposing any restrictions on the functional forms of $f\left(·\right)$ or $h\left(·\right)$. Moreover, the noise terms could exhibit nonlinear coupling with both the states and parameters. The parameters ${\mathit{\theta }}_k$ evolve in an unknown fashion, and although they are commonly assumed to remain constant in practical scenarios, they could potentially undergo time-varying alterations.

A fault can be characterized as a situation in which a component or system fails to meet a designated set of functional requirements. These requirements may be defined by a threshold, beyond which we determine that the component or system has malfunctioned. Typically, this threshold can be formulated as a function dependent on the system’s state and parameters:

$$F\left({\mathit{x}}_k,{\mathit{\theta }}_k\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if} \mathrm{a} \mathrm{requirement} \mathrm{is} \mathrm{violated}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(3)

The function $F\left(·\right)$ determines whether a fault has occurred in the component or system. It yields a value of 1 if a requirement is violated, and 0 otherwise.

The degree of fault/health, on the other hand, typically depends on the time-varying parameters. It can be quantitatively described by measuring the distance between the current parameter value and its maximum allowable value. Define vector ${\mathit{\lambda }}^{+}$ as the reciprocal of the maximum allowable value for each parameter, as illustrated in Equation (4).

$${\mathit{\lambda }}^{+}={\left[\frac{1}{{\theta }_{1,\max }},\frac{1}{{\theta }_{2,\max }},\dots ,\frac{1}{{\theta }_{n_{\theta },\max }}\right]}^T$$

(4)

Then, we can establish the degree of fault as follows:

$$DF\left({\mathit{\theta }}_k\right)={\mathit{\lambda }}^{+}\odot {\mathit{\theta }}_k$$

(5)

where $\odot$ denotes component-wise multiplication, and the degree of health can be expressed as follows:

$$DH\left({\mathit{\theta }}_k\right)={\mathit{I}}_{ones}-DF\left({\mathit{\theta }}_k\right)$$

(6)

where ${\mathit{I}}_{ones}$ is a vector of ones with the length of $n_{\theta }$.

Based on the definitions of fault and degree of fault as provided in Equations (3) and (5), the quantitative fault diagnostics problem can be transformed into a joint state-parameter estimation problem. In practical scenarios, the estimation process is subject to numerous sources of uncertainty. Uncertainty in the initial state, model, process noise, and measurement noise results in an uncertain estimation of (x_k, θ_k). Under these circumstances, to give a meaningful estimation, it is important to obtain a probability distribution of the state-parameter estimation as opposed to a single point of estimation. Thus, the goal of quantitative fault diagnostics is to estimate the probability of the state and parameter given the input-output data, i.e., $p\left(\left.{\mathit{x}}_k,{\mathit{\theta }}_k\right|{\mathit{u}}_{0:k-1},{\mathit{z}}_{0:k-1}\right)$.

2.2. Quantitative Fault Diagnostics Architecture

The proposed quantitative fault diagnostics architecture comprises a system model, a joint state-parameter estimator, and the degree of fault function, as depicted in Figure 1. The system model and the joint state-parameter estimator form the core of this architecture.

The input signal sequence ‘u’ is used as the input to the physical plant, and the noisy output measurement ‘z’ is generated as the output. When the same input signal is applied to the system model, an estimate of the output $\hat{\mathit{z}}$ can be obtained. With the output measurement, estimated output, and information about measurement noise taken into consideration, the probability distribution of the states and time-varying parameters can be calculated by the joint state-parameter estimator. The probability distribution of the states and time-varying parameters is then passed to the fault degree function $DF\left(·\right)$, where a quantitative description of the degree of fault can be obtained.

The available joint state-parameter estimators are outlined in

Table 1. Available estimators.

Estimator	Model	Assumed Distribution
Kalman filter	Linear	Gaussian
Extend Kalman filter	Locally linear	Gaussian
Unscented Kalman filter	Non-linear	Gaussian
Particle filter	Non-linear	Non-Gaussian

. Specifically, the KF family assumes that the state/parameter variables and measurements follow linear Gaussian dynamics. The standard KF is well-suited for linear systems, while the extended KF can accommodate locally linear systems. In cases where the system displays high non-linearity, the unscented KF can be employed, although it still relies on the Gaussian distribution assumption. Conversely, the particle filter can be applied to any non-linear system without necessitating assumptions about the underlying distribution. KFs are computationally efficient for linear systems and provide closed-form solutions for linear Gaussian problems. However, they may encounter difficulties with highly non-linear systems. Particle filters can handle non-linear systems but can become computationally intensive as the number of particles increases. They are more suitable for situations where other methods, like the KF, may fail. The choice between the KF and the particle filter depends on the specific characteristics of the system. If the system adheres to a linear Gaussian model, the KF is efficient and provides accurate estimates. However, if the system exhibits highly non-linear or non-Gaussian characteristics, the particle filter offers greater flexibility.

3. Modeling and Analysis of the System

Accurate state-parameter estimation necessitates a highly accurate system model. In this section, we initially develop a nonlinear dynamic model for the hydraulic cylinder. Subsequently, we validate this model by comparing its output results with those of a Simulink^® SimScape model. By utilizing the system model, we analyze the impact of time-varying parameters, enabling us to theoretically identify the fault modes of the hydraulic cylinder.

3.1. Modeling of the System

The schematic of the system is depicted in Figure 2. The force of the cylinder is regulated by a fixed pressure source at the supply, and the direction of the piston motion is controlled by a directional control valve (DCV). For the purpose of facilitating analysis, the following assumptions are adopted:

(1)

The temperature, viscosity, and bulk modulus of the working fluid remain constant.

(2)

Leakage flows are characterized as laminar.

(3)

The pressure within the working chambers of the cylinder is uniformly distributed.

Figure 2. Schematic diagram of the system.

Figure 2. Schematic diagram of the system.

Under these assumptions, we begin by examining the extension movement, with the DCV positioned to the left. The equation governing pressure buildup in the piston chamber of the cylinder can be expressed as follows:

$${\dot{P}}_1=\frac{\beta }{A_{1}x+V_{pis,dead}}\left(Q_{1,ext}-Q_{il}-A_1\dot{x}\right)$$

(7)

The pressure buildup equation for the rod chamber of the cylinder can be expressed as follows:

$${\dot{P}}_2=\frac{\beta }{A_2\left(L_s-x\right)+V_{rod,dead}}\left(Q_{il}-Q_{el}-Q_{2,ext}+A_2\dot{x}\right)$$

(8)

The flow equation for the inlet port of the piston chamber in extension motion is as follows:

$$Q_{1,ext}=C_{d1}{a}_1\sqrt{\frac{2\left|P_s-P_1\right|}{\rho }}sign\left(P_s-P_1\right)$$

(9)

The flow equation for the outlet port of the rod chamber in extension motion is as follows:

$$Q_{2,ext}=C_{d2}{a}_2\sqrt{\frac{2\left|P_2-P_b\right|}{\rho }}sign\left(P_2-P_b\right)$$

(10)

For the retraction movement, the DCV is switched to the right position, and the pressure build-up equation for the piston chamber can be expressed as follows:

$${\dot{P}}_1=\frac{\beta }{A_{1}x+V_{pis,dead}}\left(-Q_{1,ret}-Q_{il}-A_1\dot{x}\right)$$

(11)

The pressure build-up equation for the rod chamber can be expressed as follows:

$${\dot{P}}_2=\frac{\beta }{A_2\left(L_s-x\right)+V_{rod,dead}}\left(Q_{il}-Q_{el}+Q_{2,ret}+A_2\dot{x}\right)$$

(12)

The flow equation for the inlet port of the piston chamber in retraction motion is as follows:

$$Q_{1,ret}=C_{d1}{a}_1\sqrt{\frac{2\left|P_1-P_b\right|}{\rho }}sign\left(P_1-P_b\right)$$

(13)

The flow equation for the outlet port of the rod chamber in retraction motion is as follows:

$$Q_{2,ret}=C_{d2}{a}_2\sqrt{\frac{2\left|P_s-P_2\right|}{\rho }}sign\left(P_s-P_2\right)$$

(14)

For both the extension and retraction cases, the motion equation of the piston is as follows:

$$P_1{A}_1-P_2{A}_2-B\dot{x}-Kx-F_l+F_c=m\ddot{x}$$

(15)

The internal leakage can be expressed as follows:

$$Q_{il}=C_{il}\left(P_1-P_2\right)$$

(16)

The external leakage can be expressed as follows:

$$Q_{el}=C_{el}{P}_2$$

(17)

And the contact force can be expressed as follows:

$$F_c=\left\{\begin{array}{ll}-k_{c}x-B_c\dot{x},\hfill & \mathrm{if} x<0\hfill \\ 0,\hfill & \mathrm{if} 0\le x\le L_s\hfill \\ -k_c\left(x-L_s\right)-B_c\dot{x},\hfill & \mathrm{if} x>L_s\hfill \end{array}\right.$$

(18)

Define the stats vector of the system as follows:

$$\mathit{x}={\left[x_1,x_2,x_3,x_4\right]}^T={\left[x,\dot{x},P_1,P_2\right]}^T$$

(19)

Define the input vector as follows:

$$\mathit{u}={\left[u_1,u_2,u_3\right]}^T={\left[P_s,P_b,F_l\right]}^T$$

(20)

Define the time-varying parameter vector as follows:

$$\mathit{\theta }={\left[B,C_{il},C_{el}\right]}^T$$

(21)

Based on the system model given in Equations (7)–(18), the state transition model of the system can be expressed as follows:

$$\dot{\mathit{x}}=f\left(\mathit{x},\mathit{u},\mathit{\theta }\right)$$

(22)

The system output can be expressed as follows:

$$\mathit{z}=h\left(\mathit{x},\mathit{R}\right)=\mathit{x}+\mathit{R}$$

(23)

To numerically solve the non-linear continuous-time state transition model, the 4th-order Runge–Kutta method is adopted in this study. The 4th-order Runge–Kutta method can be expressed as follows:

$$\left\{\begin{array}{c}{K}_1=f\left({\mathit{x}}_k,{\mathit{u}}_k,{\mathit{\theta }}_k\right)\hfill \\ K_2=f\left({\mathit{x}}_k+\frac{T_s}{2}{K}_1,{\mathit{u}}_k,{\mathit{\theta }}_k\right)\hfill \\ K_3=f\left({\mathit{x}}_k+\frac{T_s}{2}{K}_2,{\mathit{u}}_k,{\mathit{\theta }}_k\right)\hfill \\ K_4=f\left({\mathit{x}}_k+T_s{K}_3,{\mathit{u}}_k,{\mathit{\theta }}_k\right)\hfill \\ {\mathit{x}}_{k+1}={\mathit{x}}_k+\frac{T_s}{6}\left(K_1+2K_2+2K_3+K_4\right)\hfill \end{array}\right.$$

(24)

And the discrete output measurement is as follows:

$${\mathit{z}}_k={\mathit{x}}_k+{\mathit{R}}_k$$

(25)

3.2. Verification of the Model

To validate the analytical model presented in the previous paragraphs, a Simulink^® SimScape model of the corresponding system was constructed. By employing the identical time-invariant parameters outlined in

Table 2. Main time-invariant parameters of the model.

No.	Time-Invariant Parameters	Values	No.	Time-Invariant Parameters	Values
1	m	500 kg	8	a1	7.854 × 10−5 m2
2	K	10 N/m	9	a2	7.854 × 10−5 m2
3	A1	1.76 × 10−2 m2	10	ρ	870 kg/m3
4	A2	8.45 × 10−3 m2	11	Ls	0.8 m
5	β	7 × 108 Pa	12	kc	1 × 108 N/m
6	Cd1	0.7	13	Bc	1.5 × 105 N/(m/s)
7	Cd2	0.7

, along with the designated input signals, a comparison was made between the output results of the analytical model and the SimScape model. This comparison is presented in Figure 3. The evident agreement between the analytical model and the SimScape model serves to validate the efficacy of the analytical model. The small differences are given by the numerical integrator choice.

3.3. Impact of the Time-Varying Parameters

Upon confirming the efficacy of the analytical model, it becomes applicable for analyzing the impact of time-varying parameters on the system states. As illustrated in Figure 4a, an increase in the viscous coefficient B leads to a reduction in both extension and retraction speeds. Consequently, the time required to complete a full stroke increases. This is highly comprehensible since the viscous coefficient B is directly associated with the friction force. During the extension stroke, P₁ functions as the working pressure, and P₂ acts as the back pressure. Conversely, in the retraction stroke, P₁ switches roles to become the back pressure, while P₂ becomes the working pressure. Throughout both strokes, the friction force consistently opposes the direction of motion. As a result, the impact of B on P₁ and P₂ demonstrates a contrasting pattern. During the extension stroke, an increase in B leads to a rise in P₁ and a decrease in P₂. Conversely, in the retraction stroke, an increase in B causes P₁ to decrease and P₂ to increase. In magnitude, variations in B have a notably substantial impact on x, v, and P₁, while the effect on P₂ is comparatively minor.

The internal leakage coefficient C_il is directly connected to the rate of internal leakage flow via Equation (16). This influence of C_il is depicted in Figure 4b. The effect of C_il on both x and v closely resembles that of B. An increase in C_il corresponds to a decrease in both extension and retraction speeds, consequently leading to an extended duration for extension/retraction. With an increase in C_il, both P₁ and P₂ exhibit slight increments during the extension stroke, while they decrease in the retraction stroke. Notably, the influence of C_il on P₂ is slightly more pronounced than on P₁. The variation in C_il has a more prominent effect on P₁ and P₂ when the piston is at rest compared to when it is in motion.

The external leakage coefficient C_el corresponds to the rate of external leakage flow, as shown in Equation (17). Referring to Figure 4c, an increase in C_el has a minor effect on the extension speed, but significantly reduces the retraction speed. This disparity in impact can be attributed to the occurrence of external leakage exclusively at the piston rod seal. Consequently, the influence on the pressure in the rod chamber (P₂) is notably greater than that on the pressure in the piston chamber (P₁).

It is evident that the time-varying parameters exert substantial influences on the system’s state variables, leading to diverse system performances. If the performance no longer aligns with the specified requirements, it can be inferred that a fault has occurred. With reference to the three time-varying parameters, three distinct fault modes can be identified for hydraulic cylinders: friction fault, internal leakage fault, and external leakage fault.

Note that various factors, including temperature, viscosity, bulk modulus, and the presence of contaminants within the working fluid, among others, can also significantly influence the performance of a hydraulic cylinder and potentially lead to malfunctions. It is important to recognize that these parameters are characteristics of the working fluid itself and not inherent properties of the cylinder. Consequently, in this study, we treat these fluid-related parameters as constant values over time.

Friction, internal leakage, and external leakage in a hydraulic system collectively impact both system functionality and the environment. Excessive friction within the system can lead to reduced efficiency, heat generation, component wear, and decreased precision, all of which can affect its overall performance. Internal and external leakage reduce efficiency through the loss of hydraulic fluid. Meanwhile, external leakage presents significant environmental concerns, as it can result in the release of hazardous hydraulic fluid into the surroundings, posing threats to ecosystems and human safety.

The most common causes of friction-related faults include, but are not limited to, worn seals or guides, the presence of contaminants, side-loading, and component damage. Within a hydraulic cylinder, seals, including piston and rod seals, can generate friction when they come into contact with the cylinder walls and moving components, such as the piston rod. If these seals become worn or damaged, they have the potential to increase friction levels. Particles or contaminants in the hydraulic fluid can infiltrate the space between the seals and the cylinder walls, leading to heightened friction. Additionally, improper alignment of the hydraulic cylinder or uneven distribution of the load can result in uneven pressure on the cylinder walls. In some extreme cases, this can cause the cylinder rod to bend, ultimately resulting in increased friction. Lastly, the wear and tear of components like piston rods and cylinder walls, as illustrated in Figure 5a [26], can also contribute to an increase in friction.

Both internal and external leakage can result from a variety of factors, including damaged or worn seals, contaminants, cylinder and rod damage, excessive pressure, and more. Among the primary contributors to both types of leakage are deteriorated or worn seals within the hydraulic cylinder. Seals can deteriorate over time due to factors like friction, pressure, contamination, and extreme temperatures. Physical damage to the hydraulic cylinder or rod, such as nicks, scratches, bends, or corrosion, can affect the rod seal and cause fluid to leak both internally and externally. Operating the hydraulic system at pressures exceeding its design limits can induce seal failure, contributing to both internal and external leakage. The phenomenon of hydraulic cylinder leakage is depicted in Figure 5b,c [27].

4. Joint State-Parameter Estimation

In this study, we utilize a particle filter to conduct joint state-parameter estimation. This choice is made because the hydraulic cylinder model exhibits non-linearity, and there is no evidence supporting the assumption of a Gaussian distribution of the state parameter. The particle filter method approximates the posterior distribution through Monte Carlo sampling. Specifically, a particle filter utilizes a set of independent random realizations, referred to as particles, which are directly sampled from the state space. These particles serve as representations of the posterior probability and are updated as new observations become available. The Bayesian formula is employed to appropriately position, weight, and recursively propagate the particles.

To realize joint state-parameter estimation, the state vector is augmented with the time-varying parameter vector as follows:

$$\tilde{\mathit{x}}={\left[\mathit{x},\mathit{\theta }\right]}^T$$

(26)

In this way, the state and time-varying parameters can be estimated simultaneously. Since the time-varying parameters for the hydraulic cylinder typically exhibit slow growth in practice, it is reasonable to assume that their time derivatives are zero in the state transition model.

In a particle filter, the augmented state distribution at time step k is approximated by a set of weighted particles, as follows:

$${\left\{w_k^i,{\tilde{\mathit{x}}}_k^i\right\}}_{i=1}^N$$

(27)

where the superscript i is the particle index, subscript k is the time step, N denotes the number of particles, ${\tilde{\mathit{x}}}_k^i$ is the augmented state particles, and $w_k^i$ is the non-negative numerical factors named importance weights, which sum to one. Each particle ${\tilde{\mathit{x}}}_k^i$ represents a possible realization of the augmented state. The importance weight $w_k^i$ represents the relative importance of each of the N particles.

In this study, the sequential importance sampling with resampling (SISR) particle filter is utilized. The pseudocode for a single time step of the SISR particle filter is given in Algorithm 1, and the corresponding flow chart is depicted in Figure 6. The algorithm takes in a set of weighted particles ${\left\{w_{k-1}^i,{\tilde{\mathit{x}}}_{k-1}^i\right\}}_{i=1}^N$ that represents the previous belief of the augmented states, as well as inputs from control action ${\mathit{u}}_{k-1}$ and measurement ${\mathit{z}}_k$. It then calculates and provides the updated state distribution ${\left\{w_k^i,{\tilde{\mathit{x}}}_k^i\right\}}_{i=1}^N$. Each particle is propagated forward to time step k by using the state transition model $f\left(·\right)$ established in Section 3.1. The propagated particles are then evolved by adding assumed Gaussian process noise. The importance weight is updated based on the measurement ${\mathit{z}}_k$. Specifically, the measurement model $h\left(·\right)$ is applied to the propagated particles to calculate the estimated output ${\hat{\mathit{z}}}_k$, and the likelihood of their corresponding outputs is calculated using the assumed probability density function of the measurement noise, as shown in Equation (28).

Algorithm 1: Sequential importance sampling with resampling particle filter
Input: ${\left\{w_{k-1}^i,{\tilde{\mathit{x}}}_{k-1}^i\right\}}_{i=1}^N$, ${\mathit{u}}_{k-1}$, ${\mathit{z}}_k$
Output: ${\left\{w_k^i,{\tilde{\mathit{x}}}_k^i\right\}}_{i=1}^N$
for i = 1 to N do
	${\tilde{\mathit{x}}}_k^i\sim p\left({\tilde{\mathit{x}}}_k\left|{\tilde{\mathit{x}}}_{k-1},{\mathit{u}}_{k-1}\right.\right)$# Propagate particles
	${\tilde{\mathit{x}}}_k^i={\tilde{\mathit{x}}}_k^i+{\mathit{Q}}_k^i$# Particles evolution with process noise
	$w_k^i=p\left({\mathit{z}}_k\left|{\tilde{\mathit{x}}}_k^i\right.\right)$# Update importance weight
end
$w_{sum}={\displaystyle \sum _{i=1}^N{w}_k^i}$# Cumulative weight
for i = 1 to N do
	$w_k^i=w_k^i/w_{sum}$# Normalize weights
end
if $N_{eff}<N_{threshold}$ then       # Resample if needed
	Resample N particles with replacement
	for i = 1 to N do
		$w_k^i=1/N$# Reset weights
	end
end

$$p\left({\mathit{z}}_k\left|{\tilde{\mathit{x}}}_k^i\right.\right)=\prod N\left({\mathit{z}}_k\left|{\hat{\mathit{z}}}_k,\Sigma \right.\right)=\prod \frac{\exp \left(-\frac{1}{2}{\left({\mathit{z}}_k-{\hat{\mathit{z}}}_k\right)}^T{\Sigma }^{-1}\left({\mathit{z}}_k-{\hat{\mathit{z}}}_k\right)\right)}{\sqrt{{\left(2\pi \right)}^{n_m}\left|\Sigma \right|}}$$

(28)

where n_m is the number of the measurement, equal to 4 in our case, and $\Sigma$ is the covariance matrix. The mean vector of this distribution is set to be equal to the estimated output ${\hat{\mathit{z}}}_k$, and the covariance matrix equals the covariance of the measurement noise. We evaluate this distribution at the actual measurement ${\mathit{z}}_k$. Then, we multiply the resulting values together to obtain a single likelihood value for the particle. Once the weights are updated for each particle, a normalization step is performed. This normalization involves dividing each particle’s weight by the cumulative sum of all the weights, ensuring that they collectively add up to one.

A particle filter typically begins with initialized particles, all assigned equal weights. However, only a few particles may be positioned close to the true values of the augmented state. As the algorithm progresses, any particle not aligning with the measurements will receive an extremely low weight, while those particles in proximity to the true augmented state values will carry a significant weight. When the filter encounters this scenario, it becomes degenerate. Typically, this problem is resolved by using some form of particle resampling. To minimize computational expenses, particle degeneracy is evaluated by calculating the effective sample size, as represented in Equation (29):

$$N_{eff}=\frac{1}{{\displaystyle \sum _{i=1}^N{\left(w_k^i\right)}^2}}$$

(29)

Resampling is performed whenever the effective sample size falls below a threshold, N_threshold. Following resampling, the particle weights are set to equal values of 1/N.

After updating the weights and particles using Algorithm 1, one can calculate the estimated augmented states at time step k by simply summing the weighted values of the particles:

$${\mathit{\mu }}_k=\frac{1}{N}{\displaystyle \sum _{i=1}^N{w}_k^i{\tilde{\mathit{x}}}_k^i}$$

(30)

The confidence interval (CI) can be determined by calculating the nth percentile of the particles. For instance, the lower boundary of the 95% CI corresponds to the 2.5th percentile of the particles, whereas the upper boundary corresponds to the 97.5th percentile of the particles.

For a deeper explanation of the theoretical aspects of the particle filter, readers can refer to references [28,29].

5. Performance Evaluation and Analysis

To evaluate the performance of the proposed method, we designed four experiments, as outlined in

Table 3. Design of experiments.

Experiments	Operating Condition 1	Parameters True Values 1
Exp. 1	Ps = 25, Pb = 0, FL,ext = 20, FL,ret = 0	B = 10, Cil = 20 × 10−11, Cel = 10 × 10−10
Exp. 2	Ps = 25, Pb = 0, FL,ext = 20, FL,ret = 0	B = 20, Cil = 40 × 10−11, Cel = 15 × 10−10
Exp. 3	Ps = 35, Pb = 0, FL,ext = 30, FL,ret = 0	B = 10, Cil = 20 × 10−11, Cel = 10 × 10−10
Exp. 4	Ps = 15, Pb = 0, FL,ext = 10, FL,ret = 0	B = 15, Cil = 25 × 10−11, Cel = 20 × 10−10

¹ Units: P_s, P_b: bar; F_L,ext, F_L,ret: kN; B: kN/(m/s); C_il, C_el: m³/(s·Pa).

. Specifically, Exp. 1 serves as a baseline. In comparison to Exp. 1, Exp. 2 involves altering the true values of the parameters to assess the method’s performance when dealing with variations in parameter values while keeping the operating conditions the same. In contrast, Exp. 3 maintains the same parameter values as Exp. 1 but changes the operating conditions to evaluate the generalization ability of the proposed method under different operational scenarios. Exp. 4 is used to further test the robustness of the proposed method, where both the operating condition and parameter values are changed with respect to the first three experiments. Note that varying parameter values correspond to different degrees of fault severity.

In all four experiments, the true values and measurements are obtained from the SimScape simulation, and the number of particles is consistently set at 500. The initialization of the particles in experiment 1 is depicted in Figure 7. Specifically, since the states can be measured using sensors, we initialize the state particles using a Gaussian distribution based on the mean and standard deviation values given in

Table 4. Particle initialization setup.

State [Unit]	[Mean, Standard Deviation]	Parameter [Unit]	[Minimum, Maximum]
x [m]	[0, 0.0018]	B [kN/(m/s)]	[0, 70]
v [m/s]	[0, 0.0006]	Cil [10−11 m3/(s·Pa)]	[0, 50]
P1 [bar]	[0, 0.03]	Cel [10−10 m3/(s·Pa)]	[0, 25]
P2 [bar]	[0, 0.03]

. Conversely, for the time-varying parameters, no additional information is available aside from their minimum and maximum values. Therefore, the time-varying parameter particles are initialized using a uniform distribution with the corresponding minimum and maximum values given in Table 4. The same initialization method was employed for the other experiments.

The evolution of the particles for the time-varying parameters of experiment 1 is depicted in Figure 8 at intervals of 400 time steps. It is important to recognize that the dispersion of the particles reflects the level of uncertainty. At time step 0, the particles are uniformly distributed within the respective minimum and maximum values of the time-varying parameters, signifying the highest level of uncertainty. As the iterations progress, the particles concentrate together, indicating a significant reduction in uncertainty. Meanwhile, Figure 8 serves as evidence of the particle filter algorithm’s convergence. One can observe different convergence rates among the parameters. Notably, the viscous coefficient B converges most rapidly, whereas the external leakage coefficient C_el converges more slowly, and the convergence speed of the internal leakage coefficient C_il falls in between.

The state estimation results and their corresponding 95% confidence intervals (CIs) for Exp. 1 are depicted in Figure 9. As evident from the observation, the particle filter state estimations demonstrate a high level of accuracy, closely matching the corresponding true values. Note that in Figure 9, with the exception of the initial time step, the 95% CIs are conspicuously narrow and closely aligned with the estimations, indicating a minimal level of uncertainty associated with the estimations. This result can also be supported by Figure 8, where it is evident that the particles representing the possible parameter values cluster closely together. It can also be observed that the state estimation values exhibit smaller variations in comparison to the measurement values, resulting in more accurate state values than those provided by the sensors.

Compared to state estimation, our primary focus is on the estimation of time-varying parameters. This emphasis is due to the fact that the time-varying parameter values are directly related to the degree of health of the cylinder. The time-varying parameter estimation results and the corresponding 95% CIs for Exp. 1 are illustrated in Figure 10. To quantitatively assess the accuracy of the estimations, we employ the performance metrics known as the mean absolute error (MAE) and mean absolute percentage error (MAPE), as respectively defined in Equations (31) and (32). In these equations, n represents the number of data points, y stands for the true value, and $\hat{y}$ denotes the estimated value. These values are also depicted in Figure 10 for reference. After convergence, the particle filter estimation results closely match the true values, offering a reliable estimate of the time-varying parameter values under the operating conditions of Exp. 1. The 95% CIs are notably tight, ensuring high-confidence estimations.

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|y-\hat{y}\right|}$$

(31)

$$MAPE=\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|\frac{y-\hat{y}}{y}\right|}\right)\times 100\%$$

(32)

The time-varying parameter estimation results for Exps. 2 and 3 are illustrated in Figure 11 and Figure 12, respectively. The results of Exp. 4 are demonstrated in

Table 5. Parameter estimation results of Exp. 4.

	B	Cil	Cel
MAE	0.057	0.149	0.020
MAPE	0.382%	0.595%	0.100%

. It is clear that the particle filter is capable of providing accurate parameter estimations even when faced with varied parameter values and different operating conditions. By employing parameter value estimations, Equation (5) can be applied to determine the degree of the fault, thus enabling a quantitative fault diagnostic of the hydraulic cylinder.

6. Discussion

It is demonstrated that the proposed quantitative fault diagnostic method using a particle filter is capable of providing accurate fault information beyond sensor measurements. A summary of the comparison with existing methods is presented in

Table 6. Fault diagnostics method comparison.

Methods	Qualitative/Quantitative	Deterministic/Stochastic	Generalization Ability
Data-driven	Qualitative	Deterministic/Stochastic	Poor
Analytical redundancy	Qualitative	Deterministic	Good
Adaptive robust observer	Quantitative	Deterministic	Good
Particle filter	Quantitative	Stochastic	Good

. This comparison is assessed along three dimensions: qualitative/quantitative, deterministic/stochastic, and generalization ability.

Data-driven methods based on supervised machine learning prove effective in situations where a physical model is unavailable, but comprehensive data encompassing both healthy and faulty conditions is accessible. They are predominantly qualitative, as they can solely address several predefined discrete fault severities. Depending on the specific algorithm employed, data-driven methods can be either deterministic or stochastic. However, their generalization ability is rather limited, primarily due to their dependence on training data. Analytical redundancy methods are physics-based. They excel in fault detection, but they often face challenges in diagnosing faults from the residuals. On the other hand, the adaptive robust observer method offers quantitative fault diagnostics but lacks the capability to provide confidence intervals for parameter estimations. In contrast to these approaches, the proposed method effectively addresses and overcomes these limitations.

7. Conclusions

In this study, a quantitative fault diagnostics method utilizing a particle filter for hydraulic cylinders is proposed and validated. The problem of physics-based quantitative fault diagnostics is formally formulated in a stochastic framework, and a joint state-parameter estimation-based architecture is proposed. An analysis of time-varying parameters reveals that variations in the friction coefficient significantly impact the piston position, speed, and piston chamber pressure, with minor effects on rod chamber pressure. Internal leakage coefficient variations have a more pronounced effect on piston and rod chamber pressures when the piston is stationary compared to when it is in motion. Additionally, external leakage has a notably greater influence on the pressure in the rod chamber compared to the pressure in the piston chamber.

The proposed method is evaluated via meticulously designed experiments. It demonstrates that the proposed method reliably provides accurate state-parameter estimations with accuracy exceeding 91%, even when facing varying parameter values and different operational conditions. This capability enables a quantitative assessment of the fault state. Notably, the proposed method ensures tight 95% confidence intervals after convergence. In addition, the proposed method allows for multiple fault diagnoses. These robust characteristics render the proposed method highly suitable for condition-based maintenance of hydraulic cylinders. Future work will be focused on the improvement of the computational cost of the proposed methodology and further experimental validations.