Condition Monitoring of Industrial Equipment Based on Multi-Variables State Estimate Technique

: Unexpected failures commonly occur in industrial equipment, and condition monitoring could signiﬁcantly improve the e ﬃ ciency of maintenance and failure of early alarm. A condition monitoring method using multi-variables state estimate technique (MSET) is proposed, and an improved multi-variables memory matrix construction method is employed, furthermore, an analysis of comprehensive similarity index that considering variable weights is accomplished, and incipient failure alarm thresholds are determined, which lead to e ﬀ ective early detection of failure. The method proposed in this paper is validated using actual data for blower fan in a thermal power plant, and the simulation and comparison results are discussed. The veriﬁcation results reveal that the proposed method is e ﬀ ective for failure monitoring modeling and achieve a superior accuracy, and incipient failure could be accurately detected.


Introduction
Industrial equipment is usually complicated and features multiple failure modes. As a large-scale rotate machine, once failure occurs, a long-time repairing, even shutdown, is required, resulting in a significant economic loss. On the other hand, industrial field characterized harsh operational environment, equipment would experience a high failure rate, therefore, accurate condition monitoring for industrial equipment is necessary, and plays a pivotal role in condition-based maintenance, which would be more beneficial than corrective and preventive maintenance [1,2].
Tavner [3] concluded the general condition monitoring for industrial equipment, including temperature, chemical emission, vibration, electrical current, and electrical discharging. Furthermore, the paper illustrated the statistics distribution of failure modes in industrial equipment in the past decades, according to the research, rotation related failure occupies nearly a half (44%), higher than stator related (35%) and bearing related failure (21%), and according to previous work, unexpected vibration often occurs before serious failure events [4,5].
As a typical condition monitoring parameter, vibration is chosen for example and investigated in our paper. The aim of condition monitoring is inferring the condition of equipment rapidly, and giving a clear indication for incipient failure modes using minimum measurements. In industrial fields, vibrations are directly measured, for example, vibrations data are recorded in supervisory information system (SIS) in thermal power plants. In normal operating conditions, the vibrations of different directions are controlled in a certain range. The unexpected increase of vibrations corresponds to over-load, poor lubrication, bearing loosening, rotor unbalance, or shaft bending [6]. Therefore, the actual industrial data are applied for validation, and actual failure data for blower fan have been used for revealing the effectiveness of proposed method. Finally, Section 6 concludes the work.

Introduction of Research Object
A typical industrial blower fan is discussed in our research. It is manufactured by TLT Cooperation, equipped in a thermal power plant in northern China, used for conveying air into furnace of boiler for combustion. The blades of blower fan are integrated with driving motor by shaft, and the shaft is sustained by bearings. While the motor is running, the blades rotate and air is injected into the furnace of the boiler. The structure for blower fan is shown in Figure 1, and the actual layout and location of measuring points for investigated blower fan are also exhibited. data for blower fan have been used for revealing the effectiveness of proposed method. Finally, Section 6 concludes the work.

