Yanji Sun,Yanqiu Pan, *,Zhongliang Zhou,Xin Li

1 School of Chemical Engineering,Dalian University of Technology,Dalian 116024,China

2 Petro-CyberWorks Information Technology Co.,Ltd.,Beijing 100020,China

Keywords:Compressor Condition assessment Mathematical modeling Maximal information coefficient Dynamic deterioration degree Prediction

ABSTRACT A fuzzy comprehensive assessment method of running condition was constructed and applied to a large-scale centrifugal compressor set in a petrochemical corporation aiming at the monitoring and early warning of abnormal conditions in industry.The maximal information coefficient(MIC)correlation analysis of indexes was introduced to determine the independent indexes to be assessed,and the dynamic deterioration degree was calculated using the predicted independent indexes by the second-order Markov chain model.The fuzzy membership degree weighting method was employed to assess the running condition of all units in the set.Simple and direct radar chart was used to visualize condition assessment grades.Results showed that the proposed fuzzy comprehensive assessment method successfully assessed the running condition of the set.The constructed method achieved 10 min earlier alarm than the traditional threshold alarm occurred at the specific moment of 00:44 on April 7 of 2018.The method would provide a valuable tool and have a wide engineering application in ensuring the safety and reliability of industrial production.

1.Introduction

Large-scale centrifugal compressor set is widely used in petrochemical industries.It is characterized by complex structure,high degree of automation,strict process requirements and long maintenance time.Once a failure occurs,it may cause the interruption of whole production,and even lead to a great economic loss or casualty[1].Such a set usually consists of three parts:mechanical part,process part,and auxiliary part.The mechanical part includes a compressor and a driving machine.The auxiliary part comprises an oil circuit,a gas seal,and so forth.In general,vibration sensors are installed on some key components of the mechanical part in a compressor set,with which vibration signals can be monitored.Available researches mainly concentrated on fault diagnosis of the mechanical part based on the vibration signals[2-4].There is lack of reports on the assessment of overall running condition,including all the three parts,towards a large-scale centrifugal compressor set.In addition,the vibration signal analysis was often performed after an initial failure so an early warning and prediction was difficult.Therefore,it is necessary to construct a comprehensive monitoring and assessment method of overall running condition for the large-scale centrifugal compressor set.

Recently,the overall running condition assessment for large-scale sets based on big data has become a hot spot of researches[5-7].Some approaches,such as knowledge-based expert system[8],Gauss mixture model[9,10],neural network-based intelligent method[11-13],and fuzzy comprehensive assessment[14-16],have been gradually applied to the monitoring and assessment of running condition with multi-index considered.Godwin and Matthews developed an expert system based on the prior knowledge with 14 built-in interpretable rules.This system can be used for condition assessment and initial fault diagnosis of supervisory control and data acquisition(SCADA)for the wind power set[8].However,this method exhibits its limitations for the overall assessment due to the difficulty of system update and maintenance.Dong et al.established a condition assessment model based on the multi-state feature fusion of Gauss mixture model,which employed a health degradation index to assess the condition of the wind power set[9].The method is only suitable for the case that the assessment indexes have to follow the Gaussian mixture distribution along with a complex and massive calculation.Sun et al.developed a running behavior model of the wind power set based on the neural network and stochastic process theories,and assessed the performance of 9 sets through the SCADA data analyses[11].This kind of assessment needs lots of historical data as training samples to improve the computation accuracy.Otherwise,output results are sometimes difficult to be explained,which is limitedly used in practice.The fuzzy comprehensive assessment does not need to rely too much on prior data.It can effectively deal with the fuzziness and randomness of multi-index,which has been gradually accepted for the condition assessment in recent years.Li et al.constructed their improved combination weighted fuzzy synthesis model and cloud model to optimize the maintenance strategy of wind power sets and successfully improved the operational reliability.They found that the method can reflect the real running condition better than the others[14,15].Compared with wind power sets,a large-scale centrifugal compressor set in petrochemical industries is more complex in terms of structure,process and index.Therefore,the fuzzy comprehensive assessment is pretty suitable for and enough to be used for the condition assessment of compressor sets due to its own characteristics.

