Differential Equation-Based Prediction Model for Early Change Detection in Transient Running Status

Early detection of changes in transient running status from sensor signals attracts increasing attention in modern industries. To achieve this end, this paper presents a new differential equation-based prediction model that can realize one-step-ahead prediction of machine status. Together with this model, an analysis of continuous monitoring of condition signal by means of a null hypothesis testing is presented to inspect/diagnose whether an abnormal status change occurs or not during successive machine operations. The detection operation is executed periodically and continuously, such that the machine running status can be monitored with an online and real-time manner. The effectiveness of the proposed method is demonstrated using three representative real-engineering applications: external loading status monitoring, bearing health status monitoring and speed condition monitoring. The method is also compared with those benchmark methods reported in the literature. From the results, the proposed method demonstrates significant improvements over others, which suggests its superiority and great potentials in real applications.


Introduction
Real-world systems are seldom in steady status, and they almost always operate in transient conditions that are varying over time. Detection of change(s) in running status at an early stage enables to find and locate abnormal and unexpected/undesired system behavior(s) in its successive operations. With this, corrective schedule could be made to prevent potential operation failures and ensure equipment and product reliability, safety, quality, productivity, etc. [1][2][3][4]. The problem of running status monitoring may arise in process monitoring and control applications where the system needs to make a response and/or take appropriate actions as soon as possible after a change of machine status occurs [5][6][7][8].

Differential Equation-Based Prediction Model
This section first provides an overview of the proposed prediction model, and then discusses techniques carried out in the main steps.

Overview of the Proposed Model
Given a data stream composing of continuously monitoring condition signal x t , let us assume that it has been denoted by a periodic form x t → X nT+v , where T is the period, v ∈ [1, T] is the phase or timing, and n is the cycle index. In this paper, we assume the condition signal is periodically stationary as made in some previous works, e.g., [25,39], thus we can employ a prior periodicity estimation proposed in [18] to confirm the value of T. Then, main steps of the proposed prediction model includes the following: • Model formulation, i.e., the proposed model is formulated with the considered CM signals. In this paper, a family of new time series are formed by arranging the original data at the same phase.
As such, the model is formulated so as to predict next value of each phase.

•
Parameters estimation, which estimates the parameters of the model. The numerical solution method of differential equations is used to estimate the model's parameters. The parameters of the model are constantly updated during each data prediction process. • Data prediction, i.e., the prediction of next data with the estimated model. The prediction value at each phase can be obtained with the model whose parameters have been estimated successfully.
The above steps are executed continuously for each observed cycles of CM data, allowing for residual error analysis that can quantify the error between the predicted values and the observed ones. The resulting error reflects the amount of temporal fluctuations in the CM data; more specifically, the higher residual error implies a higher probability that a change occurs and vice versa.
In the following, techniques used in the main steps of the model are discussed.

Model Formulation
For an arrived CM data stream up to the kth cycle, for each phase v, we take i continuous cycles before k to form a new series, i.e., P v = {x (k−i)T+v , ..., x (k−1)T+v }, such that the DE model can be formulated for all phases.
Conventionally, the DE model relies on the following system state function f (·) to establish the prediction, where x t is the value of differential equation in a series of time stamp t, β is parameter of the differential equation.
In this paper, since we use the past observed data at the same phase for prediction, the above-given new time series P v allows linear characteristics to emerge, based on which, we can obtain the prediction value, i.e.,x (k)T+v with the P v correspondingly. In the following, for each considered phase v, we re-denote the P v as the series {x t j }, j = 1, ..., i for simplicity.
According to the linear assumption made already, the time-series {x t j } will satisfy, where a( = 0) and b( = 0) are the parameters of differential equations that need to be estimated before prediction.

Data Prediction
Assuming that we have estimated the values of a and b, according to [40], the corresponding time response function of Equation (1) can be written aŝ where t i = t 0 + i × h, t 0 = 0, t i+1 is the time stamp of the value x t i+1 , which can also be regarded as the prediction valuex (k)T+v , and h is the distance between the two time stamps. We set h as 1 in this article. The x t i is the initial value of the response function, represented by the average value of the sequence value in order to decrease effect of errors in data collection.