Introduction of Research Object
A typical industrial blower fan is discussed in our research. It is manufactured by TLT Cooperation, equipped in a thermal power plant in northern China, used for conveying air into furnace of boiler for combustion. The blades of blower fan are integrated with driving motor by shaft, and the shaft is sustained by bearings. While the motor is running, the blades rotate and air is injected into the furnace of the boiler. The structure for blower fan is shown in Figure 1, and the actual layout and location of measuring points for investigated blower fan are also exhibited. The Supervisory Information System (SIS) for the power plant records the parameters for blower fan in a 20-s sampling solution. In each sample, related parameters include: time stamp, current air volume, load of unit, vibration of bearings, temperature of bearings, current of driving motor, operation frequency of driving motor, and power of driving motor. The parameters are shown in Table 1. Particularly, the temperature values, including temperature of bearings and end-side bearings, are measured using platinum resistance temperature detector. Before variables selection for condition monitoring, the 2 out of 3 actions are accomplished, the absolute errors between each of the two temperature values are compared, and the average of temperatures whose absolute error is smallest is regard as the actual value of measurement, which would improve the reliability of the measurement.
Except the parameters, the actual vibration alarms were also recorded in Supervisory Information System; while the vibration values exceed the thresholds setting ahead, the time stamp is recorded. The 7-day normal condition data are available for modeling. Meanwhile, actual failure alarms are acquired in our investigated blower fan, which is used for incipient failure detecting validation in our papers.  The Supervisory Information System (SIS) for the power plant records the parameters for blower fan in a 20-s sampling solution. In each sample, related parameters include: time stamp, current air volume, load of unit, vibration of bearings, temperature of bearings, current of driving motor, operation frequency of driving motor, and power of driving motor. The parameters are shown in Table 1. Particularly, the temperature values, including temperature of bearings and end-side bearings, are measured using platinum resistance temperature detector. Before variables selection for condition monitoring, the 2 out of 3 actions are accomplished, the absolute errors between each of the two temperature values are compared, and the average of temperatures whose absolute error is smallest is regard as the actual value of measurement, which would improve the reliability of the measurement.
Except the parameters, the actual vibration alarms were also recorded in Supervisory Information System; while the vibration values exceed the thresholds setting ahead, the time stamp is recorded. The 7-day normal condition data are available for modeling. Meanwhile, actual failure alarms are acquired in our investigated blower fan, which is used for incipient failure detecting validation in our papers.

Multi-Variable State Estimate Technique (MSET)
By applying MSET, considerable historical data for industrial equipment are used for constructing the normal state envelope space, which is noted as memory matrix D. While the industrial equipment operating, the measurement parameters are observed real-time, and the related n parameters construct a vector noted as X obs . In the time sample j, the observation vector is written as, It should be pointed out that for different equipment condition monitoring, only the parameters and thresholds need to be determined if reusing the proposed method, and the complete preprocess of the case study in our paper is posed in Section 4.
D is an m × n dimension matrix, which means m normal samples are strictly selected and contained in D at last. Comprehensive operating conditions, including high load, low load, change load, and other operating conditions should be contained. Memory matrix D could be written as, The vectors included in the memory matrix D construct a nonlinear normal state envelope space. After the memory matrix D constructed, dynamic state estimation could be accomplished by applying it.
The construction of estimation vector for the same time sample X est is constructed by the vectors in memory matrix D, the different of X est and X obs is applied for condition monitoring. X est could be written as, The weight matrix W stands for similarity between observation vector and components in matrix D. The optimized weight matrix W could be obtained by minimizing the error between X est and X obs , which is, e = X est − X obs (4) According to the least square algorithm, for minimizing e 2 let, where, D ij is the element for memory matrix in i line j column. Solving the partial derivative Equation for w i in S(w), which is, That is, The Equation (7) could be expressed as matrix form as, The weight matrix could be solved, that is, In particular, if components in memory matrix D are linear correlative, the matrix D T · D would be irreversible, which leads the (9) unsolvable. In order to solve the problem properly, there are many alternative operator; Argonne laboratory has also developed a unique algorithm for nonlinear calculation [15]. A nonlinear operator ⊗ is introduced, and the Equation (9) could be written as, Concise nonlinear operator that calculates the Euclidean distance arithmetic is chosen in our paper, that is, The operator we employed is easy to understand physically, and the similarity of two vectors is represented by spatial distance. If the observation vector is similar with the vectors in memory matrix D, the corresponding spatial distances would be small, and while there is a distinguishable different between observation vector and memory matrix, the distance would become large. Finally, the estimation vector could be written as, The memory matrix D established a normal nonlinear normal state envelope space, while the input observation vector X obs is located in the normal space and similar with a vector in D, the corresponding vector X est obtained by Equation (12) is rather accurate. However, if an incipient failure occurs, the characteristic of X obs would change, and the observation vector X obs would exceed the boundary of envelope space, as a result, X obs could not be represented by using the vectors included in memory matrix D, and the current observation vector X obs and the estimate vector X est would be very different.
In conclusion, the incipient failure could be detected by analyzing similarity between X obs and X est .