Three aspects need to be paid more attention in the fuzzy comprehensive assessment of running condition for a large-scale centrifugal compressor set.First,the method itself cannot solve the problem of information overlap due to correlations among indexes[17].These indexes have to be pretreated to identify and eliminate significantly correlated ones in order to ensure the accuracy of results.This kind of pretreatment cannot be found so far in literatures.A large-scale centrifugal compressor set involves numerous indexes,among which are linearly or non-linearly correlated.The commonly used Pearson correlation coefficient and least square regression coefficient methods are only suitable for linear relations,but not for non-linear relations.The Maximal Information Coefficient(MIC)method was developed on the basis of the mutual information theory[18].It can measure the dependencies among variables extensively,including both the linear and non-linear relations.It can also be effective for the non-functional dependency that cannot be expressed by a single function[19,20].Therefore,it is essential to introduce the MIC method for correlation analyses of indexes before carrying out the fuzzy comprehensive assessment to improve the accuracy of results.Next,the concept of deterioration degree has been used to describe the deviation of indexes from their normal values in order to meet the needs of different indexes to achieve an overall assessment[14].However,the deterioration degree only accounts for current values without considering the influence of changing trend on the running condition.Xiao et al.improved the traditional fuzzy comprehensive assessment method by combining a fuzzy trend prediction strategy with the deterioration degree of condition parameters to reduce the randomness and uncertainty of indexes in running conditions of sets[21].Although the above-mentioned two works conducted the twolevel fuzzy comprehensive assessments,they all ignored the correlation analyses of units,leading to the decreased accuracies of assessment results.Zheng et al.introduced the dynamic deterioration degree combined with the first-order Markov chain model to predict the changing trend of indexes so as to dynamically describe running conditions of wind power sets[16].They replaced the traditional membership function by the normal cloud model,and integrated it into the fuzzy comprehensive assessment to form a real-time assessment method.Shamshad et al.performed the first and second order Markov chain models for synthetic generation of wind speed time series[22].Carpinone et al.proposed a Markov chain model based probabilistic prediction method for an ultra-short-term prediction of wind power sets[23].They revealed that the second-order Markov chain model can improve the prediction performance and decrease the prediction error.Other researchers have also concluded that the Markov chain model is simple,fast and accurate without assumptions of prediction errors in the ultra-short-term prediction[24,25].These beneficial features fit well with such requirements as a short time interval of monitoring and a large amount of data,as well as a precise and real-time trend prediction.Finally,three levels,including target level,unit level and index level,in the fuzzy comprehensive assessment have been adopted for a large-scale set with numerous indexes[14,16,21].The eventual result of a multi-level fuzzy comprehensive assessment is presented by an assessment vector of the target level.However,a multi-level assessment may result in the concealment of some indexes through the transfer of multi-level weight factors[21].Meanwhile,a large number of indexes in a level would lead to a huge calculation capacity in correlation analyses between each two indexes.Therefore,conducting a correlation analysis between indexes within a unit and a fuzzy comprehensive assessment only for unit level would be a more preferable way for large-scale sets.This method also conforms to the daily monitoring habits of industrial operators.

The objectives of the work are to construct an improved fuzzy comprehensive assessment method and to assess a large-scale centrifugal compressor set in a petrochemical corporation.Three improvements were implemented:the MIC method for the correlation analyses between indexes,the dynamic deterioration degree based on the second-order Markov chain model for increasing the prediction performance of indexes,and the visualized radar graph for the improved display strategy of assessment results.The constructed assessment method was then applied to a case study and compared with the traditional threshold warning method.The results obtained can help ones to early discover hidden troubles and reasonably arrange maintenance plans of large-scale centrifugal compressor sets.The method provides a reliable tool in the comprehensive assessment of running condition for similar large-scale sets.

2.Construction of Fuzzy Comprehensive Assessment Method

A fuzzy comprehensive assessment method conventionally consists of four procedures:establishing an index system and an assessment grade,selecting an appropriate membership function to calculate membership matrix,L,determining the weight factor vector,W of indexes,and implementing an intended condition assessment and evaluating results.The present work introduced the correlation analysis and dynamic deterioration degree calculation of indexes after the establishment of index system and assessment grade,and then used the dynamic deterioration degrees of independent indexes as the input of the membership function for subsequent calculation.The detailed flow chart is shown in Fig.1.

