Physics-Informed Neural Networks for the Condition Monitoring of Rotating Shafts

1. Introduction

Rotating shafts are elements of engineering systems that play a paramount role in the transmission of power, encompassing speed and torque, from one point to another [1], [2]. They are typically designed to endure substantial loads and to operate at high velocities, underscoring the need for precise alignment, equilibrium, and freedom from imperfections. These considerations are significant not only for enhancing the overall system performance, but also for improving its safety and reliability [3], [4]. To attain this objective, the practice of condition monitoring (CM) for rotating shafts allows continuously evaluating the shaft condition and performance and detecting any indications of malfunction or deterioration [5], [6]. Through the application of CM methodologies, potential issues can be promptly identified, and maintenance or repair actions can be driven to prevent accidents. This approach has shown to effectively mitigate the risk of unexpected downtime, reinforcing the overall system reliability [7], [8]. Furthermore, the widespread and cost-effective availability of sensors has revolutionized the acquisition of diagnostic signals, such as accelerations, strains, and elastic waves [9]. However, the availability of big data is itself a new layer of complexity, especially in the realm of signal processing. That is, the rapid expansion of the amount of acquired data has introduced the need for (i) improved hardware and software performance, and (ii) developing tools to deal with confounding factors, including those unrelated to the system health state, such as environmental and operational conditions [10], [11].

To tackle these challenges, deep learning has stood as a pivotal technological advancement in the CM of rotating machines, offering multifaced contributions of significant importance [12]. This approach has excelled in automatically extracting intricate patterns and features from raw sensor data, enhancing the precision and reliability of fault detection and anomaly characterization. As an example, the work in [13] applied deep learning to enhance wind turbine CM, addressing the data surge from increased wind farm units. By combining convolutional neural networks (CNNs) [14] and recurrent neural networks (RNNs) [15], it efficiently extracted features, reduced dimensionality, and provided effective CM, offering both real-time unit state checks and early warning capability, even amidst accidental parameter changes. In [16], a CM model based on CNNs for automatic fault detection in rotating equipment was developed. The model, utilizing data from a single vibration sensor on the motor-drive end bearing, achieved accuracies of $99.58\%$ and $97.3\%$ when applied to two different databases under controlled ambient conditions. Another example was presented in [17], where the authors proposed a novel deep learning algorithm for detecting rotor unbalance in industrial machinery. The algorithm, extracting important vibration signatures such as fast Fourier transform (FFT) and short-time Fourier transform (STFT), combined the depth of ResNet [18] and the feature extraction capability of CNN. This hybrid approach surpassed the performance of both individual models. The study involved two analyses: binary detection of balanced vs. unbalanced cases and multilevel detection of the degree of unbalance.

However, so far, most of the methods used in this field have been based on black-box like deep learning algorithms. That is, users do not have any clues about how the input data are processed through the network for getting the predictions [19]. Moreover, it is unclear if the network behaviour adheres to the physics of the problem. Two approaches have been proposed in the literature to solve this limitation: using explainability algorithms or adopting more complex physics-based deep learning algorithms. While the former strategy is used after training to debunk the behaviour of the network, it does not bring evidence that physical laws are satisfied [20], [21], [22]. Instead, the latter solution allows making sure that the algorithm makes predictions that adhere to the physics of the problem, and typically involves including regularization terms, i.e., known physical laws, to constrain the optimization problem during training. Moreover, physics-informed algorithms stand out as flexible tools suitable for various applications: data-driven solution of partial differential equations, data-driven discovery of physical laws and parameters estimation [23], [24], [25]. In the context of CM, only a few contributions have been proposed involving the integration of the physical knowledge in the training process of deep learning models. The work in [26] presented a physics-informed deep learning method for bearing fault detection that combined a threshold model and a CNN. The approach was validated using data from bearings on an agricultural machine and a laboratory test stand in the Case Western Reserve University Bearing Data Centre. In [27], a method for identifying unbalance faults in rotary systems using physics-guided neural networks (PGNNs) was proposed. The approach involved the use of a standard neural network to localize the nodal position of the experimental fault, followed by PGNN to quantify the unbalance magnitude and phase angle. Instead, the work in [28] introduced a novel physics-informed convolution long-short-term memory (LSTM-CNN) network for rotor unbalance and shaft cracks detection and localization. In particular, the physics were taken into account through the construction of a neural network model which mimicked a finite element (FE) resolution of the problem.