Data Preprocess and Variables Selection
Data preprocess is an important stage in the data-driven modeling. Normally, missing values and outliers, should be removed or specially dealt with. In our research, the outliers and missing values are replaced by the average of the adjacent samples to guarantee the data continuity. On the other hand, the dimensional difference of input and output variables is eliminated by data standardization using Equation (13), which transforms different variables to same scale.
Before constructing the memory matrix, the variables need to be selected. Vibration variables are target vectors, and the relationship between the other variables is listed in Table 1 and target vectors need to be discussed. Grey relational analysis is an effective statistical method quantifying the relevance degrees, the "grey relational degree" of variables is calculated using consecutive data of variables, and the method is proved effective in many fields. Therefore, grey relational analysis is applied in variables selection in our paper.
Consecutive 1000 samples are used for obtaining grey relational degree. Vibration vectors are set as target sequences, and the rest 9 variables are set as comparing vectors. The grey relational degree sequence for variable could be calculated as, where, ρ is discrimination coefficient, and is set as 0.6 in our research. The average grey relational degrees for compared vector are calculated as, The results of grey relational analysis between the horizontal vibration, vertical vibration, and the compared vector are shown in Tables 2 and 3, respectively. According to the grey relational analysis results, except the vibration of horizontal and vertical itself, the temperature of end-side bearings, current air volume of unit, operation frequency of driving motor and current of driving motor (gray relational degree > 0.6), are chosen for vibration condition monitoring.

Vector Selection for Memory Matrix
While the MSET was first proposed, the D is constructed by the extreme vectors for each parameter, and the boarder of normal state envelope space is defined by the extreme vectors. In latest works, vectors of normal condition are added to memory matrix according to constant interval and access criteria, indeed, the method improved the memory matrix D and enhanced the ability for describing normal conditions. However, special normal condition vector might be omitted in sample process, which would seriously affect the prediction performance. Thus, a novel normal vectors selection method using self-adaption access criteria is proposed. After standardization, the deviation degree is reflected by the distance between normal vector with origin point in hyperplane. Because the vectors for memory matrix is selected in normal condition, the vectors are relatively intensive in the area that is close to the origin point in hyperplane, and in the area relatively far from the origin point, vectors are sparse and special normal condition often appears. The special normal condition vector might not be selected by using the constant interval and access criteria sampling method.
In order to select representative normal vector and retain some special normal vector, the sample access criteria change with probability density function, while the value of function is small, which means the vector is located in sparse area in hyperplane, the access criteria is relatively loose for introducing more special normal vectors. As shown in Figure 2, the maximum of probability density function is 0.62. The number of vector in memory matrix D is set as 200. The procedure for constructing memory matrix D is formulated in Figure 3.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 17 process, which would seriously affect the prediction performance. Thus, a novel normal vectors selection method using self-adaption access criteria is proposed. After standardization, the deviation degree is reflected by the distance between normal vector with origin point in hyperplane. Because the vectors for memory matrix is selected in normal condition, the vectors are relatively intensive in the area that is close to the origin point in hyperplane, and in the area relatively far from the origin point, vectors are sparse and special normal condition often appears. The special normal condition vector might not be selected by using the constant interval and access criteria sampling method. In order to select representative normal vector and retain some special normal vector, the sample access criteria change with probability density function, while the value of function is small, which means the vector is located in sparse area in hyperplane, the access criteria is relatively loose for introducing more special normal vectors. As shown in Figure 2, the maximum of probability density function is 0.62. The number of vector in memory matrix D is set as 200. The procedure for constructing memory matrix D is formulated in Figure 3   Where, k is an empirical parameter, and is set as 0.0015 in our simulation.

Improved Similarity Calculation
The residual of vibrations component in Xobs and Xest, including the average value for residual sequence, are usually used for judging whether the incipient failures occur in previous works. However, the information contained in the other variables, including temperature of end-side bearings, current air volume of unit, operation frequency of driving motor, and current of driving motor, would be neglected, but some vibrations failure might be first reflected by indirect parameters. Therefore, except the analysis the residual of vibrations directly, a comprehensive similarity calculation based on other components in observation vector is proposed, and the difference between components is considered.
The similarity could be formulated as, where Xobs,i and Xest,i are the component i in observation vector and estimation vector, respectively, and Ri is corresponding weight coefficient for i component in vectors, which could be obtained according to the result of grey relational analysis in Table 4. Where, k is an empirical parameter, and is set as 0.0015 in our simulation.

