Establishing a Real-Time Multi-Step Ahead Forecasting Model of Unbalance Fault in a Rotor-Bearing System

Abstract

Recently, prognostics and health management (PHM) has garnered a lot of attention in the industrial sector for its cost-effective maintenance and safe operation of the system. In this regard, vibration-based predictive maintenance using sensors plays a significant role in the diagnosis and prognosis of various faults. The need of the hour is to know when and which part must be replaced in advance for efficient and reliable operation. Unbalance is one major fault acting on any rotary system leading to excessive vibration and causing various other faults developing early failure in components directly or indirectly. In this paper, we show how a prognostic model can be built for the identification of future unbalance trend of a rotor-bearing system with the aid of a mathematical model of the system and statistical/machine learning methods. The prognostic model developed is used to forecast the unbalance time-series data of an industrial turbine rotor in real-time which forecasts the month ahead unbalance values. The proposed model is verified for prognostic analysis using datasets from a local plastic company. After careful examination of the results, it is concluded that the proposed model can aid in precisely detecting future system unbalance.

1. Introduction

Rotating machines are the crucial framework for any industry. However, numerous factors, such as varying loads or fatigue, high speeds, etc., can lead to sudden failure that may further result in shutdowns or downtime of several weeks in an industry. For economic and safety reasons, it is always necessary to ensure the safe working condition of rotary systems. Prognostics and health management (PHM) is an emerging field in this regard that deals with predictive maintenance to reduce the probability of a sudden or extreme system failure and thereby improving safety standards of any industry. Lee et al. [1] analyzed a detailed review of a systematic approach for establishing PHM based maintenance strategy for Rotating systems. Tahan et al. [2] also gave a reviewed analysis of PHM based diagnostics and prognostics for gas turbines with the help of condition-based monitoring. Generally, vibration plays a vital role in determining various types of faults such as unbalances, misalignments, cracks, etc., in a rotary system, that eventually causes a component to fail. Thus, studying vibration measurements can aid us in diagnosing the faults in the system much earlier than it takes place. Lv et al. [3] reviewed signal processing methods based on vibration for early fault diagnosis and prognosis in rotating machinery. Out of the different faults that occur in a rotary system [4], unbalance contributes the most to its failure. There is always some inherent unbalance present in a rotary system during operation due to reasons such as manufacturing defects, installation faults, etc., which ends up causing undesirable effects on the safe working of the system. Unbalance is one of the primary faults which eventually cascades into many other secondary faults. Thus, it becomes essential to continuously monitor the unbalance faults present in rotating machinery and take maintenance actions at regular intervals. In this regard, condition based monitoring (CBM) refers to continuously monitoring the condition of the system while it is in operation and formulating preventive strategies based on that [5]. However, prognostics applied to a fault like unbalance has been exceedingly difficult due to less availability of failure data and thus very little research has been done so far in this area. Therefore, it becomes essential to continuously monitor the unbalance of the system and if it gets sufficiently high, then the system needs to be balanced.

Diagnosis generally deals with the identification and localization of faults in the system [6,7,8,9] whereas prognosis refers to the development of the fault trend with time to understand its developing nature and making predictive analysis based on that. Arunthavananthan et al. [10] supplied a roadmap study of fault diagnostic methodology and process safety research for Industry 4.0. Wei et al. [11] provided a guide on the fault diagnosis approach for rotary machinery systems. Nath et al. [12], Lei et al. [13], and Liu et al. [14] contributed a detailed review of the artificial intelligence (AI) based rotor fault diagnosis (RFD) framework and described its future research direction. The focus of this analysis was to bridge the gap between real-world industrial solutions and Laboratory based solutions using machine learning and deep learning methodologies. In the past, researchers have developed different methodologies for unbalance estimation and diagnosis.

The main aim of prognostics is to predict the failure or remaining useful life (RUL) of a component which would enable more reliable operation, reduced downtime, shortened maintenance cost, and hence would provide more effective maintenance of the overall system. Rezaeianjouybari and Shang [15] discussed several ongoing research on deep learning for the development of prognostics framework. Kan et al. [16] described a detailed review of prognostics techniques on rotating machinery under non-linear and non-stationary conditions. Sikorska et al. [17] did an elaborate analysis of steps for the determination of RUL-based prognostic approach by various existing methods such as physical model, data-driven, or knowledge-based architecture and the present difficulties in its application to a real-world system. A lot of prognostics study right now generally involves diagnosing crack propagation and failure due to fatigue [18]. However, there is need of developing prognostic solutions for various faults acting in a system.

The prognostic modeling approach can be categorized into physics-based, data-driven, or hybrid i.e., combined methodologies. CBM generally involves identifying the fault indicators using data-driven methodologies and predicting the remaining useful life of components from fault trend analysis. However, there also exists a model-based prognostic strategy, in which a mathematical model is built to determine the fault’s gradual development with time. The hybrid prognostic approach is a combination of both data-driven and model-based strategies to utilize the benefits of the two. Heng et al. [19] provided a detailed discussion of prognostic models based on condition data of rotating machinery systems. An et al. [20] and Sikorska et al. [17] presented an overview of physics and the data-driven models for prognostics. In these studies, various data-driven health monitoring systems for establishing diagnostic and prognostic systems have been described for predictive maintenance. Zhong et al. [21] described an overview of data-driven methodologies for fault prognosis in Industrial systems.

Several data-driven prognostic models have been previously utilized which use time series data for data-driven prognosis [19], however they mostly evaluated on single step ahead predictions to estimate the variations of vibration-based characteristics of fault feature. Accuracy indeed decreases with the increase of prediction horizon, but it is of a more significant nature to have a more reliable long-term prediction for any important maintenance-based decision. Multi-step ahead prognosis is still a big challenge when it comes to prognosis since the accuracy of predicted values deviates increasingly from the original values as the time step increases, as this has also been stated as an important research direction in several review papers [12]. Xiao et al. [22] utilized a back propagation neural network (BPNN) over a rolling time window for multi-step ahead prediction to determine the RUL using bearing run-to-failure simulation data and PHM-2012 competition data.