Fig.1.Flow chart of centrifugal compressor unit condition assessment.

2.1.Establishment of index system and assessment grade

According to the function distributions of a set to be assessed,the set is divided into several assessment units which are expressed as{R1,R2,R3,…,Rn}where the n is the total number of units.Then,the indexes that can represent conditions of the unit are selected.For any unit Ri,they are expressed as{Ri1,Ri2,Ri3,…,Rim}where the m is the total number of indexes.There are z condition assessment grades for a unit,which are recorded as{V1,V2,V3,…,Vz}.

2.2.Correlation analysis of indexes

Relationships among numerous indexes can be linear,non-linear and even non-functional.A reliable correlation analysis method is needed to identify redundant indexes,and retain independent indexes for subsequent analysis and calculation.This pretreatment of indexes can reduce the capacity and complexity of computation and improve the accuracy of assessment.The MIC method can be sufficiently used to measure the correlation degree between indexes so as to determine independent ones[20].

The basic idea of the MIC is to segment points in the plane scatter graph by drawing a grid when two indexes are mutually correlated,for example,indexes Ri1and Ri2in the unit Ri.By gradually increasing the grid resolution and calculating the maximum mutual information value under each resolution,all the maximum values are standardized to obtain calculation results with the following equation:

where ua,bis the normalized maximum mutual information value of the segmented grid by the resolution of a multiplied b,and the a and b represent the number of segmentation points on the abscissa and ordinate axes in the scatter graph,respectively.MIC(Ri1,Ri2)is the percentage of information that can be interpreted by the index Ri1in the index Ri2,ranging from 0 to 1.The closer the MIC value to 1,the stronger the correlation between two indexes.When the two indexes are not correlated,MIC=0;when they are definitely correlated,MIC=1.

2.3.Dynamic degradation degree calculation of independent indexes

The predicted values of indexes are used as the numerical value of the dynamic deterioration degree analysis to reflect the unit condition being assessed.According to the correlation analysis results of indexes,the independent indexes are normalized as the dynamic deterioration degrees.

There are three types of calculation formulas,the smaller value better type,the middle type and the larger value better type[14,17].The smaller value better type and the middle type are selected as given respectively below:

where g(xp)is the dynamic deterioration degree which varies between 0,the best value and 1,the worst value.The xminand xmaxare the minimum and maximum values of independent indexes in the normal operating state,the xaand xbare the upper and lower values in the optimal operating state,and the xpis the predicted index value that is calculated by the second-order Markov method[23].The procedures are described as follows:

(1)Divide the interval of index variation into q subintervals which are expressed as{S1,S2,S3,…,Sq}.

(2)Count the transfer frequency of each subinterval to obtain the transfer frequency matrix P=(pijk)where the pijkdenotes the frequency of the subinterval transferred from Siat the last moment,i to Sjat the current moment,j,and then to Skat the next moment for an index.

(3)Calculate the subinterval transfer probability matrix PP=(ppijk)using Eq.(4)where the ppijkdenotes the probability of the subinterval transferred from Sito Sjand then to Skof the index.

(4)Calculate the index variation,Δxt+Δtat the next moment,t+Δt according to Eq.(5),which is the sum of the product of the median of the subinterval Skat the next moment and the corresponding transfer probability ppijkaccording to the variation,Δxt-Δtin the subinterval Siat the last moment and the variation,Δxtin the subinterval Sjat the current moment.

where Sk,upis the upper bound of subinterval Skand Sk,downthe lower bound of Sk.

(5)Obtain the predicted value xpof the index by adding the current value with the variation Δxt+Δt.

2.4.Membership matrix calculation of the unit

According to the dynamic deterioration degree of each independent index in each unit,the appropriate membership function is used to calculate the membership corresponding to each condition assessment grade to obtain the membership matrix,L for each unit.The selection of membership function should reasonably cover the whole range of dynamic deterioration.The membership function of ridge distribution has the advantages of simple shape and small error compared with other complex membership functions[21].It can therefore be selected to calculate the membership matrix for subsequent assessment.Detailed procedures can be seen in the case study.

2.5.Weight factor vector calculation of the unit