Improved Similarity Calculation
The residual of vibrations component in X obs and X est , including the average value for residual sequence, are usually used for judging whether the incipient failures occur in previous works. However, the information contained in the other variables, including temperature of end-side bearings, current air volume of unit, operation frequency of driving motor, and current of driving motor, would be neglected, but some vibrations failure might be first reflected by indirect parameters. Therefore, except the analysis the residual of vibrations directly, a comprehensive similarity calculation based on other components in observation vector is proposed, and the difference between components is considered.
The similarity could be formulated as, where X obs,i and X est,i are the component i in observation vector and estimation vector, respectively, and R i is corresponding weight coefficient for i component in vectors, which could be obtained according to the result of grey relational analysis in Table 4.
where n is the number of components that are selected in to memory matrix.

Residual Threshold
Residual threshold is determined for residuals between X obs and X est , and it is a direct vibration monitoring parameter, information contained in vibration component of observation vector is applied. In order to reflect the long-term trend of vibration monitoring and eliminate the influence of measurement errors, residuals sequence is discussed and sliding window is applied. The threshold of average value T aver is defined as, T aver = |c 1 · AVER| where AVER is the average residual in a sliding window. The current residual is compared with T aver . The coefficient c 1 is determined by operation experience and validation, in this paper, c 1 = 2.5.

Similarity Threshold
As a supplement for direct vibration monitoring, information contained in the other component in X obs and X est are applied to calculate a comprehensive similarity, which is significant for detecting faults reflecting in indirect variables. Furthermore, as shown in Figure 1 and Table 1, there are no redundancies for vibration measurement, therefore, similarity judgement is an effective way in situations where vibrating measurement fails, and it is significant for improving the accuracy of condition monitoring. The similarity threshold could be formulated as, where S(X obs , X est )min is the minimum of similarity in a sliding window, and coefficient c 2 is set as 1.1 in this paper.

Validation of Blower Fan MSET Modeling
The sample solution of data is 20 seconds, and the variables we discussed in Section 4 is available, which are vibration of horizontal x 1 , vibration of vertical x 2 , the temperature of end-side bearings x 3 , current air volume of unit x 4 , operation frequency of driving motor x 5 , current of driving motor x 6 , respectively. The typical distribution for normal data in 10/25/2015 is shown in Figure 4, the maximum of vibration in horizontal and vertical direction are 2.725 mm/s and 5.195 mm/s, respectively. In totally, 7 days normal samples were used for contrasting memory matrix D.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 17 After using the method mentioned in Section.4, memory matrix D is constructed, and 211 vectors are concluded, therein, 11 vectors are selected by extreme condition of 6 variables and remove duplication, and the rest 200 vectors are obtained by method mentioned in Figure 3. Consecutive normal samples, including steady and dynamic working condition data, are used for validating the accuracy of MEST model.

Steady Working Condition
In order to validate the effectiveness of the proposed method, traditional constant interval and access criteria method [18] are applied for comparison, which is shown as traditional method in the figures, and the size of memory matrix is also set as 200. It should be pointed out that in traditional method, the access criteria ε is a constant, and extreme condition vectors are also included in traditional method. The comparisons for proposed method and traditional method of horizontal and vertical vibration are shown in Figures 5 and 6, respectively, and the residual and similarity are compared.
The dataset for validation in Figures 5 and 6 is relative steady and chosen from a normal working condition for blower fan, named normal condition I.  After using the method mentioned in Section 4, memory matrix D is constructed, and 211 vectors are concluded, therein, 11 vectors are selected by extreme condition of 6 variables and remove duplication, and the rest 200 vectors are obtained by method mentioned in Figure 3. Consecutive normal samples, including steady and dynamic working condition data, are used for validating the accuracy of MEST model.