Trend analysis gets easier when historical failure data is given, but in a lot of scenarios it is not viable or reliable for the system to fail and necessary pre-conditioning or maintenance is done. This results in less availability of relevant historical data for any meaningful evaluation. It is important to predict the trend of how a fault feature develops over time. Guo et al. [23] proposed a method that combines errors from multiple sparse auto-encoders with Long Short-Term Memory (LSTM) method and produces a trend curve that will predict future mechanical fault variation curves.

Various sensors and measurements made at different bearing locations can enable us to predict the overall unbalance acting in a system. We can further analyze our measured data with the help of forecasting to predict future unbalance values which in return can save a lot of time and money for the effective running of the machinery.

Forecasting techniques prove to be extremely useful while predicting the future parameters of a system based on past data trends. In this, time series forecasting (TSF) refers to the process of making future data predictions based on the analysis of past and present data trends. There are several classical or statistical TSF methods such as Exponential Smoothing (ES), and Auto Regressive Integrated Moving Average (ARIMA) which are successfully being used for effective forecasting in many different sectors. These classical methods are generally simple and can give a good prediction with an acceptable error range. But as the complexity of the system increases the simplistic classical algorithms do not keep up with the sudden variations of the system parameters due to several unknown factors and can give significantly large errors. Hence, there is a need to explore emerging technologies such as machine learning techniques which can prove better in learning the complexities of a system and making more accurate long-term predictions.

Recently, several machine learning techniques have been employed to develop forecasting models for fault time trend analysis and prognostic study of the system. Pham et al. [24] used ARMA/GARCH for the forecasting of machine health to predict future failure states. They verified their model by applying it to a real compressor in a petrochemical plant. Yuan et al. [25] developed a support vector machine model optimized with an artificial immunization algorithm (SVM-AIA) for fault diagnosis of a turbo-pump rotor system. Wilson Wang [26] introduced an adaptive neuro-fuzzy approach-based forecasting model to predict real-time dynamic system behaviour in Industry. He applied his adaptive predictor to gear system condition monitoring and material fatigue testing and found his model to perform significantly better for a real-time system than classical forecasting methods. Li et al. [27] used an enhanced metabolism grey forecasting method (MGFM) using a particle filter (PF) for evaluating the remaining service life. They used the case study of a centrifugal pump and compressor to validate their results Djeddi et al. [28] used the LSTM model for prognostic analysis of Gas Turbine data. Zhang et. al. [29] employed a multi-layer LSTM model with an attention mechanism for the prediction of the RUL of a rotary machine. Xiao et al. [30] proposed a stacked long short-term memory with a multi-layer perceptron (SLSTM-MLP) method for multivariate time series forecasting of wind turbine bearing systems utilizing supervisory control and data acquisition (SCADA) data. Ma et al. [31] proposed a data fusion approach based on LSTM and ARIMA methodologies to forecast long-term fuel cell degradation. Arunthavanathan et al. [32] used a convolutional neural network and long short-term memory (CNN-LSTM) approach to develop a prognostic approach to forecast future system parameters to diagnose the faults within the system. Xia et al. [33] proposed a data-driven prognosis approach by using a hybrid deep LSTM network structure. In the model, they used the dropout technique to overcome the overfitting issue. They utilized the National Aeronautics and Space Administration (NASA) commercial C-MAPSS dataset with turbofan engine degradation data to evaluate the model performance and RUL prediction.

This study aims to develop a prognostic approach to determine the unbalance acting in the turbine rotor-bearing system to aid in maintenance-based decision-making. A mathematical model along with historical sensor data is used to generate fault data and statistical/machine learning methods like ES, ARIMA, SVR, LSTM, etc., are evaluated for developing a prognostic model of unbalance over a rolling or expanding window for a dataset from a local factory. The main idea driving the development of this approach is to set up an online monitoring and predictive maintenance strategy for unbalance faults in a rotary system. So, we formulate a real-time unbalance forecasting model which updates itself over time and gives satisfactory one-month-ahead forecasting with an acceptable error range. To validate our forecasting model, we use different traditional or classical as well as non-classical or machine learning based forecasting models to find out the most suitable model for our system and calculate its error in prediction with actual trend. The major contributions of this paper are as follows,

• A prognostic model approach is proposed with the help of moving and expanding window methodologies to help in predicting unbalance estimation in the future.
• The different forecasting models such as ES, ARIMA, SVR, LSTM are compared and evaluated for least root mean square error.
• Most of the research work for forecasting based models is developed on standard datasets available online, but less work on real datasets has been done which can provide more insight into the real time operating situation. Regarding this, our proposed model is evaluated on a real dataset from a local factory, so that it can be used in industry.

The remaining paper structure can be described as follows. Section 2 describes the mathematical derivation of the system model and unbalance estimation. Section 3 describes different forecasting models utilized to study the unbalance trend analysis. Section 4 describes the methodology for the prognostic approach developed in this study. Section 5 describes the results obtained and their subsequent discussions. Finally, Section 6 outlines the conclusions of the study. We also describe how this study can be utilized in the future for predictive maintenance purposes and the efficient reliable operation of a rotary system.

2. System Mathematical Model and Unbalance Estimation

Figure 1a represents the modeling of the rotor system with approximate dimensions and coordinates considered in this study. This model is developed from the real turbine rotor available in a local factory consisting of a flexible shaft, 3-stage impeller, two 5 tilting-pad hydrostatic fluid film bearings. The system model is mathematically derived using physics laws. Mostly, the rotor systems are derived using finite element model (FEM), but in this research assumed-mode method (AMM) is utilized to derive the equations of motion (EOMs) of the turbine rotor bearing system and the equations were formulated in a MATLAB program to solve for the parameters. The AMM based model is validated using a finite element based model and it is found to be acceptable for further analysis. Figure 1b shows the 3D model developed using Solidworks and Figure 1c represents its first natural frequency obtained using FEM. The natural frequencies obtained from both FEM and AMM models differ from each other by 2% (i.e., 84.4 and 83 Hz), under rigid support conditions. This rotor system generally runs at a speed of 11000 rpm but balanced at a lower speed of 600 rpm in each overhaul.