To the best of the authors' knowledge, still no efforts have been made for the direct estimation of multiple parameters characterizing the health state of a rotating shaft system by leveraging physics-informed neural networks (PINNs). In this work, PINNs are utilized to estimate critical health state parameters in a simple but realistic numerical case of an extended Jeffcott rotor model. This model incorporates damping effects and anisotropic supports for a more comprehensive representation. The parameters under consideration include the radial and angular position of the static unbalance caused by the disk on the shaft, stiffness along the principal axes of elasticity, and the non-rotating damping coefficient. The estimation is exclusively based on the displacement signals from the disk centre. Note that this estimation not only optimizes the performance of machineries, enhancing efficiency and reliability, but also enables predictive maintenance by identifying potential faults early on. To highlight the effectiveness and precision of the proposed methodology, various scenarios with different constant rotational speeds are examined. Furthermore, the analysis accounts for the impact of noisy input data.

The paper is organised as follows: Section 2 offers a brief overview of the necessary theoretical foundations about PINNs for parameter estimation; Section 3 shortly presents the case study and then shows in detail the implementation and the results of the PINN for the system health state characterization. Finally, Section 4 provides some concluding remarks.

2. Methodology

The proposed framework hinges upon the use of PINNs to estimate the unknown parameters characterizing the dynamics of a rotating shaft system. Notably, PINNs represent deep learning tools that combine NNs with the system governing equations, and are particularly useful when data might be limited or noisy, and where the underlying physics of the problem is well understood [23].

Assume that a generic physical system is governed by the $n$-th order ordinary differential equations (ODEs) shown in Eq. (1):

$${\mathit{u}}^{\left(n\right)}=\mathit{F}(t,\mathit{u},{\mathit{u}}^{\prime },{\mathit{u}}^{\prime \prime },\dots ,{\mathit{u}}^{\left(n-1\right)},\mathit{\lambda })$$

(1)

where $t$ refers to the system independent variable, $\mathit{u}=[u_1\left(t\right),u_2\left(t\right),\dots ,u_p\left(t\right)]$ denotes the state vector consisting of $p$ components defined in the domain $[t_0,t_f]$, and $\mathit{\lambda }=[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k]$ represents the vector made of the $k$ unknown parameters describing the system state. Subsequently, considering the NN universal approximation theorem [29], a NN can be exploited to obtain an approximation $\widehat{\mathit{u}}=\mathit{N}(\mathit{W},\mathit{b},\mathit{\lambda })$ of the state vector $\mathit{u}$, such that $\widehat{\mathit{u}}\approx \mathit{u}$. More specifically, $\mathit{W}$ and $\mathit{b}$ denote the weight and bias matrices of the NN, respectively, and their values are the result of a training process [30], as well as for the parameter vector $\mathit{\lambda }$. Note that, since $\widehat{\mathit{u}}$ is a function, its derivatives concerning the independent variable $t$ can be computed during the training process through automatic differentiation (AD) [31], [32]. Then, a function $\mathit{g}$ outlining the approximation of Eq. (1) can be defined, as reported in Eq. (2):

$$\mathit{g}={\widehat{\mathit{u}}}^{\left(n\right)}-\mathit{F}(t,\widehat{\mathit{u}},{\widehat{\mathit{u}}}^{\prime },{\widehat{\mathit{u}}}^{\prime \prime },\dots ,{\widehat{\mathit{u}}}^{\left(n-1\right)},\mathit{\lambda })$$

(2)

To enable the neural network to fine-tune the parameters $\mathit{W},\mathit{b}$ and $\mathit{\lambda }$ in order to (i) fulfil the underlying ODEs describing the system behaviour and (ii) to fit the available data (i.e., gathered measurements, in which the state vector $\mathit{u}$ is known), two different loss functions are considered, as shown in the following Eq. (3) and Eq. (4).

$${\mathcal{L}}_f=\frac{1}{N_e}·\sum _{i=1}^{N_e}{\left|\mathit{g}\left(t_e^i\right)\right|}^2$$

(3)

$${\mathcal{L}}_u=\frac{1}{N_e}·\sum _{i=1}^{N_e}{\left|\mathit{u}\left(t_e^i\right)-\widehat{\mathit{u}}\left(t_e^i\right)\right|}^2$$

(4)

where ${\mathcal{L}}_f$ denotes the loss for the ODEs fulfilling while ${\mathcal{L}}_u$ is the loss of the observed data; $t_e^i$ indicates the generic $i$-th element of the vector ${\mathit{t}}_e$, made of $N_f$ elements inside the domain $[t_0,t_f]$, in which ${\mathcal{L}}_f$ and ${\mathcal{L}}_u$ are evaluated. These loss functions are subsequently integrated to yield the loss term $\mathcal{L}$, as presented in Eq. (5):