Parameters Estimation
For the purpose of prediction ofx (k)T+v , the first step is to estimate the optimal parameters of a and b in the DE model given in Equation (2). Because the signal is usually sampled through observation, we re-define the above variables at sampling time t j (j = 1, ..., i) with subscription j (e.g., x t j = x j ) for convenience in the following of this paper. The differential Equation (1) is calculated by simpsons method [40] by in which f j = f (x j , u, t j ), j = 1, ..., i − 2 and R j+2 is local truncation error [40]. According to Equation (2), Equation (4) can be transformed into which is thus re-written as Then, Equation (6) can be re-expressed as where x 2 . When ∑ i−2 j=1 E 2 j+2 has a minimum value for the observed value x j , the estimated parameters a and b can be finally obtained by in which The predicted valuex (k)T+v can be obtained by substituting the estimated parameters a and b into the time response Equation (3), where i = 4 is used in the following experiment.

Residual Error Analysis
The residual error q (k)T+v can be defined by the absolute value of the difference between prediction datax (k)T+v and the actual monitoring data x (k)T+v as The cumulative value of {q (k)T+v } for a whole cycle is gathered, allowing for analyzing the machine condition periodically, which is calculated by where Q k is the resulting value for the kth cycle. Subsequently, we further standardize the obtained sequence {Q k } by s k = (Q k − Q)/σ, where Q and σ are the sample mean and standard deviation of {Q k }.
The anomaly score {s k } describes the dynamic characteristics of machine condition monitoring over time. That is, when the condition of machine is stable, this implies that the condition change does not occur, and the values of {s k } will be relatively small; otherwise, the values will be large.