Steady Working Condition
In order to validate the effectiveness of the proposed method, traditional constant interval and access criteria method [18] are applied for comparison, which is shown as traditional method in the figures, and the size of memory matrix is also set as 200. It should be pointed out that in traditional method, the access criteria ε is a constant, and extreme condition vectors are also included in traditional method. The comparisons for proposed method and traditional method of horizontal and vertical vibration are shown in Figures 5 and 6, respectively, and the residual and similarity are compared.
The dataset for validation in Figures 5 and 6 is relative steady and chosen from a normal working condition for blower fan, named normal condition I. Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 17    As shown in Figures 5 and 6, the ranges for horizontal and vertical vibrations are 0.589 mm/s to 0.752 mm/s and 0.717 mm/s to 1.253 mm/s, respectively, which are normal vibration ranges. The vibrations in the dataset are steady, and in that condition, the residuals comparing results reveal the advantage of the proposed method, the horizontal vibration or vertical vibration. The average absolute values of residuals for the proposed method for horizontal and vertical vibrations are 0.003 mm/s and 8.73 × 10 −4 mm/s, respectively, that for traditional method are 0.0143 mm/s and 0.0238 mm/s, respectively. The similarity obtained according to Equation (17) also reflects the advantage of the proposed method, the similarity for the proposed method is above 0.97, including horizontal and vertical vibrations, that for most of samples are over 0.999, and the results for similarity comparisons prove the advantage of proposed method.
For comprehensively evaluating the ability of characterizing normal state envelope space, the root mean square error (RMSE) and MAPE (mean absolute percentage error) are applied for comparison, the calculating formula are, where y i andŷ i are observation and estimation values, respectively. The RMSE and MAPE comparison of proposed method and traditional method for normal condition I is shown in Table 4. In normal condition I, the RMSE and MAPE for proposed method of horizontal vibration x 1 are 0.0045 mm/s and 0.45%, respectively, that for traditional method are 0.0265 mm/s and 2.12%, respectively. For vertical vibration x 2 , the RMSE and MAPE for proposed method are 0.0017 mm/s and 0.11%, that for traditional method are 0.0177 mm/s and 2.79%, respectively. The results reveal that while in a relative steady normal condition, the proposed method achieves a better estimation performance.

Dynamic Working Condition
In order to further validate the effectiveness of the proposed method, the horizontal and vertical vibrations for blower fan in dynamic work condition, named condition II, are applied. The comparison for proposed method and traditional method in condition II is shown in Figures 7 and 8.    As shown in Figure 7, the horizontal vibration of test data set is change in the range of 1.144 mm/s to 1.767 mm/s, which is the normal vibrating range according to engineering handbook. Comparing the residual results, the proposed method achieves a higher accuracy, the residual is contained in −0.0243 mm/s to 0.0274 mm/s, which for traditional method is −0.0411 mm/s to 0.0756 mm/s.
It should be pointed out that the condition changed obviously at certain samples, such as sample A and sample B marked in the Figure, and the condition changes do not reflect in residuals obviously, however, the similarity reveals the condition changes accurately, at the sample A and sample B, the vibrate condition changes and similarity descends correspondingly. Owning to the calculation of the other variables in observation vector, information except direct vibration variables is contained in similarity calculation, which is significant for a comprehensive vibration condition monitoring.
As shown in Figure 8, the vertical vibration of test dataset is changed in the range of 1.573 mm/s to 3.273 mm/s, which is also a normal vibrating range, and the trend for vertical vibration is similar with horizontal vibration. The ranges for residual for proposed method and traditional method are −0.048 mm/s to 0.1434 mm/s, −0.5495 mm/s to 0.192 mm/s, respectively, and the results reveal that the proposed method achieve a higher estimation accuracy.
The RMSE and MAPE for condition II is shown in Table 5. For horizontal vibration x 1 , the RMSE and MAPE of the proposed method are 0.0077 mm/s and 0.36%, respectively, that for traditional method are 0.0189 mm/s and 0.85%, respectively. For vertical vibration x 2 , the RMSE and MAPE for the proposed method are 0.0293 mm/s and 0.94%, that for the traditional method are 0.1323 mm/s and 3.33%, respectively. The results reveal that while in dynamic normal condition, the proposed method keeps a superior estimation performance.