$$\mathcal{L}{=\alpha ·\mathcal{L}}_f+{\beta ·\mathcal{L}}_u$$

(5)

where $\alpha$ and $\beta$ denote two coefficients employed to assign greater weight either to the contributions derived from the accessible data or those related to the system physics. Consequently, the objective of minimizing $\mathcal{L}$ is enforced, enabling the PINN to infer the unidentified parameters that define the system dynamics. A scheme showing how the PINN is trained is presented in Figure 1.

3. Case study

3.1. Extended Jeffcott rotor with unknown system parameters

The PINN approach described in the previous Section is here applied on a numerical case study of a rotating shaft system in which the dynamics is simulated with the use of an extended Jeffcott rotor model [33]. The system, illustrated in Figure 2, consists of a $1m$ long rotating shaft (i.e., $l=1m$) made of aluminium (Young’s modulus $E$ set to $70000MPa$, and density $\rho$ equal to $2700\frac{kg}{m^3}$) and supported at both ends. It has a circular cross section with a diameter of $20mm$, and it incorporates a disk (representing, for instance, a flywheel, fan, turbine, gear, etc.) that is mounted at distances $l_1$ and $l_2$ from the respective supports. The shaft rotates at a velocity $\Omega =\frac{d\vartheta \left(t\right)}{dt}$, in which $\vartheta \left(t\right)$ denotes the angle defined with respect to the $x$ axis of the right-handed $xyz$ reference frame which is fixed in space. In this reference frame, the disk lies within the $xy$ plane, and the $z$ axis is aligned with the line connecting the two supports. Without any loss of generality, the disk is here considered to be positioned at the midpoint of the shaft, i.e., $l_1=l_2=l/2$. Moreover, the disk centre of mass $P$ is displaced from the axis of rotation, whose trace in the disk plane is identified with the point $C$, generating a static unbalance defined by the distance $\epsilon$ and the angle $\phi =\vartheta \left(0\right)$. This unbalance causes the point $C$ to displace from the line joining the supports leading the shaft to whirl around it.

The numerical model used to compute the rotor dynamics consists of a system of two second-order ODEs with a state vector of two components (making reference to Section 2, $n=2$ and $p=2$, respectively), as reported in Eq. (6) and Eq. (7):

$$\left[\begin{array}{cc}m& 0\\ 0& m\end{array}\right]\left\{\begin{array}{c}{\ddot{x}}_c\\ {\ddot{y}}_c\end{array}\right\}+\left[\begin{array}{cc}{c}_n+c_r& 0\\ 0& c_n+c_r\end{array}\right]\left\{\begin{array}{c}{\dot{x}}_c\\ {\dot{y}}_c\end{array}\right\}+\left[\begin{array}{cc}{k}_x& {\Omega ·c}_r\\ -\Omega {·c}_r& k_y\end{array}\right]\left\{\begin{array}{c}{x}_c\\ y_c\end{array}\right\}=\left\{\begin{array}{c}{f}_x\\ f_y\end{array}\right\}$$

(6)

$$\left\{\begin{array}{c}{f}_x\\ f_y\end{array}\right\}=\left\{\begin{array}{c}{m}_d·\epsilon ·\left({\Omega }^2·cos\left(\Omega ·t+\phi \right)+\dot{\Omega }·sin\left(\Omega ·t+\phi \right)\right)\\ m_d·\epsilon ·\left({\Omega }^2·sin\left(\Omega ·t+\phi \right)-\dot{\Omega }·cos\left(\Omega ·t+\phi \right)\right)-m·g\end{array}\right\}$$

(7)