Based on the AMM, the displacements of rotor’s center line in x and y directions, measured with respect to the stationary coordinates, can be expressed as u(z, t) and v(z, t),

$$u(z,t)={\displaystyle \sum _{i=1}^m{p}_i(t){\mathsf{\Phi }}_i(z)}={\mathbf{p}}^{\mathrm{T}}\mathsf{\Phi }$$

(1)

$$v(z,t)={\displaystyle \sum _{i=1}^m{q}_i(t){\mathsf{\Phi }}_i(z)}=q^{\mathrm{T}}\mathsf{\Phi }$$

(2)

Wherein, assumed modes i.e., Φ_i(z), consists of the natural modes of a free-free beam plus two rigid body modes, translation, and rotation. p and q are generalized coordinates and m is the number of modes deemed required accuracy.

The EOMs, via the Lagrange’s equation, are derived to be,

$$\left[\begin{array}{cc}M+N& 0\\ 0& M+N\end{array}\right]\left\{\begin{array}{c}\ddot{p}\\ \ddot{q}\end{array}\right\}+\left[\begin{array}{cc}{C}_x& 2\mathsf{\Omega }N\\ -2\mathsf{\Omega }N& C_y\end{array}\right]\left\{\begin{array}{c}\dot{p}\\ \dot{q}\end{array}\right\}+\left[\begin{array}{cc}{K}_x-{\mathsf{\Omega }}^{2}N& 0\\ 0& K_y-{\mathsf{\Omega }}^{2}N\end{array}\right]\left\{\begin{array}{c}p\\ q\end{array}\right\}=\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}$$

(3)

$$M_{ij}={\displaystyle {\int }_0^l\rho A{\mathsf{\Phi }}_i}(z){\mathsf{\Phi }}_j(z)dz+{\displaystyle \sum _{r=1}^{n_d}{m}_r^d}{\mathsf{\Phi }}_i(z_r){\mathsf{\Phi }}_j(z_r)$$

(4)

$$N_{ij}={\displaystyle {\int }_0^l\rho I_A{{\mathsf{\Phi }}^{\prime }}_i}(z){{\mathsf{\Phi }}^{\prime }}_j(z)dz+{\displaystyle \sum _{r=1}^{n_d}{I}_r^d}{{\mathsf{\Phi }}^{\prime }}_i(z_r){{\mathsf{\Phi }}^{\prime }}_j(z_r)$$

(5)

$$C_{x_{ij}}={\displaystyle \sum _{n=1}^{n_b}{C}_{xn}{\mathsf{\Phi }}_i(z_n){\mathsf{\Phi }}_j(z_n)}$$

(6)

$$C_{y_{ij}}={\displaystyle \sum _{n=1}^{n_b}{C}_{yn}{\mathsf{\Phi }}_i(z_n){\mathsf{\Phi }}_j(z_n)}$$

(7)

$$K_{x_{ij}}={\displaystyle {\int }_0^{l}E{I}_A}{{\mathsf{\Phi }}^{″}}_i(z){{\mathsf{\Phi }}^{″}}_j(z)dz+{\displaystyle \sum _{n=1}^{n_b}{k}_{xn}{\mathsf{\Phi }}_i(z_n){\mathsf{\Phi }}_j(z_n)}$$

(8)

$$K_{y_{ij}}={\displaystyle {\int }_0^{l}E{I}_A}{{\mathsf{\Phi }}^{″}}_i(z){{\mathsf{\Phi }}^{″}}_j(z)dz+{\displaystyle \sum _{n=1}^{n_b}{k}_{yn}{\mathsf{\Phi }}_i(z_n){\mathsf{\Phi }}_j(z_n)}$$

(9)

Here M, N, C, and K, respectively, represent the mass, rotary inertia, damping, and stiffness matrix respectively. Also, $(\dot{})=d/dt, (\prime )=d/dz$; n_d, m^d, I^d, n_b, and I_A denote the number of disks, disk’s mass, disk’s diametric inertia, number of bearings, and shaft’s area moment of inertia, respectively. $\mathsf{\Omega }$ is the rotational speed in rad/s, $\rho$ is the mass density in kg/m³ and Q_x, Q_y are generalized force vectors.

The unbalanced (centrifugal) forces when decomposed into x and y components can be expressed as,

$$F_x=U_g{\mathsf{\Omega }}^2\cos (\mathsf{\Omega }t+{\alpha }_g)\delta (z-z_g)$$

(10)

$$F_y=U_g{\mathsf{\Omega }}^2\sin (\mathsf{\Omega }t+{\alpha }_g)\delta (z-z_g)$$

(11)

where $U_g=(m\cdot e)$ is the customarily specified unbalance with a unit of gmm, e is the eccentricity i.e., the distance of displaced mass centre from the geometric center due to unbalance forces. α_g is the unbalance phase angle relative to the key phasor. Here, the subscript g denotes the gravity center.

The generalized force vectors are calculated as,

$$\begin{array}{l}{Q}_x=U_g{\mathsf{\Omega }}^2\cos {\alpha }_g\mathsf{\Phi }(z_g)\cos (\mathsf{\Omega }t)-U_g{\mathsf{\Omega }}^2\sin {\alpha }_g\mathsf{\Phi }(z_g)\sin (\mathsf{\Omega }t) \\ \hspace{1em}=F_c(z_g)\cos (\mathsf{\Omega }t)-F_s(z_g)\sin (\mathsf{\Omega }t)\end{array}$$

(12)

$$\begin{array}{l}{Q}_y=U_g{\mathsf{\Omega }}^2\sin {\alpha }_g\mathsf{\Phi }(z_g)\cos (\mathsf{\Omega }t)+U_g{\mathsf{\Omega }}^2\cos {\alpha }_g\mathsf{\Phi }(z_g)\sin (\mathsf{\Omega }t) \\ \hspace{1em}=F_s(z_g)\cos (\mathsf{\Omega }t)+F_c(z_g)\sin (\mathsf{\Omega }t) \end{array}$$