Incipient Failure Detection
According to the operation handbook, the vibration conditions often evaluated by measuring values directly, the vibration alarm is not triggered until the vibration measure values exceed the threshold, therefore, an early alarm could not be realized, and the condition judgment criteria is shown in Table 6, which means once the vibration measurements exceed 5.1 mm/s, the alarm occurs. The exceed-threshold events could be detected by an effective failure early detection, and the incipient failure thresholds for residual and similarity proposed in 4.4 are applied. As shown in Figure 9, the typical vertical vibration alarm events are discussed, and the process is exhibited, residual and similarity are shown correspondingly. There are two vibration alarms in the process, noted as alarm 1 and alarm 2, respectively. According to the data records, the timelines for the event are:  The size of sliding window is set as 10, and the effectiveness of residual threshold and similarity threshold are validated in the alarms. The thresholds, current value, and corresponding sample for alarms are compared in Table 7. In the alarm 1, the actual alarm is NO.124 sample, and the failure alarm is earlier detected in NO.120 sample and NO.121 sample by using residual threshold and similarity threshold, respectively, and the failure could be detected 1 minute ahead. Furthermore, in the alarm 2, the actual alarm is NO.276 sample, in comparison, the failure alarm is detected in NO.263 sample and NO.260 sample, respectively, and the failure is detected over 5 minutes ahead.
For a further discussion, the other variables in vector used for similarity calculation are investigated, and the results indicated that the deviations of the temperature of end-side bearings x3 The size of sliding window is set as 10, and the effectiveness of residual threshold and similarity threshold are validated in the alarms. The thresholds, current value, and corresponding sample for alarms are compared in Table 7. In the alarm 1, the actual alarm is NO.124 sample, and the failure alarm is earlier detected in NO.120 sample and NO.121 sample by using residual threshold and similarity threshold, respectively, and the failure could be detected 1 minute ahead. Furthermore, in the alarm 2, the actual alarm is NO.276 sample, in comparison, the failure alarm is detected in NO.263 sample and NO.260 sample, respectively, and the failure is detected over 5 minutes ahead.
For a further discussion, the other variables in vector used for similarity calculation are investigated, and the results indicated that the deviations of the temperature of end-side bearings x 3 and current of driving motor x 6 lead to the similarity descend in the early stage of failure, as shown in Figure 10.
The results indicate that failure information reflected by the other variables indicates the prognosis of failure.

Conclusions
An industrial equipment condition monitoring method based on MEST is proposed in this paper, and a novel method for construction of the memory matrix is put forward. The improved similarity that combines relative information in other variables is applied in monitoring. For different equipment and monitoring parameters, the method could be applied. The actual operational data for a blower fan in thermal power plant is used for validation, and the condition estimation performances are discussed. The actual data for failure event are used for validating the effectiveness of incipient failure detection employing proposed condition monitoring and early alarm method, the performance is compared with traditional exceed threshold method.
The method has advantages such as, it is concise and easy to understand conceptually, which should be attractive to industrial engineers. Irrespective of the steady or dynamic work condition, the proposed memory matrix construction method improves the accuracy of condition monitoring.

Conclusions
An industrial equipment condition monitoring method based on MEST is proposed in this paper, and a novel method for construction of the memory matrix is put forward. The improved similarity that combines relative information in other variables is applied in monitoring. For different equipment and monitoring parameters, the method could be applied. The actual operational data for a blower fan in thermal power plant is used for validation, and the condition estimation performances are discussed. The actual data for failure event are used for validating the effectiveness of incipient failure detection employing proposed condition monitoring and early alarm method, the performance is compared with traditional exceed threshold method.
The method has advantages such as, it is concise and easy to understand conceptually, which should be attractive to industrial engineers. Irrespective of the steady or dynamic work condition, the proposed memory matrix construction method improves the accuracy of condition monitoring. The proposed method using residual and similarity enhances the alarm timeliness, and proposed method could identify the incipient failure ahead.
Although the proposed method has been employed in a specific blower fan, the model can be further explored with suitable parameters selection and threshold determination. The method should be further proved for different industrial equipment, and further research is planned for taking the method into account in other industrial equipment.