where $x_c$ and $y_c$ denote the $x$ and $y$ coordinates of the point $C$, respectively, and represent the components of the state vector (i.e., $\mathit{u}=[x_c,y_c]$); $c_n$ and $c_r$ are the equivalent viscous damping terms of the stationary and rotating parts of the system, respectively; $m$ represents the summation of the disk mass $m_d$ and the shaft equivalent mass at the disk location $m_s^e$; $g$ indicates the gravity acceleration; $k_x$ and $k_y$ are the system stiffnesses along the $x$ and $y$ axes, respectively, that are assumed to coincide with the axes of the ellipse of elasticity, i.e., the principal axes of elasticity of the supporting structure. Notably, the rotating damping $c_r$ is here considered to be coincident to the contribution given by the shaft material properties, thus neglecting any potential additional term, and it is computed exploiting the approximation of a linear system [33], i.e., $c_r=2{\xi }_r\sqrt{k_s·m_s}$. Here, ${\xi }_r$ denotes the rotational damping coefficient, that is assumed to be $0.001$, $m_s$ represents the mass of the shaft, while $k_s=\frac{48·E·I}{l^3}$ denotes the shaft flexural stiffness, with $I$ indicating the moment of inertia of the shaft cross section. The disk mass $m_d$ is considered to be equal to $2kg$, while the shaft equivalent mass at the disk location $m_s^e$ is computed as $m_s^e=\frac{k_s·\delta }{g}$, in which $\delta =\frac{5·m_s·g·l^3}{384·E·I}$ denotes the static displacement of the shaft at the disk location due to the shaft weight. The stiffnesses $k_x$ and $k_y$ are the combination of the shaft stiffness $k_s$ and those of the supports $k_b^x$ and $k_b^y$ along the $x$ and $y$ axes, respectively. That is, $\frac{1}{k_x}=\frac{1}{k_s}+\frac{1}{k_b^x}$ and $\frac{1}{k_y}=\frac{1}{k_s}+\frac{1}{k_b^y}$. However, determining the stiffness values of the supports can be challenging for different reasons, e.g., when they have complex geometries and interactions, due to misalignments or imperfections in installation, when the support materials is not well-defined or uniform, and they can change due to degradation over time [34]. A similar reasoning applies for the non-rotating damping and for the static unbalance. Hence, in this work $\epsilon$, $\phi$, $c_n$, $k_x$ and $k_y$ represent the components of the parameter vector $\mathit{\lambda }$ , and their values are estimated through a PINN. Such parameters are selected because tracking their value is essential for maintaining the performance, reliability, and safety of rotating machinery.

Figure 3 shows the model responses $x_c$ and $y_c$ obtained by solving Eq. (6) and Eq. (7) with a Runge-Kutta 4/5 integration method [35] in the time range $\left[\mathrm{0,10}\right]s$. In the simulated scenario a constant rotational speed of $\Omega =30\frac{\mathrm{r}\mathrm{a}\mathrm{d}}{\mathrm{s}}$ is imposed, and $k_x=7.76\frac{N}{mm}$, $k_y=6.71\frac{N}{mm}$, $c_n=7.0·{10}^{-3}\frac{N·s}{mm}$, $\epsilon =8mm$, $\phi =10deg$. No artificial noise is added. Figure 4 shows the scatter plot of the position of the disk centre $C$ in the $xy$ plane over time.

3.2. System state characterization through physics-informed neural networks

A PINN is then exploited to estimate the unknown parameters of the analysed rotating shaft system, i.e., to estimate the parameter vector $\mathit{\lambda }=[\epsilon ,\phi ,c_n,k_x,k_y]$. Notably, the employed NN architecture consists of 1 input neuron that takes in the generic time instant $t$, 1 hidden layer made of 200 neurons, and 2 output neurons to predict the values of $x_C$ and $y_C$. Moreover, the hyperbolic tangent $tanh$ activation function [36] is used in the hidden layer, while the output layer embeds a linear activation function. The true system responses $x_C$ and $y_C$ required for training the PINN are numerically obtained with the Runge-Kutta 4/5 integration method in the time range $\left[\mathrm{0,1}\right]s$ (i.e., $t_0=0s$ and $t_f=1s$). Within this range, the true solution is sampled with a sampling frequency of $10kHz$, which means that 10001 equally spaced points are considered in time. The same time instants are also considered for building the vector ${\mathit{t}}_e$ that is used for training (i.e., $N_f=10001$). It is worth noting that various sampling frequency values were examined in this study. Specifically, the investigation covered a sampling frequency range of $[1,10]kHz$ with increments of $3kHz$. Due to the similarity in outcomes from the analysis, these results are omitted here for the sake of conciseness.