(13)

Substituting the values of Q_x and Q_y from (12), (13) in the EOM (3), and solving for cosine and sine components of generalized components which can be evaluated from the following equation,

$$\left[\begin{array}{cccc}{K}_x-{\mathsf{\Omega }}^{2}M-2{\mathsf{\Omega }}^{2}N& \mathsf{\Omega }{C}_x& 0& 2{\mathsf{\Omega }}^{2}N\\ -\mathsf{\Omega }{C}_x& K_x-{\mathsf{\Omega }}^{2}M-2{\mathsf{\Omega }}^{2}N& -2{\mathsf{\Omega }}^{2}N& 0\\ 0& -2{\mathsf{\Omega }}^{2}N& K_y-{\mathsf{\Omega }}^{2}M-2{\mathsf{\Omega }}^{2}N& \mathsf{\Omega }{C}_y\\ 2{\mathsf{\Omega }}^{2}N& 0& -\mathsf{\Omega }{C}_y& K_y-{\mathsf{\Omega }}^{2}M-2{\mathsf{\Omega }}^{2}N\end{array}\right]\left\{\begin{array}{c}{p}_c\\ p_s\\ q_c\\ q_s\end{array}\right\}=\left\{\begin{array}{c}{F}_c\\ -F_s\\ F_s\\ F_c\end{array}\right\}$$

(14)

The above equation gives complete solution to an unbalance excitation. The rotor’s unbalance responses in x and y directions, after solving Equation (14), can be expressed as

$$\begin{array}{ll}u(z,t)& =p_c^{\mathrm{T}}\mathsf{\Phi }(z)\cos (\mathsf{\Omega }t)+p_s^{\mathrm{T}}\mathsf{\Phi }(z)\sin (\mathsf{\Omega }t)\\ & =u_c(z)\cos (\mathsf{\Omega }t)+u_s(z)\sin (\mathsf{\Omega }t)=A_x\cos (\mathsf{\Omega }t-{\varphi }_x)\end{array}$$

(15)

$$\begin{array}{ll}v(z,t)& =q_c^{\mathrm{T}}\mathsf{\Phi }(z)\cos (\mathsf{\Omega }t)+q_s^{\mathrm{T}}\mathsf{\Phi }(z)\sin (\mathsf{\Omega }t)\\ & =v_c(z)\cos (\mathsf{\Omega }t)+v_s(z)\sin (\mathsf{\Omega }t)=A_y\cos (\mathsf{\Omega }t-{\varphi }_y)\end{array}$$

(16)

The calculated unbalance vibration from above equation is compared with the measured response from the sensors to generate a residual and, model parameters is estimated by minimizing the residual value.