During the calculation of weight factor vector,expert evaluation method and analytic hierarchy process(AHP)are often used in practical applications.The AHP still relies on expert knowledge and experience,and then uses mathematical methods to eliminate subjective components as far as possible.It also needs a consistency test of the comparison matrix to guarantee the rationality of weight factors,as well as the objectivity and authenticity of a fuzzy comprehensive assessment result[21].

The AHP method was used to determine the weight factor of each independent index.The first step of the AHP is to establish a hierarchical structure model,and then construct a pairwise comparison matrix.The second step is to calculate the weight factor vector and conduct its consistency test,and then calculate the combined weight vector and conduct its combined consistency test[26].The weight factor can be constant and variable.Compared with the wind power set,the large-scale centrifugal compressor set runs more stably,and the correlation between indexes basically remains unchanged.Therefore,the constant weight is in line with the characteristics of the compressor set.The weight vector is denoted by W for a unit.Detailed procedures can be seen in the case study.

2.6.Fuzzy comprehensive assessment of running condition

The weighted membership vector of the ith unit is the dot product of the weight factor vector,Wiand the membership matrix,Lias expressed below[27]:

where the m is the membership value of the unit under different condition grades.

The rule of the“maximum membership”is adopted as the assessment criterion for each unit[14].A simple and direct radar graph is employed to visualize the eventual result of the fuzzy comprehensive assessment for the set.

3.Case Study

With the constructed fuzzy comprehensive assessment method,the present work performed an assessment of running condition for a BCLtype centrifugal compressor set in a petrochemical corporation as a case study.

3.1.Index system and assessment grade

The index system of the centrifugal compressor set was constructed based on the on-site operating parameters provided by the real-time database(RTDB)and the S8000 on-line condition monitoring system.The data were collected ranging from January 1,2018 to April 30,2018 with an interval of 1 min for each record.

According to the unit functions,the centrifugal compressor set was divided into five units,a compressor unit,a turbine unit,a process unit,an oil circuit unit and a gas seal unit.50 main indexes representing its running condition were selected in total for the set as listed in Table 1.Measuring points for each index are located at the P&ID diagrams of the three parts as shown in Figs.2-4.

Based on the actual management requirements of the centrifugal compressor set,the running condition was divided into four grades:V1for excellent,V2for good,V3for medium and V4for poor.There is a condition warning under the poor grade.

3.2.Correlation analysis result of indexes

The correlation analysis of indexes based on the MIC method was conducted towards each unit.Due to the similarity of the calculation process,only the analysis of the gas seal unit is exhibited here as anexample.Detailed procedures of the correlation analysis for the two indexes are described as follows.In the open access software Python 3.7.0,load the Numpy numerical calculation package and the Mine package,and input the monitored values of two indexes.By running the MIC function,i.e.,Eq.(1)in the Mine package,the maximum mutual information values of all segmentation modes in different grids were calculated and outputted.For the two indexes,FI24674 and FI24675 within the gas seal unit,for example,the MIC value was 0.616.The same procedures were repeated for other pairs of the indexes,and the results were summarized in Fig.5,which was a symmetrical matrix.

Table 1 Indexes corresponding to five units of the set to be assessed

During the correlation analysis,if the MIC value is more than 0.5,the two indexes are considered highly correlated.Otherwise,they are less correlated[20].The PDIC24671(a pressure difference of the unit)was found to be poorly related to the other indexes from the row or column values as shown in Fig.5.This indicates that the PDIC24671 is an independent index in the unit,and it should therefore exist in each combination.The next step was to seek combinations of the other five indexes.The FI24674 was correlated with all other indexes except for the PDIC24671 from column one.Hence,the FI24674 and the PDIC24671 were considered as combination C1.The FI24675 was not correlated with the PDI24670 and the PDI24672 from column two in Fig.5.However,the latter two indexes were correlated with each other.Therefore either the PDI24670 or the PDI24672 can be combined with the FI24675 together with the PDIC24671 as combinations C2 and C3.Combinations C4 and C5 were obtained similar to combination C1.Table 2 lists all five combinations representing the condition of the gas seal unit.According to the rule of less independent indexes that can sufficiently represent the running condition of the unit,along with operator experiences,combination A containing the FI24674(the flow rate of exhaust)was selected as the independent index system for further assessment.

In the same way,the independent indexes of the other units can be obtained.Twenty-five independent indexes of all the units are listed in Table 3.