A representative scenario is examined to assess the efficacy of the proposed methodology. The scenario involves the imposition of a constant rotational speed of $\Omega =45\frac{\mathrm{r}\mathrm{a}\mathrm{d}}{\mathrm{s}}$, alongside the following unknown parameters: $\epsilon =8.00mm$, $\phi =10.00deg$, $c_n=7.00·{10}^{-3}\frac{N·s}{mm}$, $k_x=5.53\frac{N}{mm}$, $k_y=7.25\frac{N}{mm}$. No artificial noise is added to the system responses. The PINN is trained on an AMD Ryzen 9 5900HX 3.30 $GHz$ processor using a limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) optimization algorithm [37] with a learning rate of $\eta =0.01$, and encompassing $10000$ training iterations. Moreover, the arbitrary initial guess set for the unknown parameters are $\epsilon =0.00mm$, $\phi =0.00deg$, $c_n=0.001·{10}^{-3}\frac{N·s}{mm}$, $k_x=1.00\frac{N}{mm}$, $k_y=2.00\frac{N}{mm}$. The training process, as depicted in Figure 5, illustrates the progressive reduction of the training loss as the number of iterations increases. Instead, Figure 6 shows how the trained PINN is able to fit the available data to estimate the unknown parameters, showing a comparison between the true solution and the one obtained with the PINN. The outcome reveals that the approximation of the state vector offered by the PINN closely aligns with the actual vector, showcasing marginal disparities primarily observed in the state variable $y_C$ during the initial time instants. Finally, the true value and the correspondent PINN estimation for all the system unknown parameters is shown in Table 1. What emerges is that all the parameters are estimated by the PINN with a remarkable accuracy, presenting the best performance in identifying the system stiffness values $k_x$ and $k_y$, where the relative error in the estimation remains below $0.70\%$. Notably, among the parameters under estimation, the angle $\phi$ of the static unbalance and the non-rotating damping $c_n$ prove to be the most challenging. This observation finds potential justification in the relatively subdued impact these parameters exert on the state variables compared to their counterparts. Nevertheless, even in these instances, the relative error remains constrained within $5.80\%$.

In order to assess the robustness of the PINN to external influences affecting input data (e.g., measurement noise, environmental vibrations, etc.), the same scenario is revisited while maintaining consistent training process conditions (i.e., algorithm used, learning rate, initial parameters guess, etc.). In this context, the input variables $x_C$ and $y_C$, employed for the estimation of system parameters, are subjected to a perturbation through the introduction of supplementary numerical noise with a signal-to-noise ratio ($SNR$) of $40dB$. The data fitting performed by the PINN and the parameters estimation are reported in Figure 7 and Table 2, respectively. The outcomes demonstrate that the PINN continues to function as a dependable tool for characterizing the health state of the system, even when confronted with external disturbances within the input data. Notably, no indications of diminished algorithm performance are discernible, with all parameter estimations retaining a relative error of less than $6.00\%$.

Lastly, an additional example is presented below to assess the estimation capabilities of the PINN when fed with data pertaining to a scenario marked by distinct imposed conditions and health states. The scenario encompasses the application of a constant rotational speed of $\Omega =35\frac{\mathrm{r}\mathrm{a}\mathrm{d}}{\mathrm{s}}$, accompanied by the following unknown parameters: $\epsilon =12.00mm$, $\phi =20.00deg$, $c_n=5.00·{10}^{-3}\frac{N·s}{mm}$, $k_x=7.76\frac{N}{mm}$, $k_y=6.14\frac{N}{mm}$. Moreover, an artificial noise with an $SNR$ of $30dB$ is added to the input data, and the same training conditions of the previous scenarios are kept. The true and PINN solutions of the state variables $x_C$ and $y_C$ are reported in Figure 8, and the estimation of the unknown parameters is shown in Table 3. The PINN approximation of the state variables $x_C$ and $y_C$ appears to closely match the true solution. That is, the PINN manages to smooth out all the perturbances introduced by the added numerical noise, thus acting as a filter. The unknown parameters are satisfactorily estimated even in this scenario. The lowest relative estimation error characterises the system stiffnesses $k_x$ and $k_y$, i.e., $0.52\%$ and $0.16\%$, respectively, while the non-rotating damping $c_n$ is identified with a $9.80\%$ error, which may be mainly related to the considerable amount of noise added.

4. Discussion

The paper has introduced a novel approach employing PINNs for estimating unknown parameters characterizing the health state of rotating shaft systems. The investigation has focused on a realistic numerical case study involving an extended Jeffcott rotor model, which has incorporated damping effects and anisotropic supports. The parameters considered have encompassed the radial and angular position of static unbalance caused by a shaft-mounted disk, stiffness values along the principal axes of elasticity, and the non-rotating damping coefficient. The estimation has relied exclusively on displacement signals from the disk centre, and various scenarios, incorporating different constant rotational speeds, have been thoroughly examined. Results have revealed the implemented PINN accuracy in estimating these parameters, demonstrating minimal relative errors even in the presence of substantial data noise. Beyond optimizing machinery performance and enhancing efficiency and reliability, the proposed estimation method has facilitated predictive maintenance by early fault identification.

Future works may explore (i) extending the approach to more complex rotor systems, (ii) experimental validation, and (iii) conducting comparative studies with other cutting-edge methods and machine learning approaches for parameter estimation in rotating systems.