$$\{\begin{array}{c}u(z_i,t)=u_c(z_i)\cos (\mathsf{\Omega }t)+u_s(z_i)\sin (\mathsf{\Omega }t)\\ v(z_i,t)=v_c(z_i)\cos (\mathsf{\Omega }t)+v_s(z_i)\sin (\mathsf{\Omega }t)\end{array} , z_i=z_1 \mathrm{or} z_2$$

(17)

The subscript 1 and 2 at z denotes the two different bearing locations. The unbalance mostly affects vibration components at 1× (rotational speed’s frequency), hence real-time 1× components dispalcement values are extracted from the sensor data at the two bearing locations to generate a residual with the above calculated value. Thus, the feature vector for this is given by,

$$f={\{u_{c1},u_{s1},v_{c1},v_{s1},u_{c2},u_{s2},v_{c2},v_{s2}\}}^T$$

(18)

and residual error is given by,

$$E={\displaystyle \sum \parallel (f_{\mathrm{T}}-}{f}_{\mathrm{m}})\parallel$$

(19)

where, $f_{\mathrm{T}}$, and $f_{\mathrm{m}}$, respectively denote the calculated and measured feature vector and the unbalance vector d to be identified at the two bearing planes is.

$$d={\{U_1,{\alpha }_1,U_2,{\alpha }_2\}}^T$$

(20)

The optimization problem can hence be formulated as

$$\{\begin{array}{l}\min \mathrm{E}={\displaystyle \sum _{i=1}^8|}|f_{\mathrm{T}}(d)-f_{\mathrm{m}}(d\left)\right||\hfill \\ \mathrm{subject} \mathrm{to} \hfill \\ U_1\le 1.2U_1(\mathrm{latest} \mathrm{one})\hfill \\ \begin{array}{l}{U}_2\le 1.2U_2(\mathrm{latest} \mathrm{one})\hfill \\ \begin{array}{l}{\alpha }_1\in [{\alpha }_1^L,{\alpha }_1^R]\hfill \\ {\alpha }_2\in [{\alpha }_2^L,{\alpha }_2^R]\hfill \end{array}\hfill \end{array}\hfill \end{array}$$

(21)

$U_1$, $U_2$ are then calculated as the sum total unbalance acting at bearing locations 1 and 2 respectively and ${\alpha }_1$, ${\alpha }_2$ is their corresponding phase angles. An optimization problem was created to determine the bearing parameters C and K for the system as shown in Equation (21). Note that the constraint of unbalance quantities $U_1$ and $U_2$ less than 20% increment of the latest ones is deliberately set up for computational efficiency, because unbalance deterioration hardly exceeds 20% over one-day operation unless a sudden brake down occurs. The boundaries for phase angles are calculated according to the measured vibrations and details can be referred to [34].

Equation (21) is then solved for unbalances at two bearing planes d by the particle swam optimization (PSO) [35]. The PSO method was found to give the most robust and fast results when compared to five other optimization methods. As long as d is identified, the unbalances can be configured into the so-called static ($U_g$, ${\alpha }_g$) and dynamic ($U_d$, ${\alpha }_d$) unbalance, which respectively denote the offset of gravity center and the angle of principal rotary inertia axis away off the rotation axis. More detailed analysis on the methodology can be referred to [34]. In this paper, a physics-based model is built based on mathematical derivations and parameter identification by PSO, the model was further validated using real data from the industry.

A rotor’s balance would gradually deteriorate with time with load acting on it and its continuous operation [36]. Both static and dynamic unbalance forces will keep on increasing. A continuous monitoring approach can then be applied to collect the data and calculate the unbalances to determine the trend for prognosis.

3. Forecasting Methods

In this paper, we propose various classical and non-classical/machine learning-based forecasting methods that can be employed further in combination with our mathematical model to establish a real-time unbalance prognosis system, which would be helpful in determining the maintenance strategy of the rotor system. Assesment of the onset of fault, is one of the imporant aspects of PHM and efficient condition monitoring of the overall system. After the identification of fault, understanding its trend from past data and prediction of failure in the future, also has significant consequence. Given below are the brief descriptions of all the methods that have been used.

3.1. Classical/Statistical Methods

3.1.1. Simple Exponential Smoothing (SES)

This is one of the oldest and simplest statistical forecasting technique which has been in use for a long time, the reason being the calculations are easier and gives a relatively good result in no time [37]. In this type of forecasting method, weights are assigned on past value which decreases exponentially over time. The equation describing this method is given below,

$$f_t=\alpha d_{t-1}+(1-\alpha )f_{t-1}$$

(22)

where the value of α varies from 0 to 1. $f_t$ is the forecast value for time t and $f_{t-1}$ is the forecast value for time t − 1.

3.1.2. Auto Regressive Integrated Moving Average (ARIMA)

In real life, datasets are mostly non-stationary in nature and in order to forecast using ARIMA model, it is first required to make them stationary. ARIMA model [38] is an extension of ARMA model wherin there is an additional integration step which is achieved by differencing the time-series. Differencing refers to subtracting the previous value from the current value. ARIMA model is based on its own lags and lagged forecast errors, so that its equation can be used to forecast future values, given if the time series is stationary. Since ARIMA modeling deals with a stationary time series, if the time series is not stationary, they are made stationary first with help of differencing. ARIMA Model Equation can be described as,

$$Y_t=\alpha +{\beta }_1{Y}_{t-1}+{\beta }_2{Y}_{t-2}+\dots +{\beta }_p{Y}_{t-p}{\epsilon }_t+{\phi }_1{\epsilon }_{t-1}+{\phi }_2{\epsilon }_{t-2}+\dots +{\phi }_q{\epsilon }_{t-q}$$

(23)

where $Y_{t-p}$ represents the lags of series and ${\beta }_p$ is the coefficient of lag, $\alpha$ is the intercept term, $\epsilon$ is the error of autoregressive models of the respective lags and ${\phi }_q$ is the coefficient of errors.

$$\mathrm{Predict} Y_t=\mathrm{Constant}+\mathrm{Linear} \mathrm{combination} \mathrm{Lags} \mathrm{of} Y (\mathrm{up} \mathrm{to} p \mathrm{lags})+\mathrm{Linear} \mathrm{Combination} \mathrm{of} \mathrm{Lagged} \mathrm{forecast} \mathrm{errors} (\mathrm{up} \mathrm{to} q \mathrm{lags}).$$

3.2. Machine Learning Methods

3.2.1. Support Vector Regression (SVR)

To solve non-linear regression problem, inputs are non-linearly mapped into a high dimensional feature space where the outputs are related linearly, given by,

$$y=f(x)=w\cdot \phi (x)+b$$

(24)

Here y is the model output corresponding to input x which is mapped into feature space by a non-linear function $w\cdot \phi (x)$, where w is the weight vector and $\phi (x)$ is the mapping function in the feature space. The regression problem for this, can be expressed as the following optimization problem given by,

$$\begin{array}{c}\underset{w,\epsilon ,{\epsilon }_i}{\min }\frac{1}{2}{w}^2+C{\displaystyle {\sum }_{i=1}^l({\xi }_i}+{\xi }_i^{*})\\ \mathrm{Subject} \mathrm{to}:\\ y_i-(w\cdot x_i+b)\le \epsilon +{\xi }_i\\ (w.x_i+b)-y_i\le \epsilon +{\xi }_i^{*}\\ {\xi }_i,{\xi }_i^{*}\ge 0,i=1,2,\dots ,l\end{array}$$

(25)

For a more detailed explanation of mathematical formulations in SVR can be referred to from the article [39].

3.2.2. Long Short-Term Memory (LSTM)

LSTM is an advancement to recurrent neural networks (RNN) where it addresses the problem of vanishing gradient and is able to remember long term dependencies of data, initially proposed by Sepp Hochreiter and Jurgen Schmidhuber in 1997 [40]. LSTM generally consists of memory cells which are handled by three special gates, input gate, output gate and forget gate. LSTM is found to be useful when building sequence to sequence architectures [41].

For a time series network at time t, if $x_t$ is the input data for LSTM cell, $h_{t-1}$ is the output of previous cell, $c_t$ is the memory value of current memory cell, $c_{t-1}$ is the memory value of previous memory cell and $o_t$ is the output of LSTM cell, then the calculations are given by following equations,

$$c_t=\tanh (w_c[h_{t-1},x_t]+b_c)$$

(26)

$$i_t=\sigma (w_i[h_{t-1},x_t]+b_i)$$

(27)

$$f_t=\sigma (w_f[h_{t-1},x_t]+b_f)$$

(28)

$$c_t=f_t\ast c_{t-1}+i_t\ast c_t$$

(29)

$$o_t=\sigma (w_o[h_{t-1},x_t]+b_o)$$

(30)

$$h_t=o_t\ast \tanh (c_t)$$

(31)

where $w_i$, $w_c$, $w_f$, $w_o$ are the weights, $i_t$,$f_t$ and $o_t$ represents input gate, forget gate and output gate respectively. $\tanh$ and $\sigma$(sigmoid) are the two activation functions used for calculations in a simple LSTM cell. Figure 2 shows a unit cell representation of LSTM network with various gates.

4. Proposed Prognosis Methodology Framework

As technology advances with Industry 4.0, there has been a need of establishing systems with early fault detection and failure prognosis. The main reason for prognostics system generally is to determine or predict with accuracy the exact time when the component will fail, i.e., determining its RUL. However, in the case of unbalance prognosis of a rotor bearing system, for example the model given in Figure 1, it is difficult to determine the RUL due to unbalance itself, as the turbine shaft will not normally fail only due to the unbalance forces; even though it may lead to secondary faults as well as increase the vibration of the system and reducing overall system efficiency. Thus, normally no sufficient historical data exists for failure due to unbalance, hence, it becomes absolutely necessary to balance the rotary system from time-to-time for safe, reliable and efficient operation of the system. If we can evaluate the trend of unbalance and determine the time or threshold limit when the balancing has to be done, then the overall system can perform efficiently and more reliably as well as it can prove to be more significant for the industry.

As we know, balancing of a rotating system improves its service life, makes the system more efficient with reliable operation, aids in smooth running and reduces overall vibration and unwanted noise in the system. Hence, it becomes essential to continuously monitor the unbalance in the system and check for its trend analysis. So, our objective is to evaluate a good prognostic approach with the aid of monitoring, which can give month ahead predictions of unbalance with acceptable accuracy, so that further balancing-based maintenance steps can be taken according to the forecast analysis. The proposed methodology employed is described in Figure 3.

The main steps to employed in our prognostic approach are described as follows,

Step 1: Data Acquisition: Sensors are employed for the collection of data from the bearing locations of the system during its operation.

Step 2: Model-Based Unbalance Calculation: The static and dynamic unbalances are calculated from the derived mathematical model methodology as stated above in Section 2 and can refer further in details in paper [34]. The unbalance dataset (U_g, α_g, U_d, and α_d) evaluated from the model-based calculations consists of 489 days of data which varies with respect to time. This data is considered raw data for further processing of trend analysis.

Step 3: Data Preprocessing: The raw data generated is processed by moving an average filter over it to remove the noise and smoothen the data. The smoothening of data is done using the function “smooth” (MATLAB Curve Fitting Toolbox, function “smooth” with “rlowess” which is a robust linear regression function over a specified window). A window of 15 is used in this case for smoothening. Figure 4 describes the plots of smoothened dataset evaluated from raw static and dynamic unbalance (U_g, α_g, U_d, and α_d) datasets, calculated from Equation (21). The curves in Figure 4 look jagged due to some factors, such as loading fluctuations in operation, noises in measurement etc., but mostly the fluctuations caused by shaft torsional vibration, which has not been considered in bending vibration analysis. Torsional vibration will affect the key-phasor timing and subsequently result in variations of vibration measurements. These factors, however, do not influence the unbalance in a sense of long-term running. In order to prevent convergence issues and eliminate the effect of different scales, input and output variables were then normalized. At each run, the mean and standard deviation of the training set were calculated, and normalized data was fed to the models for evaluation. Equation (32) shows the equation for normalization where X is the input and X_t is the normalized output. μ and σ are mean and standard deaviation respectively.

$$X_t=\frac{(X-\mu )}{\sigma }$$

(32)

Step 4: Models’ Hyperparameters Selection: There are several parameters in each model that were varied to arrive at the best model parameter for each case. The proposed framework is modeled in MATLAB software. The description of different hyperparameters varied in all the models is shown in

Table 1. Hyperparameters taken into consideration for different models.

Model	Hyperparameters	Values
ES	Smoothing Parameter	[0.1, 0.2, …, 0.9]
ARIMA	(p, d, q)	D = 0, 1, and p, q = 0, 1, 2, 3
SVR	ε	0.1111
	Kernel	gaussian
	n-fold cross validation	10
LSTM	Hidden units	10
	Epochs	200
	Dropout	0.2
	Optimizer	Adam
	Initial Learn Rate	0.005
	Mini-Batch Size	128
	Learn Rate Drop Factor	0.1
	L2 Regularization factor	0.0001
	Gradient decay factor	0.1
	Gradient threshold	1

.

The smoothing parameter in ES method is varied from 0.1 to 0.9, to determine which one performs better for the given dataset.

The variables for ARIMA, i.e., p, d, q were also varied for evaluation. For modeling of ARIMA, the MATLAB function arima is used with gaussian conditional probability distribution.

For SVR modeling, MATLAB function fitrsvm is used with the gaussian kernel function to optimize the SVR. A ten-fold cross-validation is performed to remove overfitting in the model.

Hyperparameter tuning in LSTM: Figure 5 shows the different layer structures in developing LSTM forecasting model. The LSTM model is sequence to sequence structure, where preprocessed data after normalizing is used as a sequential input, LSTM layer with hidden units of 10 is used to map the output and a dropout probability of 0.01 is used to avoid overfitting in the training, finally a fully connected layer, regression layer is added to extract the output from the trained network. Adam optimizer is used since it performed better than sgdm and rmsprop optimizer. The network was run for a maximum of 200 epochs. To avoid the problem of gradient exploding, the gradient threshold is set to 1. The loss function in the regression layer is mean squared error, MSE, and training stops when MSE is minimum. An L2 regularization factor of 0.0001 is added to the model. The initial learning rate was set at 0.005.

Training architecture for machine learning models: the dataset is divided into two XTrain and YTrain; in which XTrain ranges from start to end –− lag of time-series and YTrain ranges from lag + 1 to end of time-series.

Step 5: Forecasting using the recursive method: Several forecasting methods were employed to perform the forecasting analysis using a recursive strategy. Recursive strategy refers to the estimation of the next day’s forecast given by the past values and then feeding back the forecasted value for another day’s forecast. Both rolling and expanding window were evaluated and we found expanding window to have better forecast results, so we have used expanding window to develop our model. Thus, a p-step ahead forecast using recursive strategy can be represented by Figure 6. In this figure, X_i represents data values at time i, shown in yellow color and blue color block is the model which is used to forecast day ahead values.

Step 6: Numerical Evaluation: To validate our forecasting model, we also need to check for its accuracy which will further determine whether our model can be effectively used for a real-world scenario with precision or not. For this reason, we have used root mean square error for our evaluation.

Root mean squared error (RMSE): RMSE is one the most popular performance measurement techniques commonly employed to check the error between predicted and original values. The RMSE error can be defined by the following mathematical expression,

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{(y_i^{obs}-y_i^{pred})}^2}}$$