3.3.Dynamic degradation degree calculation result of independent indexes

The second-order Markov chain model was employed to calculate the predicted value and dynamic deterioration degree of each independent index.The index XI24666 of the process unit was used as an example.According to the collected operating data of indexes,the XI24666 was ranging between-2.422 and 2.527,which was denoted as[a,b]=[-2.422,2.527]of the index variation sequence.Based on the central limit theorem,the average value distribution of the sample can be approximated as the normal distribution when the sample size is sufficiently large(≥30)[28].Since the variations of the selected independent indexes all satisfy the normal distribution,the central limit theorem was therefore employed for subinterval division,i.e.,taking the sample mean value as the center and the standard deviation as the standard for dividing[29].Procedures of the dynamic degradation degree calculation are briefly described as follows[30].The arithmetic mean value of the XI24666 wasand the standard deviation(S=was S=0.140.Five subintervals of the XI24666 were divided with the formulas in the second column of Table 4.

Based on the variation sequence and the subintervals in Table 4,the target subinterval at which the variation value was located was able to be determined.The subinterval transfer frequency matrix P=(pijk)was obtained by counting the subinterval transfer frequencies with i,j and k from 1 to 5.Table 5 lists the results for i=1 and j=1 as a representative.

Fig.2.P&ID diagram of mechanical part(compressor and turbine units).

The subinterval transfer probability matrix PP=(ppijk)was calculated with Eq.(4)based on Table 5.The results are shown in Table 6 for the case of i=1 and j=1.

The index variation,Δxt+Δtat the next moment,t+Δt was calculated with Eq.(5)based on the subinterval transfer probability matrix PP=(ppijk).The predicted value,xpwas obtained by adding the current value with the variation,Δxt+Δt,and then substituted into a calculation formula,either Eq.(2)or Eq.(3),of the dynamic deterioration degree to obtain the dynamic deterioration degree,g(xp)of the indexes.According to index characteristics,temperature,vibration,displacement and pressure difference of PDI24601 adopt to the smaller value better type,and the others to the middle type.The index XI24666 represents the throttle openness in the process unit,which definitely belongs to the middle type.Its xpand the g(xp)with Eq.(3)at moments t1-6are listed in Table 7.

At the moment of 00:34 on April 7 of 2018,the dynamic deterioration degrees of each independent index for the five units are listed in Table 8.

The dynamic deterioration degrees of predicted indexes were then substituted to the membership functions for the calculation of fuzzy membership matrix as seen in the next subsection.

3.4.Calculation results of fuzzy membership matrix for all units

The ridge distribution membership function was used to calculate the dynamic deterioration degrees of the indexes and the fuzzy membership relationship of each condition grade in each unit[21].The membership functions are written below according to the 4 condition grades of the unit.

Fig.3.P&ID diagram of process unit.

Fig.4.P&ID diagram of auxiliary part(oil circuit and gas seal units).

Fig.5.MIC values of indexes in the gas seal unit.

Table 2 Independent index combinations of the gas seal unit

Table 3 Independent indexes of the five units

The dynamic deterioration degree,g(xp)of the index for each unit was substituted into the above four membership functions(Eqs.(7a)-(7d)).The membership degree matrix,Liof the dynamic deterioration degree for the ith unit with total independent index number of h was therefore obtained as follows:

By substituting the dynamic deterioration degrees in Table 8 into Eqs.(7a)-(7d)and(8),the membership degree matrices corresponding to the 5 units were obtained respectively as below:

3.5.Calculation results of weight factor vector for all units

The AHP was used to calculate the weight factor vector for all units[31].Only the calculation procedures of the process unit are given here as an example to avoid repeated descriptions.

Table 4 Subintervals of the index XI24666

Table 5 Subinterval transfer frequency matrix P=(p11k)

Table 6 Subinterval transition probability matrix PP=(pp11k)

A pairwise comparison matrix A=(aij)4×4was first constructed according to the importance intensity of the independent indexes,where aijrepresents the importance ratio between the two indexes,and aij=1/ajiwith i,j=1,…,4.The importance intensity was regulated according to the rating scale method of AHP[32],as well as the engineering design and maintenance records of the unit along with the fault statistics of similar sets[33].The rule of the rating scale method is that each of the judgments towards the importance intensity is assigned a number on a scale.The importance intensity can be“equal importance”,“somewhat more important”,“much more important”,“very much more important”and“absolutely more important”,which corresponds to 1,3,5,7 and 9,respectively.Among the four independent indexes,FIC24033,PIC24226,TI24140 and XI24666,the former three were equally important while the XI24666 was somewhat more important than the others.Therefore,the ratios of two importance intensities among the FIC24033,PIC24226 and TI24140 were all 1,and the ratios of the XI24666 to each of the other three ones were all 3.The comparative matrix A was then constructed as

Secondly,the consistency of the matrix A was tested with the following equation.

where the CR is the consistency ratio,the RI is the averaged random consistency value that was 0.9 from the matrix A,and the CI is the consistency value calculated from the following equation.

where the λmaxis the maximum eigenvalue of the matrix A.

Table 7 Predicted values and dynamic deterioration degrees of index XI24666

Table 8 Dynamic deterioration degrees of each index for all units

If the consistency ratio value is less than 0.1,the matrix A is consistent.The matrix A passed the consistency test in the present work.

Finally,the weight factors,wiof the four indexes were obtained with the following equation.

The results were w1-3=0.167 and w4=0.499.The weight factor vectors of the other units are listed in Table 9.

3.6.Fuzzy comprehensive assessment result of all units

With Eq.(6),the weighted membership vectors of the five units were calculated as(1,0,0,0)for the compressor as M1,the turbine as M2,the oil circuit as M4,and the gas seal as M5,and(0.454,0.046,0,0.500)for the process as M3.

According to the rule of the“maximum membership”,the condition grades of the compressor,the turbine,the oil circuit and the gas seal were judged as all“excellent”while the process was as“poor”at 00:34 on April,7 of 2018.This indicated that the proposed fuzzycomprehensive assessment method successfully assessed the running condition of the large-scale centrifugal compressor set in a petrochemical corporation.Based on the above calculation results,the radar chart of the fuzzy comprehensive assessment of running condition is shown in Fig.6 at this specific moment.It can be seen intuitively from the chart that the process unit occurred an abnormal situation at the moment.

Table 9 Weight factor vectors of each index for all units

Fig.6.Radar chart of running condition for the centrifugal set at 00:34 on April 7 of 2018.

In practice,the index XI24666 indicating the throttle openness in the process unit of the set issued a threshold alarm at 00:44 on April 7 of 2018,which indicated that the process unit was suffering the assessment grade of“poor”.Fig.7 shows the time comparison of independent indexes between the proposed method and traditional way in the process unit.The proposed method achieved a warning of 10 min earlier than the traditional warning.

To further examine the constructed fuzzy comprehensive assessment method,the running condition at 03:08 on April 7 of 2018 was picked up as a comparative example as shown in Fig.8.It can be seen that the gas seal unit was graded as“good”,the process unit was as“medium”,and the others were as“excellent”,and the set was still in a normal running condition without warning.

4.Conclusions

The present work constructed a fuzzy comprehensive assessment method of running condition for a large-scale set.The MIC correlation analysis of indexes was introduced after the establishment of index system and condition assessment grade to eliminate highly correlated indexes.The second-order Markov chain model was used to calculate the predicted values of independent indexes with which the dynamic deterioration degrees were then determined.Afterwards,the membership matrix and the weight factor vector were calculated and the unit condition grade was attained using the rule of the“maximum membership”.The proposed method successfully assessed the running condition of a large-scale centrifugal compressor set in a petrochemical corporation.A radar chart was used to visualize the assessment results of all units in the set to overcome the deficiency of multilevel fuzzy comprehensive assessment.Comparison of warning time showed that the constructed method can achieve 10 min earlier alarm than the traditional way at the specific moment of 00:44 on April 7 of 2018.

Compared with the fault diagnosis,the real-time condition assessment and prediction are more helpful for practical operation.On-site operators are more concerned about the key issues such as whether or not the set can run,how long the set runs and how the variation trend of indexes changes.The running condition assessment and prediction technology is quite suitable for large-scale centrifugal compressor sets.The constructed method provides a useful tool for the continuous and safe operation of compressor sets and has a high engineering application value for avoiding unnecessary shutdown,shortening maintenance time and maximizing economic benefits towards corporations.