Simulation Validation
For the purpose of evaluation of the effectiveness of the proposed model, we applied it on synthetically generated testing data, where different levels of Gaussian White Noise are added for simulating the real engineering scenarios. Considering that the status change(s) in real machine operations can result in the change(s) of relevant CM variables in terms of amplitude, frequency, or both of them [7,18,[41][42][43], the testing signals are formulated, respectively, as where c is the change time, m is the length of data, and ω is set as π. The results are given in Figure 1, where the Figure 1a-c shows the simulated signals generated by Equations (11)-(13), respectively. For each subfigure, the SNR are set as 10, 30 and 50 dB, respectively. It can be clearly seen that the anomaly scores are stable before the status change, and then the occurrence of status change causes an abrupt increase of the anomaly score, which implies that it can be detected successfully. Here, we also found that, for amplitude change and amplitude and frequency change cases, the residuals of 30 dB and 50 dB are different from the one of the case 10 dB. The main possible reason is that, with a lower SNR, i.e., 10 dB, the noise will be more taken into account for prediction of the future value with an amplitude change; however, with a higher SNR, the noise will have less influence on the perdition process.

Hypothesis Testing for Decision-Making
On the basis of the resulting anomaly score s k , we test a null hypothesis using the 3σ criterion in order to detect whether a change occurs at the current kth cycle by where H 0 means that no change occurs on the kth cycle as long as | s k − s k−1 |< 3σ , and H 1 indicates that a change occurs when | s k − s k−1 |≥ 3σ (Note that the beginning nine cycles of data are used for initialization and the change detection begins with the 6th cycle after initialization). Here, s k−1 and σ are the mean and standard deviation of an assumed Gaussian distribution, respectively, and they are calculated by Here, it is worth mentioning that there exist other alternatives such as Gaussian Mixed Model (GMM) [44] and other non-Gaussian assumptions [11,45] for change detection. However, since the focus of this paper is on the prediction model, we do not investigate these alternatives in this study.

Proposed Machine Running Status Monitoring Framework
Together with the DE model described in Section 2, an analysis of continuous monitoring of CM signals by means of a null hypothesis testing (Section 3) is proposed to inspect/diagnose whether an abnormal running status change occurs or not in successive machine operations. Total four steps are included in the framework as given below: (1) Collect CM data from the considered machine in a continuous manner; (2) Compute the prediction value using the proposed DE model; (3) Calculate anomaly scores based on residual error analysis at the current inspection time; (4) Make the change decision by testing a null hypothesis. Report an alarm to the user; otherwise, go to Step 2 to continue.
The flowchart of the proposed framework is provided in Figure 2.

Continuous signal collection
Signal source

Data acquisition card
On-line signal process

Experimental Validation
To evaluate the effectiveness of the proposed method, we applied it to three representative industrial applications/tasks which are listed as below: • External loading status monitoring: External loading status is essential for condition monitoring during unsteady machine operations because a piece of equipment under operation may be exposed to a series of varying loads according to the user's needs [46]. Moreover, the load is a critical operating condition factor which has significant impact on machine health [47]. Detection of changes in load condition makes it possible for the machine to adjust itself once a load change occurs for safety protection [48].
• Bearing health status monitoring: As we all know, functional degeneration of machine components during the lifetime is common and unavoidable. The component degeneration will cause some undesired/unexpected consequences [48][49][50]. Based on CBM, maintenance can be scheduled in an optimal way with respect to cost, reliability, availability, or other logistic metrics of interest. Thus, automatic detection of changes in bearing health status can serve as a starting point for fault diagnosis or prediction of functional failure at an early stage. • Rotational speed monitoring: The rotational speed in machine operations may fluctuate due to condition variations or unsteady environments [51]. Speed condition monitoring helps to find the unexpected running behaviors for operation maintenance [52,53], thus highly desired in online process monitoring of industrial manufacturing, numerical controlled machining, ect.
The proposed method is compared with the ARIMA model proposed in [18]. In addition, we also compare it with two methods: RMS and kurtosis, which are widely used in the literature [36][37][38]. In the following experiments, we utilize the first detected alarm to determine whether the detection is successful or false. We then use the precision indicator which is defined as ratio of the number of correctly detected changes over the total number of changes to quantify the detection performance. However, for the case only including several testing data, we directly provide the detection result.

Case Study I: External Loading Status Monitoring
The testing data used in this section are the motor bearing data provided by CWRU [33]. In the experiment, the vibration data were collected from the drive end and the fan end of the motor driving mechanical system with the sampling frequency of 12 kHz [54].
The testing data for loading change detection was formed by concatenating each two load condition signals with fixed other relevant variables. We use the L → L + ∆L to represent a loading status change, where L hp is the initializing load and ∆L hp is the changing load. More specifically, ∆L > 0 represents an increasing load change and ∆L < 0 represents a decreasing load change. The combination of L = {0, 1, 2, 3} hp and ∆L = {−3, −2, −1, 1, 2, 3} hp (indicating different loading conditions in the employed CWRU data set) are used to simulate the loading changes. For example, when the load L begins at 3 and decreases by ∆L of −3, the change can be expressed as 3 → 3 + (−3). In total, we have 12 combinations of simulated working load changes, as given in Table 1. In each combination, the vibration signals collected from drive end and fan end are used for analysis, respectively; thus, 24 signals are used for change detection. Table 1. Simulated working load change of L → L + ∆L hp. Figure 3 shows an example of change detection for a working load change from 1 hp to 2 hp, where the original signal was collected from drive end with 1:10 down sampling and the change is labeled by a human instructor. From Figure 3a-d, we show the detection result by our method, ARIMA model, Kurtosis and RMS. It can be seen that our method can successfully detect the change cycle with a slight time delay. As seen in Figure 3a, the anomaly scores are relatively small and are distributed irregularly when the motor load is set around 1 hp. Then, the occurrence of loading status change causes an abrupt increase of the anomaly score which can be successfully detected through the 3σ criterion based hypothesis testing; however, the method using the ARIMA model did not detect any alarm because the anomaly score is always in the confidence area, as shown in Figure 3b; the detected change cycle using the Kurtosis and RMS are much earlier than the actual change cycle, which belong to false detection as seen in Figure 3c,d; the main reason is that the proposed method has better power to detect weak changes than the other three methods.   Table 2 gives the comparison results. In the table, the N/3 represents that the number of samples N in which change cycle is correctly detected in the three samples, N ∈ {1, 2, 3}. It can be clearly observed that our method achieves the best performance and the detection precision is 100%. In other words, whether the vibration signal is collected from the driven or the fan end, the proposed method can detect all the state changes precisely without generating any false alarms.

Case Study II: Bearing Health Status Monitoring
This section considers two real-engineering scenarios of bearing health status: a gradually degenerated process and a sharply degenerated process. Consequently, this case study is employed for this investigation.
The testing data set comes from the publicly available PRONOSTIA platform [34]. This system is designed to test and validate bearings fault detection, diagnostic, and prognostic approaches. The main purpose of PRONOSTIA platform is to provide real experimental data that characterize the degradation of ball bearings along their whole operational life (until their total failure). In order to reduce the bearing's life duration, a radial force is set to the bearing's maximum dynamic load of 4 kN and applied on the tested bearings. The force is generated by a cylinder pressure, and the pressure is delivered through a pressure regulator.
During the tests, accelerometers are fixed on the outer race of the bearing, and vibration signals are captured. The sampling rate is 25.6 kHz, the length of every sampled acceleration waveform is 2560 data points, i.e., 0.1 s, whereas recordings are repeated every 10 s. The vibration signals are transmitted into a PC for data visualization and storage through a National Instruments data acquisition (DAQ) card. Based on experience, it is considered as the end of life when the vibration level is greater than 20 g.
The experiment is carried out under three different operating conditions. In our case study, four representative bearing signals during the first operating condition are used for bearing early failure detection. The operating condition is under 1800 r/min (resolution per minute) with a 4000 N load. Figure 4 shows detection results of the four bearings with different methods. From the original data of each Figure 4a,b, we can see that the vibration amplitude of bearing 1 and bearing 2 have gradually increasing trends, which indicates that the failure gets severe gradually. On the contrary, the vibration amplitude of bearing 3 and bearing 4 increases sharply at the end of lifetime as seen in the original data of each Figure 4c,d, which means a quick degradation processes. Table 3 gives the specific comparison results accordingly, where the number in the table represents the first early failure alarm point, and "N/A" represents false detection status change because there are too many alarm points and the first alarm point is far earlier than the early failure point. Thus, it can be clearly seen that kurtosis fails to distinguish the normal and the abnormal states, which can be verified by subfigure (3)s in all subfigures a, b, c and d in Figure 4. In addition, for bearing 3, as shown in subfigure (c), ARIMA and RMS have false alarms in a normal state. To summarize, the detection performance order is: our method > RMS > ARIMA > Kurtosis, revealing that our method has a superior ability compared to the others.

Case Study III: Speed Condition Monitoring
We used the experimental setup as shown in Figure 5 to collect testing data. In this setup, an accelerometer sensor mounted on the gearbox was used to acquire vibration signals with a sampling frequency of 1000 Hz, and then the collected signal data were transmitted to PC. During the data collection, we first set the initial speed of motor at ν, and then changed with an interval of ∆ν to simulate a speed status change (i.e., ν → ν + ∆ν). The testing values of ν and ∆ν were given as follows: • ν {250, 300, 350} rpm and ∆ν {50, 100, 150, 200, 250} rpm.
Therefore, we can obtain 15 parameter combinations with different speed condition changes as seen in Table 4. For each combination, we collected 10 data samples to form testing data set. Thus, we have a total number of 150 data sequences in our testing dataset.

Motor Accelerometer
Gear box Load Figure 5. Experimental setup. Figure 6 shows an example of change detection for a speed change from 350 rpm to 400 rpm, where the change is labeled by a human instructor. From Figure 6a-d, we show the detection result by our method, ARIMA model, Kurtosis and RMS. It can be seen that our method can successfully detect the change cycle. As seen in Figure 6a, the anomaly sores show a remarkable increase at the change time such that the change was detected successfully. However, the method using the ARIMA model did not detect any alarm because the anomaly score is always in the confidence area, as shown in Figure 6b; the detected change cycle using the Kurtosis and RMS is much earlier than the actual change cycle, which is attributed to false detection, as seen in Figure 6c,d. Table 5 gives the comparison results of the different methods under different settings of initial speed ν rpm. We can see that the proposed method also achieves a perfect detection performance on speed status change detection.

Results Summary and Discussion
To summarize, from the above experiments, we notice that the proposed differential equation-based prediction model has prior and remarkable detection performances in three different real applications, especially in weak state change. However, in the proposed method, we assume that the collected condition monitoring signal is periodically stationary. Although this assumption is reasonable and widely used in the literature [7,18,41], we have also noted that the actually collected signal in some cases may vary in periods, i.e., cycle non-stationary. With the purpose of extending the method to these cases, we would employ a preprocessing, e.g., angle re-sampling [55] and phase estimation using dynamic time warping [56,57], in order to suppress such temporal non-stationary before applying it to monitoring non-stationary condition signals.

Conclusions
In this paper, we have proposed a new method for monitoring and diagnosing the running condition of considered rotating machines. The base of the method is an early detection of condition changes based on the DE model that can realize one-step-ahead prediction of machine status. Together with this model, an analysis of continuously monitoring condition signal by means of sum residual error is presented to inspect/diagnose whether an abnormal condition change occurs or not based on a null hypothesis testing. The method was evaluated with three representative condition monitoring applications: external load status monitoring, bearing health status monitoring and speed condition monitoring. The method was also compared with those benchmark methods reported in the literature. The results demonstrated significant improvements of the proposed method over others, indicating its superiority and great potential in real engineering applications.
In the future, we will optimize the method to achieve a higher computational efficiency, and apply the method in the workshop for practical usages.
Author Contributions: X.W. conceived of the presented idea and developed the framework, performed the experiments and wrote the paper; G.C. contributed to the experiment; G.L. conceived the structure, provided guidance and modified the manuscript; Z.L. provided important guidance for this paper; and P.Y. gave a detailed revision. All authors have read and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.