(33)

where, $y_i^{obs}$ is the current observation and $y_i^{pred}$ is the forecasted value from model. This RMSE is used at last after the entire prognostic model is built for comparision about which model will perform better under different training criterion, to better evaluate the real industrial situation.

5. Results and Discussions

We have derived four different time-series datasets for our prognostic framework evaluation namely, U_g, α_g, U_d, and α_d utilizing the sensor’s data and mathematical model which span over 489 days. Here, the g subscript refers to static unbalance acting at its center of gravity, while the d subscript refers to dynamic unbalance acting at its center of gravity. The datasets are univariate time series i.e., they vary with time and are extracted from a local factory recorded in the year 2020 to 2021 out of the company’s database.

For establishing a forecasting problem and training the dataset in real-time, we first selected to train an initial window of certain days of data and continuously forecasting 30 days at a time and expanding the training set simultaneously i.e., expanding window method, or moving the training set with a fixed window i.e., sliding window method and forecasting 30 days again while tuning the different model parameters to achieve as much accuracy as possible. Finally, evaluating with different forecast modeling methodologies to check which works best for the rotating system dataset extracted from the industry. Figure 7 illustrates the different approaches where training is the input time series, and the test is the forecast computed from the training set.

We chose 30 days to make monthly ahead unbalance predictions as it might be more useful for consideration in the future for predictive maintenance purposes.

Datasets used: unbalance values from the calculations are,

U_g (static unbalance)—489 days of data,

α_g (phase angle of U_g)—489 days of data,

U_d (dynamic unbalance)—489 days of data,

α_d (phase angle of U_d)—489 days of data.

Different forecasting models are built to address the forecasting problem. All the methods were modeled with the help of MATLAB R2022a programmed on an intel core i7-8700 CPU (3.1GHz), 16 GB RAM, NVIDIA GeForce GTX1650 GPU Computer.

For the selection of exponential smoothing parameter in the simple exponential model, it is varied from 0 to 1 and the one with least RMSE is selected for modeling. Figure 8 shows the variation of RMSE on U_g when smoothing parameter, α is varied. From this figure, we can observe that at smoothing parameter 0.1 gives the least RMSE value for exponential smoothing model.

For the selection of the ARIMA model, autocorrelation function (ACF) and partial autocorrelation function (PACF) figures of the original time series are evaluated for the stationarity of our time-series data. Autocorrelation refers to the degree of correlation the time series has with the lagged version of itself over subsequent time intervals. The lag values help to understand which lag has a more significant correlation in a time series, which in turn helps in determining the more fitting model parameters. PACF on the other hand, gives partial correlation with its own lagged values and regressively controls the value of shorter lags. Figure 9 and Figure 10 show the original ACF and PACF curves of time-series data. From Figure 9, we have found that there is a slight trend in the data, also to confirm that the series is not stationary, augmented dickey fuller (ADF) test is carried out. ADF test is a kind of unit root test, in which if there is no root then the time-series is stationary otherwise it is non-stationary. The ADF test result justifies that the time series is not stationary, since the p-value of ADF test for d = 0, for unbalance datasets (U_g, α_g, U_d, α_d) is less than 0.05. To make the time-series stationary, data is transformed using first-order regular differencing of the series. After the transformation, it is checked for stationarity again using the ADF test and it showed the p-value more than 0.05 which means time series is now stationary.

Table 2. ADF test results.

Time–Series	p Value	Result
Ug (d = 0)	0.6235	Not Stationary
Ug (d = 1)	0.001	Stationary
Ud (d = 0)	0.4042	Not Stationary
Ud (d = 1)	0.001	Stationary
αg (d = 0)	0.9648	Not Stationary
αg (d = 1)	0.001	Stationary
αd (d = 0)	0.8474	Not Stationary
αd (d = 1)	0.001	Stationary

shows the ADF test values. Figure 11 and Figure 12, show ACF and PACF curves of differenced time series. From the Figure 11, upto lag 5 is more significant, whereas in Figure 12 upto lag 4 is significant. To calculate for p, q parameters and reduce the comprexity of modeling, we chose lags 1, 2, 3, to check for the best p, q values. Since, the first-order difference of the time series is established, to choose the best p and q parameter of ARIMA (p, d, q) modeling, akaike information criterion (AIC), bayesian information criterion (BIC) values are calculated for different combinations as represented in

Table 3. AIC, BIC values of different ARIMA models.

Arima Models	AIC	BIC
ARIMA (2,1,1)	1668.5	1671.2
ARIMA (2,1,2)	1641.4	1644.1
ARIMA (2,1,3)	1638.2	1640.8
ARIMA (3,1,2)	1638.6	1641
ARIMA (3,1,3)	1636.4	1638.8

. AIC is used to measure the model for analytical fitting whereas BIC is used to prevent the model for becoming too complex. The expressions of AIC, BIC are described below,

$$\mathrm{AIC}=-2\log (\mathrm{maximumlikelihood})+2k$$

(34)

$$\mathrm{BIC}=-2\log (\mathrm{maximumlikelihood})+\frac{k\log n}{n}$$

(35)

where k is independently adjusted number of parmeters and n is number of data points.

The lowest criteria value of AIC, BIC generally corresponds to the best ARIMA fit model. Since, AIC, BIC value of ARIMA (3,1,3) is found to be the lowest, so it was selected as the forecasting model for the given dataset.

Similarly, for the selection of LSTM model several variations of hidden sizes, number of layers, dropout layer probabilities, epoch numbers are tried using trial and error to establish a perfectly suiting forecast method. RMSE is selected as the performance criteria to choose the best model out of others. Figure 13 shows the training performance curve of LSTM, it shows that the LSTM converges after a suitable number of iterations with less loss. Hence, the LSTM model is fit to make future predictions.

After the selection of various models, the initial window length of training is varied both for sliding window and expanding window methodologies, to understand the effect of windowing for training in the forecasting trend analysis as well as perform the real-time trend analysis and it was evaluated using the least RMSE criterion. It is observed that expanding window gives less RMSE values in all the models. This could be due to the fact that our data is initially less, and it needs more data for better future predictions. Figure 14 and Figure 15 show the effect of windowing in the overall RMSE values both for sliding and expanding windows of the respective models. The windowing does not affect the simple exponential smoothing method because the forecast value in this directly depends on single past value. In both the figures, red line has been used for ARIMA with asterisk symbol at different windows, blue line is for SVR with circle symbol, yellow line LSTM with plus symbol, purple line is for ES with square symbol respectively. From the figure, it can be seen that most of the time, ARIMA performs better in forecasting performance according to root mean square error (RMSE) evaluation.

Figure 16a–d shows the overall variations of different models forecasting unbalance trend predictions between the real and ES, ARIMA, SVR, and LSTM models. The red line corresponds to the real values whereas ES model is shown by grey dotted, ARIMA model with green, SVR model with orange dashed and LSTM model with a blue line. Note that, the curves here are different model’s 30 days ahead moving forecast values using expanding window method and considering the initial window length of 200. The training data starts at 200 and keeps on increasing as it moves forward by taking the 30 days original data as the data becomes available. The overall comparison of model performances in terms of RMSE is shown in

Table 4. Overall RMSE performance comparison on different models (red color represents the best value).

	Window	ES	ARIMA	SVR	LSTM
Ug	Sliding	44.8961	38.5293	39.5695	43.4528
	Expanding	44.8961	37.8010	37.7116	38.3049
αg	Sliding	25.1188	24.3614	29.3558	26.3688
	Expanding	25.1188	22.0240	28.4961	26.5326
Ud	Sliding	20.6870	16.4772	18.9522	20.3374
	Expanding	20.6870	14.5949	17.7584	20.5
αd	Sliding	2.7578	2.386	3.2505	3.05
	Expanding	2.7578	2.3345	2.8	3.137

. From the figures and table, it is clear that ARIMA model gives best prediction. Even though ARIMA is a simple model, it generally gives a good performance in many scenarios where complex machine learning based systems fail to give good prediction due to highly non-linear data structure. However, it is always reasonable to construct machine learning models which have the capability to study non-linear relationships of a model parameters and perform reasonably well in the long run and evaluate long time dependencies on itself. In Table 4, the overall RMSE evaluation with sliding and expanding windows of different models is represented, here the red color represents the lowest RMSE values. However, in our observation, SVR also gives satisfactory results and is the second best model. Although LSTM is recently a popular machine learning forecasting tool, yet the LSTM results highly depend on good hyperparameter tuning and availability of larger dataset to learn the long term dependencies. Since there is no set rule to determine the best parameter setting for any dataset, a lot of times, trial and error or optimization algorithms are utilized to determine those parameters. Even though, we have used several different parameter combinations to find the most suitable LSTM forecasting model for determining our prognostic model system, there is still a probability of finding a better machine learning model for predictions since it has the capability to learn non-linear relationships of a trend very well. Figure 17 shows the monthly-ahead forecasts using ARIMA methodology and expanding window. An initial window of 200 days is selected for training and the window keeps expanding by 30 days data with every month as time progresses. This prognostic methodology can be utilized for further setting up of a threshold limit with the help of continuous monitoring using this strategy to evaluate the necessary maintenance based actions.

6. Conclusions

The accuracy of future predictions can significantly impact the overall decision-making of maintenance-based actions. For this, many state-of-the-art forecasting models based on statistical methods, or machine learning can play a pivotal role in analyzing real-time data and estimating future failure from the obtained data. In this research, we have successfully established a methodology to solve a forecasting problem that can benefit the industrial sector by monitoring the health conditions of the rotor-bearing system.

To achieve this, we have developed a prognostic approach for identifying the development of unbalance forces in the system with the aid of forecasting models such as ARIMA, ES, SVR, and LSTM, etc., and compared them to see what best fits our system with time-varying data. After that, we also developed our multi-step ahead forecasting algorithm that can update itself in real-time and give future predictions of unbalance with an acceptable range of prediction horizons. From all the observations we found that the ARIMA model with gaussian error performs very well to solve our problem of a prognostic framework, however SVR—a machine learning model—also gave a satisfactory performance. Furthermore, we tested the promising LSTM methodology which highly depends on the good tuning of hyperparameters and larger datasets, thus increasing model complexity. Even though in our present case, it didn’t provide satisfactory results; tuning the hyperparameters properly or combining the model with other methodologies may provide better forecasting results. The main purpose of this research is to propose a methodology for a prognostic framework of unbalance based on model-based and statistical/machine learning methodologies and compare the different methodologies to derive a real-time framework which can help in evaluating future fault trend in the system with the least RMSE values. Even though, our main conclusion is that ARIMA is the best methodology for the present case, which is a rotor-bearing system, it might not be so for every case as the dataset can be different on a different system or different kind of fault trend evaluation. We would like to contemplate more in detail about this in the future.

Also, better tuning of parameters and combinations of models can be explored to extract a better prediction model. Further modifications can be made to the modeling to make the predictions better. One of the main interests in the development of this prognostic system is to aid the decision-making process for effective and reliable maintenance in a rotary system by conducting necessary balancing of the system after it achieves a threshold for more efficient operation of the overall system with less downtime. The model developed can be used in the future to determine the timeline when the balancing of the rotor must be done for a smooth operation of the overall system as well as what implications do unbalance forces hold in the overall fault diagnosis. With the availability of failure data, more accurate evaluations for scheduling time for rotor balancing can be established and that is also a motivation for future research direction.