# Dynamic Early Warning Method for Major Hazard Installation Systems in Chemical Industrial Park

Authors: Yaguang Kong; Chenfeng Xie; Song Zheng; Peng Jiang; Meng Guan; Fang Wang
Journal: Complexity
 (2019)

---

## Abstract

The production and storage of major hazard installations (MHIs) bring potential risks to chemical industrial park (CIP). In the production system of MHIs, its dangerous degree is mainly determined by key parameters, and abnormal key parameters often lead to accidents. To predict the real-time risk values of MHIs and improve accident prevention ability of CIP, we need a method that can combine dynamic prediction and assessment. Quantitative risk assessment (QRA) is not capable of modelling risk variations during the operation of a process. Therefore, this paper adopts the data-driven approach. Inspired by visual qualitative analysis and quantitative analysis, a dynamic early warning method is proposed for MHIs. We can get the future trend of these key parameters by using strongly correlation variables to predict key parameters. Fuzzy evaluation analysis is performed on the risk levels of key parameters, and the dynamic evaluation index of these MHIs is obtained. This method can be applied to the dynamic evaluation of MHIs system in CIP. It can contribute to the safety of CIP in some aspects.


---

## Body

## 1. Introduction

The demand for everyday life materials has stimulated the vigorous development of the chemical industry. Large CIP have been built in various countries and regions. CIP have the characteristics of high complexity, strong association, and long production process. It is a process with multiple variables, strong coupling. Industrial process variables deviate from target values and cause a series of failure due to equipment failure, production equipment aging, personnel operation error, raw material characteristics, and changes in external environment [1]. It may also result in industrial accidents, such as fire, explosions, and leaks. In severe cases, this deviation may lead to casualties, major economic losses, and serious environmental pollution [2]. On December 3, 1984, the release of methyl isocyanate gas killed some 3000 people in Bhopal, India. This incident is considered to be the world's worst industrial chemical accident, on July 10, 1976, as a result of a reactor runaway accident in a chemical plant in Seveso, Italy. Dispersed dioxin into the local atmosphere led to skin disease chloracne in about 250 people and caustic soda burns in about 450 people [3]. On August 12, 2015, a hazardous chemical warehouse in Tianjin Binhai New Area exploded, leading to the death of 165 people, 8 missing, 798 injured, and 6.87-billion-yuan direct economic loss [4]. Therefore, accident prevention ability of CIP should be improved.QRA is a conventional risk modeling method. Although QRA methods are useful in modeling accident scenarios and quantifying risk, but their static structure is not capable of modeling risk variations during the process operation [5]. Dynamic early warning is a way to improve the preventive ability of the production system of MHIs in the chemical industry. In the area of early warning and assessment, with the use of the Process Resilience Analysis Framework, Jain et al. [6] proposed a resilience-based approach to managing uncertainties to better predict process upsets. Meel and Seider [7] proposed that the failure probabilities of safety systems and end-states are estimated using copulas and Bayesian analysis to ensure better predictions in plant-specific dynamic failure assessment. Cheng et al. [8] presented a method based on Bayesian and Vine Copula to predict the frequency of chemical accidents based on variations in different groups. Although Bayes’ theorem helps to obtain posterior probabilities, it necessitates identifying likelihood function, which is a difficult task if it is not a conjugate distribution to prior probability. Abdolhamidzadeh et al. [9] proposed the concept of ‘local domino effect’ to deal with a chain of accidents occurring within a process unit, achieving a good result to handle the uncertainty and the complexity. Pariyani et al. [10] introduced a distributed control system through bias analysis to calculate the probability of system failure, the risks associated with low probability, severe consequence incidents being predicted more accurately. Zhou et al. [11] established a security risk assessment method based on weighted fuzzy Petri net. This approach can easily model different relationships between the risk factors as well as their importance. In terms of social factor, Jain et al. [12] proposed a novel framework, Process Resilience Analysis Framework for incorporating both technical and social factors in an integrated approach resulting in a reduced frequency of loss of containment events. Aven and Ylönen [13] proposed a closer integration of the risk analysis and management approach and the sociotechnical perspective on safety can be used to improve risk and safety regulations. These methods take human factors, social factors, and domino effect into account, but they cannot dynamically assess the risk changes of hazards caused by abnormal parameter changes.The prediction of parameters is an important part of dynamic early warning. Other fields have been applied to investigate parameter prediction models [14–16]. However, these traditional methods face the background of big data with large computational complexity, multiple processes, and multivariable chemical processes [17]. For instance, due to the ease in configuring the alarms in control systems, the number of alarms in a plant has also gone up. This has led to frequent system problems, increase in the operator workload due to alarm overload, and industrial accidents. [18]. Hence, such methods are generally inapplicable. Today, big data have become an important strategic resource for a country [19]. Goel et al. [20] discuss the potential of big data analytics in the area of process safety and risk management in the energy industry. A large amount of data is also available in the CIP. In the early warning of CIP, conducting a reasonable analysis of big data and fully considering the correlation between key parameters of various process industries is essential for prediction analysis [21]. If all variables in the production system are not filtered, then the calculation is highly complex. Visual dynamic risk assessment method is an instrument that can qualitatively analyze the correlation between high-dimensional variables. Saatci [22] used Parallel Coordinates to explore the correlation between respiratory signals and airflow, relative temperature, relative humidity, and other factors and achieved satisfactory results. Azhar and Rissanen [23] evaluated an application of Parallel Coordinates for interactive filtering of alarm data by comparing its user performance against typical alarm lists. The results showed that Parallel Coordinates reduce alarm filtering time and reduce human mistakes. Visualization method can only simultaneously perform qualitative analysis, and the correlation between variables should be quantitatively analyzed. Berthold and Höppner [24] investigated clustering time series based on Euclidean distance and Pearson correlation and contributed to future research.The prediction of MHIs should be performed rapidly. Excessive delay will result in unexpected consequences. Zhao et al. [25] used GPU to conduct an accelerated operation and effectively improve the prediction speed. Ruiz et al. [26] employed a genetic algorithm to optimize neural network weights, enhance prediction accuracy, and reduce prediction time. Chemical process industry is associated with many variables. To handle the high-dimensional variables, Erkmen and Yıldırım [27] used principal components analysis as a feature extraction method to promote classification performance. Support vector machine [28], back propagation neural network (BP neural network) [29], and probabilistic models are commonly used in traditional prediction models [30]. However, these methods can be implemented to static prediction only. When these methods are applied to the prediction of industrial process dynamics, a dynamic time modeling problem is transformed into a stationary space modeling problem, which inevitably leads to a decline in prediction accuracy. We need to construct the model; it has a feedback link, which can later memorize the state information of the key parameters of the chemical industry, allowing the system to adapt to time-varying characteristics. For example, Qin et al. [31] used an Elman neural network to improve the smooth transition period of an autoregressive model and accurately predicted wind energy and speed.The present paper is divided as follows. The overall process of Section2 briefly introduced a structure of dynamic early warning method. Section 3 presents the dynamic prediction and evaluation method, which combines variable correlation analysis, and feature analysis and feature vector acquisition, to conduct strong correlation analysis and dimensionality reduction. Construct a recurrent neural network to predict MHIs. Finally, we use the MHIs risk level to defining a method for the dynamic evaluation. Section 4 takes Quzhou City CIP as an example to perform a case study of the method. Section 5 provides some concluding remarks.


## 2. Dynamic Prediction and Evaluation Method for Major Hazard Installations

Given that the chemical accidents are usually caused by some parameter fluctuations; therefore, it is necessary to construct a dynamic prediction model. The applications of neural network in various fields have continuously emerged recently. Compared with support vector machine [28] and nonlinear regressions, neural network possesses the robust multivariable prediction performance and can constantly adjust the weight among input layer, output layer, and hidden layer to achieve high-precision prediction [29]. In traditional neural network, we add a state layer to construct the dynamic recurrent network, which can better adapt to the dynamic change characteristics of chemical production parameters. Analyzing the strong correlation variables of key parameters is crucial due to numerous parameters and high complexity of the chemical process. Correlation analysis work, such as that of Cheng et al. [8], uses the Copula function to establish the asymmetric correlation model of anomalous number of events, considering the correlation among the four teams. Under the multivariable condition of chemical production engineering, the process of establishing Copula function is considerably complicated, and a fast and accurate method of correlation analysis is necessary. The visualization method can be used for fast qualitative analysis of relevant variables and, combined with quantitative analysis method, effectively implement the selection of strong correlation variables of key parameter.The structural block diagram of the dynamic early warning method is presented in Figure1. The first part of the algorithm is variable correlation analysis, the second part is feature analysis and feature vector acquisition, and the third part is dynamic recurrent prediction. And the last part is the dynamic risk assessment. The implementation process of the algorithm is as follows.Figure 1
Structure block diagram of dynamic early warning method.(1) All correlation variables of the key parameters of MHIs in the CIP use the qualitative analysis of variable correlation by data visualization [32].(2) The average value and standard deviation of key parameters and related variables are calculated.(3) The correlation coefficient r value between two variables is obtained and the strong correlation variable. (4) The covariance matrix of the strong correlation variable is also calculated, and its eigenvalue and the eigenvector are obtained.(5) The order of the eigenvalue from large to small is selected, and the maximum is chosen according to the established rule. K eigenvalue constitutes the corresponding eigenvector matrix.(6) The samples of all relevant variables are projected into the selected K feature vectors.(7) The strong correlation matrix is used for dynamically recursive prediction.(8) The strong correlation variables used dynamic recurrent prediction and evaluated the risk level using the dynamics prediction results of each parameter. The assessment standard is divided into two levels ultimately determining the risk level of CIP.


## 3. The Proposed Methodology

### 3.1. Correlation Analysis of Variables

Firstly, there is a qualitative analysis of the correlative variables of key parameters of MHIs, and these correlation variables are composed into a high-dimensional data set. Each variable is divided into a coordinate axis that is equal in distance and parallel to each other. Each axis represents an attribute dimension, and the value of the variable corresponds to the corresponding position on the axis. In this way, every chemical process variable can be divided into a broken line on the parallel axis of n according to its attribute value. The implementation steps are described in the following.The associated variables of production key parameters in the chemical industry MHIs system are combined into the following scoring matrix S:(1)S=S11S12⋯S1pS21S22⋯S2p⋮⋮⋱⋮Sn1Sn2⋯Snp.In the plane with the Descartes coordinate system, each row in the score matrix S,Sn=Sn1Sn2⋯Snp corresponds to a fold line connected by a k-1-line segment in a parallel coordinate system. These data must be converted from the Descartes coordinate system to the parallel coordinate system because each coordinate axis is equal in length; thus, the scoring matrix is homogenized 0–1. The relationship between the processed data sn,p′ and the original data sn,p is sn,p′=Sn,p-Sn,min/sn,max-⁡Snmin, where Sn,min and Sn,max are the maximum and minimum values of the correlation variable, respectively.Parallel axes represent the correlation variables of each key parameter in MHIs, which are connected by the parameter values of these correlation variables. If the line between the two variables is crossed to show the X type, then the two variables are negatively correlated. If the line between the two variables is parallel, then the two variables indicate positive correlation. If the lines between variables are randomly crossed, then no unique relationship exists between the two variables.Then, quantitative analyzed the correlation variables of key parameters of MHIs [33]. Strong and medium correlation variables are further divided to verify the visualization method. The correlation coefficient is reflected in the symbol “r,” and the Pearson correlation coefficient is as follows:(2)rA,B=1n-1∑i=1nAi-AvgAδABi-AvgBδB,where n is equal to the number of selected data items.Avg(A) and Avg(B) are the average values of selected data on attributes A and B, respectively, and δA,δB are their corresponding standard deviations.


### 3.2. Feature Analysis and Feature Vector Acquisition

Concerning the majority of chemical processes, the key parameters have numerous strong correlation variables. These strongly correlation variables constitute a high-dimensional data, which must be reduced in order to improve the efficiency and accuracy of prediction of key parameters [34], n-dimensional features to K dimensions (k<n). The steps are as follows.(i) The strong correlation variables are selected from variable correlation analysis constituting the data matrix X:(3)X=x11x12⋯x1mx21x22⋯x2m⋮⋮⋱⋮xn1xn2⋯xnm.(ii) The dimensionality of each data matrix dimension has a serious influence on the principal component, which may cause a considerable difference between the number and the value of the principal component and the actual situation. The strong correlation variables constitute the data matrix X, and each variable is standardized reprocessed to obtained the standardized data matrix Z:zn,m=xnm-xn-/Ym2. The overall mean of xm is x-m=1/n∑n=1nxn,m, and the total variance of xm is Ym2=1/n-1∑n=1nxn,m-x-n2. The correlation matrix, C=1/nZTZ, is calculated, in which ZT is the transpose of the standardized data matrix Z, n represents the row of matrix Z, m represents the column of matrix Z, and the calculated correlation matrix C comprises the D principal component.(iii) Jacobi’s method is employed to find the eigenvectorwi,i=1,2,…,m of the correlation matrix C and eigenvalue matrix Λ=diag⁡λ1,λ2,…,λm, where λ is the eigenvalue and diag represents a diagonal matrix.(iv) The eigenvalues are arranged in order from large to smallλ1>λ2>⋯>λm. The order of the eigenvector columns is adjusted correspondingly, allowing the maximum variance to be the first principal components, the sublarge variance as the second principal components, and the least variance as the D principal component.(v) The K principal component of the maximum variance is selected to allow the K principal component to contain most of the original data information. The cumulative variance contribution of the commonly selected K principal components is 85% larger than that of the total variance. That is,∑j=1kλj/∑j=1dλj≥85%.(vi) By selecting the eigenvectorswi,i=1,2,…,k, of K principal components, k independent linear combination of new variables is obtained:(4)ξ1=x1w11+x2w21+⋯+xdwd1ξ2=x1w12+x2w22+⋯+xdwd2⋮ξk=x1w1k+x2w2k+⋯+xdwdk.(vii)ξ1,ξ2,…,ξk is the rebuilt K feature after the reduction of dimension. These features represent the principal components of all strong correlation variables, and the matrix T of the K principal components of the strong correlation variable is as follows:(5)T=x1w11x2w21⋯xdwd1x1w12x2w22⋯xdwd2⋮⋮⋱⋮x1w1kx2w2k⋯xdwdk.


### 3.3. Dynamic Recurrent Prediction

In the prediction model, in traditional neural network, we add a state layer to construct the dynamic recurrent network. As in Figure2, which shows four layers, namely, input layer, hidden layer, state layer, and output layer [35].Figure 2
Recurrent neural network structures.The mathematical model is as follows:(6)xt=fw1xct+w2ut-1+θ1,xct=xt-1,yt=w3xt+θ2,where y(t) is the output of the output layer, u(t) is the external input of the input layer, the matrix T of the K principal component of the strong correlation variable is input,xc(t) is the output of the state layer, x(t) is the output of the hidden layer, w(1) is the connection weight of the state layer and the hidden layer, and w(2) is the connection weight between the input layer and the hidden layer. w(3) is the connection weight of the hidden layer and the output layer, the θ(1) is the hidden layer threshold, and the θ(2) is the output threshold value, where f represents the activation function of the neural network using the Sigmoid function presented in formula (7):(7)fx=11+exp⁡-x,0<fx<1.Formula (6) can be introduced as follows:(8)xct=fwk-11xct-1+wk-22xct-2,wherewk-1(1) and wk-2(2) represent the connection weight at previous different times. Formula (8) indicates that xc(t) is related to the weight of the connection at the previous time and realizes the characteristic of dynamic recursion [36].The input and output of input layer are respectively(9)Iα1t=uαt-1,(10)Oα1t=Iα1t.α = 1 ,  2 ,  ⋯  E 1, where the input and output of the hidden layer are respectively(11)Iβ2t=∑α=1E1wβα2Oα1t+∑β=1E2wβγ1Oγct+θ1,(12)Oβ2t=fIβ2t=11+ⅇIβ2t.β = 1 ,  2 ,  ⋯  E 2. The input and output of the state layer are respectively as follows:(13)IγCt=Oβ2t-1,(14)OγCt=IγCt=xCt.γ = 1 ,  2 ,  ⋯  E 2. The input and output of the output layer are respectively(15)Iμ3t=∑j=1s2wμβ3Oβ2t+θ2,(16)yμt=Oμ3t=gI˙μ3t.whereμ=1,2,⋯E3. Among them, E1, E2, and E3 are the input, hidden, and output layers, respectively, and the number of layers of the state layer is the same as that of the hidden layer. The final yt is the predictive value of the key parameters.


### 3.4. Dynamic Risk Assessment and Early Warning

Fuzzy comprehensive evaluation is a comprehensive decision-making methodology of a multivariable problem-solving complex decision process [37]. It is developed from fuzzy sets. We through the first-level and second-level fuzzy comprehensive evaluation to determine the risk R value of CIP. Hierarchical chart is presented in Figure 3. The specific steps are as follows.Figure 3
Hierarchical chart of dynamic assessment of MHIs.Step 1 (construction of second-level fuzzy comprehensive evaluation).
The second-level fuzzy comprehensive evaluation defines an evaluation matrix of key parametersQpσp=[Q1σ1,Q2σ2,…,Qpσp], p is the number of MHIs, and σ1,σ2,…,σp is the number of key parameters for each MHIs; they may be different values. The percentage of the predicted value of the key parameter y(t) exceeding the high alarm value G: Cpσp=yt-G/G×100% multiplies the corresponding elements between Qpσp and Cpσp. Then, the second-level fuzzy comprehensive evaluation index matrix is obtained. (17)Apσp=Cpσp·Qpσp.Step 2 (construction of the first-level fuzzy comprehensive evaluation).
According to the risk degree of each key parameter in MHIs, the weight of the first-level fuzzy comprehensive evaluation is determined:Dp=[D1,D2,…,Dp]. The second-level fuzzy comprehensive evaluation matrix Apσp=[A1σ1,A2σ2,…,Apσp] obtained by (6) and the weight of the first-level fuzzy comprehensive evaluation are conducted by point multiplication, and the risk grade index of the MHIs in the CIP is obtained as follows:(18)Bp=Apσp·Dp,
wherep is the total number of MHIs in CIP.Finally, the dynamic level of MHIs in CIP is determined by the evaluation result. The score of dynamic evaluation is decided by the scoring method, and the range of the grade is provided in Table1. R, which is specified in Identification and Control of MHIs (GB18218), is used as a grading index [38]. Here, R obtains the value of hazard installations dynamic index Bp. The historical factors of MHIs may affect the overall risk. The data of all dynamic risk grades are derived from the last time the true value of the historical key parameter exceeds the high alarm value to return to the safety value or the current value.Table 1
Dynamic risk level.
comment grade range Disastrous I R ≥ 100 High risk II 50 ≤  R <  100 Moderate risk III 10 ≤ R < 50 Low risk IV R < 10Step 3.
The final dynamic assessment of MHIs and early warning results are used to urge security personnel to exclude the potential risks of MHIs for the first time.


## 3.1. Correlation Analysis of Variables

Firstly, there is a qualitative analysis of the correlative variables of key parameters of MHIs, and these correlation variables are composed into a high-dimensional data set. Each variable is divided into a coordinate axis that is equal in distance and parallel to each other. Each axis represents an attribute dimension, and the value of the variable corresponds to the corresponding position on the axis. In this way, every chemical process variable can be divided into a broken line on the parallel axis of n according to its attribute value. The implementation steps are described in the following.The associated variables of production key parameters in the chemical industry MHIs system are combined into the following scoring matrix S:(1)S=S11S12⋯S1pS21S22⋯S2p⋮⋮⋱⋮Sn1Sn2⋯Snp.In the plane with the Descartes coordinate system, each row in the score matrix S,Sn=Sn1Sn2⋯Snp corresponds to a fold line connected by a k-1-line segment in a parallel coordinate system. These data must be converted from the Descartes coordinate system to the parallel coordinate system because each coordinate axis is equal in length; thus, the scoring matrix is homogenized 0–1. The relationship between the processed data sn,p′ and the original data sn,p is sn,p′=Sn,p-Sn,min/sn,max-⁡Snmin, where Sn,min and Sn,max are the maximum and minimum values of the correlation variable, respectively.Parallel axes represent the correlation variables of each key parameter in MHIs, which are connected by the parameter values of these correlation variables. If the line between the two variables is crossed to show the X type, then the two variables are negatively correlated. If the line between the two variables is parallel, then the two variables indicate positive correlation. If the lines between variables are randomly crossed, then no unique relationship exists between the two variables.Then, quantitative analyzed the correlation variables of key parameters of MHIs [33]. Strong and medium correlation variables are further divided to verify the visualization method. The correlation coefficient is reflected in the symbol “r,” and the Pearson correlation coefficient is as follows:(2)rA,B=1n-1∑i=1nAi-AvgAδABi-AvgBδB,where n is equal to the number of selected data items.Avg(A) and Avg(B) are the average values of selected data on attributes A and B, respectively, and δA,δB are their corresponding standard deviations.


## 3.2. Feature Analysis and Feature Vector Acquisition

Concerning the majority of chemical processes, the key parameters have numerous strong correlation variables. These strongly correlation variables constitute a high-dimensional data, which must be reduced in order to improve the efficiency and accuracy of prediction of key parameters [34], n-dimensional features to K dimensions (k<n). The steps are as follows.(i) The strong correlation variables are selected from variable correlation analysis constituting the data matrix X:(3)X=x11x12⋯x1mx21x22⋯x2m⋮⋮⋱⋮xn1xn2⋯xnm.(ii) The dimensionality of each data matrix dimension has a serious influence on the principal component, which may cause a considerable difference between the number and the value of the principal component and the actual situation. The strong correlation variables constitute the data matrix X, and each variable is standardized reprocessed to obtained the standardized data matrix Z:zn,m=xnm-xn-/Ym2. The overall mean of xm is x-m=1/n∑n=1nxn,m, and the total variance of xm is Ym2=1/n-1∑n=1nxn,m-x-n2. The correlation matrix, C=1/nZTZ, is calculated, in which ZT is the transpose of the standardized data matrix Z, n represents the row of matrix Z, m represents the column of matrix Z, and the calculated correlation matrix C comprises the D principal component.(iii) Jacobi’s method is employed to find the eigenvectorwi,i=1,2,…,m of the correlation matrix C and eigenvalue matrix Λ=diag⁡λ1,λ2,…,λm, where λ is the eigenvalue and diag represents a diagonal matrix.(iv) The eigenvalues are arranged in order from large to smallλ1>λ2>⋯>λm. The order of the eigenvector columns is adjusted correspondingly, allowing the maximum variance to be the first principal components, the sublarge variance as the second principal components, and the least variance as the D principal component.(v) The K principal component of the maximum variance is selected to allow the K principal component to contain most of the original data information. The cumulative variance contribution of the commonly selected K principal components is 85% larger than that of the total variance. That is,∑j=1kλj/∑j=1dλj≥85%.(vi) By selecting the eigenvectorswi,i=1,2,…,k, of K principal components, k independent linear combination of new variables is obtained:(4)ξ1=x1w11+x2w21+⋯+xdwd1ξ2=x1w12+x2w22+⋯+xdwd2⋮ξk=x1w1k+x2w2k+⋯+xdwdk.(vii)ξ1,ξ2,…,ξk is the rebuilt K feature after the reduction of dimension. These features represent the principal components of all strong correlation variables, and the matrix T of the K principal components of the strong correlation variable is as follows:(5)T=x1w11x2w21⋯xdwd1x1w12x2w22⋯xdwd2⋮⋮⋱⋮x1w1kx2w2k⋯xdwdk.


## 3.3. Dynamic Recurrent Prediction

In the prediction model, in traditional neural network, we add a state layer to construct the dynamic recurrent network. As in Figure2, which shows four layers, namely, input layer, hidden layer, state layer, and output layer [35].Figure 2
Recurrent neural network structures.The mathematical model is as follows:(6)xt=fw1xct+w2ut-1+θ1,xct=xt-1,yt=w3xt+θ2,where y(t) is the output of the output layer, u(t) is the external input of the input layer, the matrix T of the K principal component of the strong correlation variable is input,xc(t) is the output of the state layer, x(t) is the output of the hidden layer, w(1) is the connection weight of the state layer and the hidden layer, and w(2) is the connection weight between the input layer and the hidden layer. w(3) is the connection weight of the hidden layer and the output layer, the θ(1) is the hidden layer threshold, and the θ(2) is the output threshold value, where f represents the activation function of the neural network using the Sigmoid function presented in formula (7):(7)fx=11+exp⁡-x,0<fx<1.Formula (6) can be introduced as follows:(8)xct=fwk-11xct-1+wk-22xct-2,wherewk-1(1) and wk-2(2) represent the connection weight at previous different times. Formula (8) indicates that xc(t) is related to the weight of the connection at the previous time and realizes the characteristic of dynamic recursion [36].The input and output of input layer are respectively(9)Iα1t=uαt-1,(10)Oα1t=Iα1t.α = 1 ,  2 ,  ⋯  E 1, where the input and output of the hidden layer are respectively(11)Iβ2t=∑α=1E1wβα2Oα1t+∑β=1E2wβγ1Oγct+θ1,(12)Oβ2t=fIβ2t=11+ⅇIβ2t.β = 1 ,  2 ,  ⋯  E 2. The input and output of the state layer are respectively as follows:(13)IγCt=Oβ2t-1,(14)OγCt=IγCt=xCt.γ = 1 ,  2 ,  ⋯  E 2. The input and output of the output layer are respectively(15)Iμ3t=∑j=1s2wμβ3Oβ2t+θ2,(16)yμt=Oμ3t=gI˙μ3t.whereμ=1,2,⋯E3. Among them, E1, E2, and E3 are the input, hidden, and output layers, respectively, and the number of layers of the state layer is the same as that of the hidden layer. The final yt is the predictive value of the key parameters.


## 3.4. Dynamic Risk Assessment and Early Warning

Fuzzy comprehensive evaluation is a comprehensive decision-making methodology of a multivariable problem-solving complex decision process [37]. It is developed from fuzzy sets. We through the first-level and second-level fuzzy comprehensive evaluation to determine the risk R value of CIP. Hierarchical chart is presented in Figure 3. The specific steps are as follows.Figure 3
Hierarchical chart of dynamic assessment of MHIs.Step 1 (construction of second-level fuzzy comprehensive evaluation).
The second-level fuzzy comprehensive evaluation defines an evaluation matrix of key parametersQpσp=[Q1σ1,Q2σ2,…,Qpσp], p is the number of MHIs, and σ1,σ2,…,σp is the number of key parameters for each MHIs; they may be different values. The percentage of the predicted value of the key parameter y(t) exceeding the high alarm value G: Cpσp=yt-G/G×100% multiplies the corresponding elements between Qpσp and Cpσp. Then, the second-level fuzzy comprehensive evaluation index matrix is obtained. (17)Apσp=Cpσp·Qpσp.Step 2 (construction of the first-level fuzzy comprehensive evaluation).
According to the risk degree of each key parameter in MHIs, the weight of the first-level fuzzy comprehensive evaluation is determined:Dp=[D1,D2,…,Dp]. The second-level fuzzy comprehensive evaluation matrix Apσp=[A1σ1,A2σ2,…,Apσp] obtained by (6) and the weight of the first-level fuzzy comprehensive evaluation are conducted by point multiplication, and the risk grade index of the MHIs in the CIP is obtained as follows:(18)Bp=Apσp·Dp,
wherep is the total number of MHIs in CIP.Finally, the dynamic level of MHIs in CIP is determined by the evaluation result. The score of dynamic evaluation is decided by the scoring method, and the range of the grade is provided in Table1. R, which is specified in Identification and Control of MHIs (GB18218), is used as a grading index [38]. Here, R obtains the value of hazard installations dynamic index Bp. The historical factors of MHIs may affect the overall risk. The data of all dynamic risk grades are derived from the last time the true value of the historical key parameter exceeds the high alarm value to return to the safety value or the current value.Table 1
Dynamic risk level.
comment grade range Disastrous I R ≥ 100 High risk II 50 ≤  R <  100 Moderate risk III 10 ≤ R < 50 Low risk IV R < 10Step 3.
The final dynamic assessment of MHIs and early warning results are used to urge security personnel to exclude the potential risks of MHIs for the first time.


## 4. Practical Application

CIP is composed of many MHIs; the air separation distillation section is one of the typical representatives of MHIs. This practical application relies on the industrial data of a large CIP in Quzhou City from July 2017, and a series of model establishment, correlation analysis, prediction, and assessment is completed to realize the dynamic early warning function of MHIs in the CIP.This practical application selects an industry of air separation distillation section in the CIP as a dynamic early warning object of MHIs [39], a simplified graphic sketch of the air separation distillation section in Figure 4. This section simplifies and removes some industrial process from the air separation distillation section. To make it easier to understand the industrial process, we are based on data-driven method and can neglect to describe the mechanism of the chemical process industry. The section can separate air into oxygen and nitrogen. Risk degree of the section is affected by several key links, such as air separation tower system, air purification system, and air precooling system. These production lines are strictly dependent on pressure and temperature, if pressure and temperature are abnormal, they may result in an explosion and fire. Air separation tower system consists of Turboexpander and distillation column, etc. We removed several weakly parameters and links for ensuring factories data privacy, but did not affect our application result. We take liquid air level in lower column LI101 as the key parameter in air separation distillation section. Since most of the top feed comes from the bottom of the tower, the working conditions at the bottom of the tower play a decisive role in the air separation distillation section. The abnormal liquid air level in lower column will result in the decrease of oxygen purity, which poses a serious threat to the normal production of the system. The section has 38 liquid levels, pressure, temperature, and other sensor variables, as shown in Table 2, in which the liquid air level in lower column of the key parameter is LI101. In order to predict LI101, we need to find the strong correlation variables in 38 variables. The change trend of LI101 is predicted by the subtle changes in these strong correlation variables. It should be noted that LI101 is only one of many key parameters, analysis process of other key parameters is similar enough to analysis process of LI101 so that they are not repeated here.Table 2 item Description PI101 Air intake chamber pressure (MPa) PI102 Nitrogen Outlet Tower Pressure (MPa) PI103 Oxygen-enriched air outlet pressure (MPa) PI1201 Air outlet left adsorption bucket pressure (MPa) PI1202 Air outlet right adsorption barrel pressure(MPa) PI1203 Instrument air pressure (MPa) PI01 Pressure of air outlet main heat exchanger(MPa) PI02 Lower column pressure (MPa) PI03 Evaporator condenser pressure (MPa) PDI01 The resistance of lower column (MPa) PI401 Turboexpander1# Intake pressure (MPa) PI402 Turboexpander1#Exhaust pressure (MPa) PI403 Turboexpander2#Intake pressure (MPa) PI404 Turboexpander2#Exhaust pressure (MPa) PI2001 Turboexpander1#Bearing air pressure (MPa) PI2002 Turboexpander2# Bearing air pressure (MPa) PI2003 Turboexpander1#Sealing pressure (MPa) PI2004 Turboexpander2#Sealing pressure (MPa) TI101 Air temperature before entering tower (°C) TI102 Outlet Tower Temperature of Contaminated Nitrogen (°C) TI103 Outlet Tower Temperature of Contaminated Nitrogen (°C) TI1201 Air outlet left adsorption bucket temperature (°C) TI1202 Air outlet right adsorption bucket temperature (°C) TE1204 Temperature of regenerated gas outlet adsorption bucket (°C) TI01 Temperature of main air outlet heat exchanger (°C) TI02 Turboexpander1# Intake temperature (°C) TI03 Turboexpander2# Intake temperature (°C) TI04 Turboexpander1# Exhaust temperature (°C) TI05 Turboexpander2# Exhaust temperature (°C) TI07 Exit temperature of precooler (°C) FI1201 Regenerative gas flow rate (Nm3/h) FI101 Air flow rate (Nm3/h) FI102 Nitrogen flow rate (Nm3/h) LI101 liquid air level in lower column (mm) LI01 Water cooled liquid level (mm) LI02 Liquid oxygen level in the main condenser (mm) SI402 Expander speed (RPM) LV02 Liquid-air throttle valve (%)Figure 4
Air separation distillation section.The 38 variables of air separation distillation section are analyzed by visual correlation analysis. The experimental datasets comprise 1000 sets, as shown in Figure5. The value of LI101 exceeds the high alarm value when the broken line between two adjacent variables is red. Conversely, the blue line indicates that the LI101 is not exceeding the high alarm value.Figure 5
Visual correlation analysis of 38 variables.The 38 variables are not convenient to observe. So, analyze the correlation between 22 variables on the local graph. As shown in Figure6, the lines between LI101 and adjacent variables PI404 and PI403 are randomly crossed, indicating that no special correlation exists between these variables.Figure 6
Visual correlation analysis local graph A.As shown in Figure7, lines between LI101 and PI2002, LI101 and PI2003 are overall parallel, indicating a positively correlated relationship. This method can intuitively and rapidly select strong correlation of key parameter.Figure 7
Visual correlation analysis local graph B.After visual correlation analysis, we filtered 27 correlation variables and eliminated 11 weakly correlation variables.Next, in order to verify the accuracy of visual correlation analysis and further filter strong correlation variables, Pearson correlation analysis of 27 correlation variables and the results of the filtered strong correlation variables are shown in Table3.Table 3
Strong correlation variable table.
item PI02 PI03 PI401 PI402 PI1201 PI2003 PI102 correlation 0.674∗∗ 0.661∗∗ 0.645∗∗ 0.621∗∗ 0.670∗∗ 0.669∗∗ 0.674∗∗ conspicuousness 0.000 0.000 0.000 0.000 0.000 0.000 0.000 N 1000 1000 1000 1000 1000 1000 1000  item PI103 FI102 FI101 FI1201 TI03 TI04 TI102  correlation 0.659∗∗ 0.659∗∗ 0.602∗∗ 0.697∗∗ 0.617∗∗ -0.611∗∗ -0.633∗∗ conspicuousness 0.000 0.000 0.000 0.000 0.000 0.000 0.000 N 1000 1000 1000 1000 1000 1000 1000The table of strong correlation variable has three rows. The first row is the correlation value, in which 0 to 0.33 is weak correlation, 0.33 to 0.67 is medium correlation, and 0.67 to 1 is strong correlation [40]. The second row is the significant test result, that is, sig bilateral test. The asterisk behind the correlation row values represents the degree of saliency: ∗ representing 0.01 < sig < 0.05 indicates significance; ∗∗ representing sig < 0.01 denotes a high degree of significance; and the third row is the sample capacity.Take the data of PI02 in Table3 as an example, where correlation coefficient r = 0.674, significant P = 0, and sample capacity N = 1000. The results indicate that PI02 is the lower column pressure, and this pressure has a strong correlation with LI101.The strong correlation variable in Table3 is used to reduce the dimension.The data matrix X is composed of strong correlation variables, and the matrix is normalized. The normalized data matrix Z is obtained, and the correlation matrix C is calculated. The correlation matrix C comprises 14 principal components. The correlation coefficient matrix of solid correlation variables is shown in Table4.Table 4
Correlation coefficient matrix.
Z-Score PI02 PI03 PI401 PI402 PI1201 PI2003 PI102 PI103 FI102 FI101 FI1201 TI03 TI04 TI102 PI02 1.000 0.354 0.744 0.376 0.720 0.805 0.519 0.884 -0.468 0.396 0.250 0.798 0.114 0.749 PI03 0.354 1.000 0.854 0.976 0.694 0.714 0.699 0.477 -0.623 0.984 0.866 0.156 0.543 0.850 PI401 0.744 0.854 1.000 0.856 0.878 0.933 0.727 0.779 -0.675 0.871 0.750 0.505 0.463 1.000 PI402 0.376 0.976 0.856 1.000 0.637 0.695 0.659 0.462 -0.512 0.993 0.862 0.218 0.444 0.855 PI1201 0.720 0.694 0.878 0.637 1.000 0.958 0.798 0.828 -0.902 0.693 0.650 0.379 0.677 0.868 PI2003 0.805 0.714 0.933 0.695 0.958 1.000 0.785 0.873 -0.766 0.734 0.599 0.525 0.511 0.930 PI102 0.519 0.699 0.727 0.659 0.798 0.785 1.000 0.674 -0.742 0.696 0.670 0.233 0.531 0.716 PI103 0.884 0.477 0.779 0.462 0.828 0.873 0.674 1.000 -0.695 0.509 0.331 0.726 0.288 0.780 FI102 -0.468 -0.623 -0.675 -0.512 -0.902 -0.766 -0.742 -0.695 1.000 -0.590 -0.611 -0.144 -0.796 -0.657 FI101 0.396 0.984 0.871 0.993 0.693 0.734 0.696 0.509 -0.590 1.000 0.866 0.211 0.515 0.869 FI1201 0.250 0.866 0.750 0.862 0.650 0.599 0.670 0.331 -0.611 0.866 1.000 -0.038 0.672 0.738 TI03 0.798 0.156 0.505 0.218 0.379 0.525 0.233 0.726 -0.144 0.211 -0.038 1.000 -0.313 0.519 TI04 0.114 0.543 0.463 0.444 0.677 0.511 0.531 0.288 -0.796 0.515 0.672 -0.313 1.000 0.441 TI102 0.749 0.850 1.000 .855 0.868 0.930 0.716 0.780 -0.657 0.869 0.738 0.519 0.441 1.000The correlation coefficient matrix shows a strong correlation between the indexes. For example, the correlation coefficient between PI03 and PI402 and FI101 is large, as in Table4. This result shows that an overlap exists between their index information. Next, the Jacobi method is used to find the eigenvector wi,i=1,2,…,m of the correlation matrix Λ=diag⁡λ1,λ2,…,λm, where λ is the eigenvalue and diag is a diagonal matrix. The eigenvalues are sorted by the order of λ1>λ2>⋯>λm, and the order of the eigenvector column is adjusted correspondingly, making the maximum variance the first principal components, the second largest variance as the second principal components, and the minimum variance as the thirteenth principal components.As shown in Table5, the principal component characteristic root and contribution rate are characteristic root λ1=9.543, characteristic root λ2=2.328, and characteristic root λ3=1.23. The cumulative variance contribution rate of the first three principal components reached 93.597%, which covers most of the information. This finding suggests that the first three principal components can represent the 13 principal components; therefore, the first three indexes can be extracted. The principal component is taken as ξ1,ξ2, and ξ3. The principal component load vector is calculated by the three principal components, and the result is presented in Table 6.Table 5
Principal component characteristic root and contribution rate.
assembly Initial Eigenvalues Extraction of load square sum total variance cumulative% total variance cumulative% percentage percentage 1 9.543 68.165 68.165 9.543 68.165 68.165 2 2.328 16.631 84.796 2.328 16.631 84.796 3 1.232 8.800 93.597 1.232 8.800 93.597 4 0.341 2.434 96.031 5 0.204 1.454 97.484 6 0.139 0.990 98.474 7 0.076 0.543 99.017 8 0.061 0.432 99.450 9 0.045 0.322 99.771 10 0.021 0.149 99.920 11 0.007 0.051 99.971 12 0.002 0.017 99.988 13 0.002 0.012 100.000 14 5.197E-5 0.000 100.000Table 6
Principal component load vector.
assembly 1 2 3 PI02 0.7 0.654 −0.06 PI03 0.874 -0.338 0.302 PI401 0.97 0.093 0.154 PI402 0.851 -0.284 0.429 PI1201 0.936 0.074 -0.314 PI2003 0.948 0.212 -0.118 PI102 0.833 -0.075 -0.165 PI103 0.801 0.515 -0.199 PI102 -0.805 0.135 0.512 FI101 0.884 -0.289 0.34 FI1201 0.79 -0.5 0.155 TI03 0.419 0.841 0.224 TI04 0.596 -0.552 -0.617 TI102 0.965 0.109 0.174The indexes PI02, PI03, PI401, PI402, PI1201, PI2003, PI102, PI103, PI102, FI101, FI1201, and TI102 have high load on the first principal component, and the correlation is high. TI03 has high load on the second principal components and a strong correlation. TI04 has high load on the third principal components and a strong correlation.The load matrix of the principal component is the eigenvector of the principal component. Principal component load vector is divided by the arithmetic square root of the principal component characteristic root, that is, the load matrix of the principal component:(19)Ui=Aiλi.The results are shown in Table7. The expression of the principal component coefficient is(20)ξ1=x1U1+x2U2+⋯+xdUd.Table 7
Load matrix of the principal component.
assembly 1 2 3 PI02 0.23 0.43 -0.05 PI03 0.28 -0.22 0.27 PI401 0.31 0.06 0.14 PI402 0.28 -0.19 0.39 PI1201 0.3 0.05 -0.28 PI2003 0.31 0.14 -0.11 PI102 0.27 -0.05 -0.15 PI103 0.26 0.34 -0.18 PI102 -0.26 0.09 0.46 FI101 0.29 -0.19 0.31 FI1201 0.26 -0.33 0.14 TI03 0.14 0.55 0.2 TI04 0.19 -0.36 -0.56 TI102 0.31 0.07 0.16The eigenvector is the load matrixUi,i=1,2,…,d, and the eigenvalue is the value of other variables at the same time.Table7 uses the three principal component to dynamic recurrent prediction.


### 4.1. Prediction Model

BP neural network is a feedback-free feedforward neural network, which is composed of a set of interconnected neurons, in which each neuron is connected with the corresponding weight [41]. The extreme learning machine (ELM) is a single hidden layer feedforward neural network. It is made up by a simple structure with strong learning and generalization ability. This model has been extensively used in the fields of pattern recognition and regression estimation [42].Through filtered of key parameter correlation variables and dimensionality reduction, 14 groups of strong correlation variables were filtered from 38 groups of correlation variables. After dimensionality reduction, 3 groups of principal components representing 14 groups of strong correlation variables were obtained.The following is the prediction of key parameter. The method is compared with BP neural network, ELM, and dynamic recurrent prediction model without dimensionality reduction. For fairness, the selected variables are also analyzed by visual qualitative analysis and quantitative correlation analysis. The predicted experimental data are derived from 14 sets of solid correlation variables, which remove abnormal values data and acquire 1000 sets of effective data as experimental samples. The distribution of the sample set is presented in Table8.Table 8
Sample set allocation table.
Data sets prorate% amount Training 70% 700 Validation 15% 150 Testing 15% 150In this study, the configuration of hardware involves CPU of Intel Core i5-7300HQ, with main clock frequency of 2.50 GHz and GPU is GeForce GTX1050Ti by NVIDIA. This study adopts the following international evaluation index to measure the prediction efficiency of the model.(1) Maximum Error (ME)(21)ME=max⁡fi-yi,i=1,2,3,…,n.Among them, n is the sample size,fi is the predicted value, and yi is the true value.(2) Mean Absolute Error (MAE) is the average of the absolute value of the deviation between all observed and predicted values.(22)MAE=1n∑i=1i=nfi-yi.(3) Root Mean Square Error (RMSE) represents the discreteness of the error distribution. A small RMSE leads to concentrated error distribution and high prediction accuracy.(23)RMSE=1n∑i=1i=nfi-yi2.(4) Mean Absolute Percentage Error (MAPE) can be used to measure the prediction results of a model, and its calculation formula is as follows:(24)MAPE=sumfi-yi×100/yin.For the dynamic recurrent prediction model, the determination of the number of hidden layer neurons is an empirical problem that must be constantly tested. The number of hidden layer neurons in the dynamic recurrent network structure, which is not reduced by dimensionality reduction, is set to 2–10. The prediction error proportionality to the number of neurons in different hidden layers is shown in Figure8. In the dynamic recurrent network model, MAPE gradually decreases when the number of hidden layer neurons increases from 2 to 9. When the number of hidden layer neurons increases from 9 to 10, MAPE also gradually increases, indicating that the optimal number of hidden layer neurons was 8. For the Dimensionality Reduction-Dynamic Recurrent model, the optimal number of hidden layer neurons is 5. When the number of hidden layer neurons is larger than 5, the increasing trend of MAPE gradually flattens and eventually stabilizes with the increase in the number of hidden layer neurons, indicating that the stability of the Dimensionality Reduction-Dynamic Recurrent model is strong.Figure 8
MAPE value and the number of neurons.According to the analysis results of the two preceding models, the network structure of dynamic recurrent model is 14-9-1, and the network structure of Dimensionality Reduction-Dynamic Recurrent model is 3-5-1.Four methods were tested 30 times using the sample data with a time step of 1 min to visualize their prediction effect. When the dynamic recurrent network is trained, the transfer function of the hidden layer unit takes the S tangent function Tansig, the output unit also uses the S tangent function Tansig, and the training function adopts the momentum and adaptive gradient decreasing training function, traingdx.An excerpt fragment from the 150 sets of test data is shown in Figure9. Results show that the four prediction modeling methods have achieved good prediction results for the time series of key parameter. However, there is a high correspondence between fitting curve of predicted data and actual real values and have greater prediction accuracy than the three other prediction models.Figure 9
Prediction results of key parameters.Several kinds of evaluation indexes of the four prediction models are calculated. Figure10 indicates the prediction error chart of the four prediction methods, and Table 9 shows the performance indexes of several prediction methods, and the result is presented. The MAE and the RMSE of this method are all excellent among the four methods, which show that this method has higher precision than the other methods in the chemical key parameter prediction. From the mean value of the maximum error and the absolute value of the percentile relative error, the method demonstrates more stable than the other methods in predicting the dynamic time series of the chemical industry.Table 9
Performance evaluation of different prediction models.
Method ME MAE RMSE MAPE (%) TIME (s) BP 74.641 9.982 129.060 1.023% 337.414 ELM 70.414 7.916 99.137 0.793% 0.77 Dynamic Recurrent 49.381 5.768 69.117 0.577% 40.341 Proposed Method 31.360 4.879 55.958 0.498% 34.732Figure 10
Prediction error of four prediction models.Among the four methods of this application, ELM has the shortest prediction time but has weaker prediction accuracy and stability than the proposed method in this paper, and the BP neural network is poor in all aspects. After dimensionality reduction, the performance of this method is superior to the single dynamic recurrent model method.


### 4.2. Evaluation and Discussion

In this application, the risk degree of air separation distillation section is assumed to be affected by several key links, such as air separation tower system, air purification system, and air precooling system. Moreover, the risk degree of air separation tower system is considerably influenced by several key links. The evaluation index and evaluation coefficient are presented in Table10.Table 10
Second-level fuzzy evaluation index and evaluation coefficient.
key production links key parameters evaluation labeling evaluation coefficient Air precooling system Circulation water level Q 11 880  Air separation tower system (contain Distillation column) Liquid air level in lower column Q 21 1110 Distillation column pressure Q 22 1050 Liquid oxygen level in the main condenser Q 23 1010 Distillation column top temperature Q 24 1040  Air purification system Adsorption barrel temperature Q 31 890 Adsorption barrel pressure Q 32 870The second-level fuzzy comprehensive evaluation is used to define the evaluation matrixQ2σ2=[Q21,Q22,Q23,Q24] of the fine air separation distillation section, and the key parameter y(t) exceeding the high alarm value G: (25)C2σ2=yt-GG×100%.The following is only for liquid air level in lower column.Multiply the elements corresponding toQ21 and C21. The second-level fuzzy comprehensive evaluation index matrix of liquid air level in lower column:(26)A21=C21·Q21.The 0–25 period of real-time prediction of the excerpt key parameter is presented in Figure11. The yellow line at 1030 mm is LI101 high alarm value and purple line at 1040 mm is High-High alarm value. The current predictive value of time is 24 min. The corresponding value is 960.236 mm.Figure 11
Real-time prediction and display of key parameters.The effect of the fuzzy comprehensive evaluation of historical key parameter becomes smaller with the increase time. This study considers an algorithm for generating an attenuation factor to accurately evaluate the influence of the fuzzy comprehensive evaluation of historical key parameter on MHIs.Assuming that the second-level fuzzy comprehensive evaluation of key parameter isA, allowing A to assign a timeliness weight tp to the numerical adjustment to eliminate the effect of time factors on the overall MHIs rating results, L is the number of abnormal after the last time the real value of the key parameter exceeds the high alarm.(27)∑p=1Ltp=1,1≤p≤L.Time weighttp is called time attenuation factor. In this paper, the attenuation trend of second-level fuzzy comprehensive evaluations with time is presented by decreasing the ratio of equal ratio, where tp-1≠tp≠tp+1.(28)tp-1tp=tptp+1,2≤p≤L-1.Input. L historical key parameter matrix A21 and corresponding time attenuation factor t1,t2,…,tL.Output. Second-level fuzzy comprehensive evaluation value of eliminating time factor is A′(29)A′=∑p=1LA×tp.The interval L is 7 for the last time the segment data exceeds the high alarm value to the last return to the security value, in which the key parameter exceeds the percentage of the normal threshold. The second-level fuzzy comprehensive evaluation and the attenuation factor are shown in Table11.Table 11
Percentage of key parameter over normal threshold, second-level fuzzy comprehensive evaluation, and attenuation factor.
key parameter C 21 A 21 t p 1037.760 0.753% 8.283 0.038 1035.590 0.543% 5.973 0.055 1038.630 0.838% 9.218 0.078 1038.190 0.795% 8.745 0.112 1042.530 1.217% 13.387 0.160 1044.700 1.427% 15.697 0.229 1042.530 1.217% 13.387 0.327The second-level fuzzy comprehensive evaluation value excluding time factor isA′21, which is calculated as follows:(30)A′21=8.2835.9739.2188.74513.38715.69713.3870.0380.0550.0780.1120.1600.2290.327=12.4558.Then, the weight set of the first layer is determined: the effect of each index on the air separation distillation section risk grade, that is, the weight allocation, is shown in Table12.Table 12
Level weight allocation.
Hazard installation Key production number link Weight(D) Air separation distillation section B Air precooling systemB1 0.07 Air separation tower systemB2 0.76 Air purification systemB3 0.17Multiply the second-level fuzzy comprehensive evaluation matrixA2′ and the weight of the first-level fuzzy comprehensive evaluation, so the risk grade index of the air separation distillation section in CIP is as follows:(31)B2=A2σ2′·D2.We assume that the risk grade index of the air precooling system and the air purification system isB1=4.796,B3=6.163, the second-level fuzzy comprehensive evaluation values A′22, A′23 and A′24 of distillation column pressure, liquid oxygen level in the main condenser, and distillation column top temperature are 10.637, 9.58, and 6.37, respectively, which are the key parameters of air separation tower system; therefore, the risk grade indexR = B 1 × 0.07 + B2×0.76 + B3×0.17 = 4.796×0.07 + (10.637+9.58+6.37+12.4558) × 0.76+6.163×0.17 = 31.056 of the air separation distillation section. Combined with the dynamic evaluation level in Table 1, the dynamic grade of the air separation distillation section in CIP is determined to be 31.056, which is a moderate risk grade. The division of the dynamic risk level fully considers the influence of the historical factors of MHIs on the overall risk, and the historical data will have decreasing impact on the overall risk over time.


## 4.1. Prediction Model

BP neural network is a feedback-free feedforward neural network, which is composed of a set of interconnected neurons, in which each neuron is connected with the corresponding weight [41]. The extreme learning machine (ELM) is a single hidden layer feedforward neural network. It is made up by a simple structure with strong learning and generalization ability. This model has been extensively used in the fields of pattern recognition and regression estimation [42].Through filtered of key parameter correlation variables and dimensionality reduction, 14 groups of strong correlation variables were filtered from 38 groups of correlation variables. After dimensionality reduction, 3 groups of principal components representing 14 groups of strong correlation variables were obtained.The following is the prediction of key parameter. The method is compared with BP neural network, ELM, and dynamic recurrent prediction model without dimensionality reduction. For fairness, the selected variables are also analyzed by visual qualitative analysis and quantitative correlation analysis. The predicted experimental data are derived from 14 sets of solid correlation variables, which remove abnormal values data and acquire 1000 sets of effective data as experimental samples. The distribution of the sample set is presented in Table8.Table 8
Sample set allocation table.
Data sets prorate% amount Training 70% 700 Validation 15% 150 Testing 15% 150In this study, the configuration of hardware involves CPU of Intel Core i5-7300HQ, with main clock frequency of 2.50 GHz and GPU is GeForce GTX1050Ti by NVIDIA. This study adopts the following international evaluation index to measure the prediction efficiency of the model.(1) Maximum Error (ME)(21)ME=max⁡fi-yi,i=1,2,3,…,n.Among them, n is the sample size,fi is the predicted value, and yi is the true value.(2) Mean Absolute Error (MAE) is the average of the absolute value of the deviation between all observed and predicted values.(22)MAE=1n∑i=1i=nfi-yi.(3) Root Mean Square Error (RMSE) represents the discreteness of the error distribution. A small RMSE leads to concentrated error distribution and high prediction accuracy.(23)RMSE=1n∑i=1i=nfi-yi2.(4) Mean Absolute Percentage Error (MAPE) can be used to measure the prediction results of a model, and its calculation formula is as follows:(24)MAPE=sumfi-yi×100/yin.For the dynamic recurrent prediction model, the determination of the number of hidden layer neurons is an empirical problem that must be constantly tested. The number of hidden layer neurons in the dynamic recurrent network structure, which is not reduced by dimensionality reduction, is set to 2–10. The prediction error proportionality to the number of neurons in different hidden layers is shown in Figure8. In the dynamic recurrent network model, MAPE gradually decreases when the number of hidden layer neurons increases from 2 to 9. When the number of hidden layer neurons increases from 9 to 10, MAPE also gradually increases, indicating that the optimal number of hidden layer neurons was 8. For the Dimensionality Reduction-Dynamic Recurrent model, the optimal number of hidden layer neurons is 5. When the number of hidden layer neurons is larger than 5, the increasing trend of MAPE gradually flattens and eventually stabilizes with the increase in the number of hidden layer neurons, indicating that the stability of the Dimensionality Reduction-Dynamic Recurrent model is strong.Figure 8
MAPE value and the number of neurons.According to the analysis results of the two preceding models, the network structure of dynamic recurrent model is 14-9-1, and the network structure of Dimensionality Reduction-Dynamic Recurrent model is 3-5-1.Four methods were tested 30 times using the sample data with a time step of 1 min to visualize their prediction effect. When the dynamic recurrent network is trained, the transfer function of the hidden layer unit takes the S tangent function Tansig, the output unit also uses the S tangent function Tansig, and the training function adopts the momentum and adaptive gradient decreasing training function, traingdx.An excerpt fragment from the 150 sets of test data is shown in Figure9. Results show that the four prediction modeling methods have achieved good prediction results for the time series of key parameter. However, there is a high correspondence between fitting curve of predicted data and actual real values and have greater prediction accuracy than the three other prediction models.Figure 9
Prediction results of key parameters.Several kinds of evaluation indexes of the four prediction models are calculated. Figure10 indicates the prediction error chart of the four prediction methods, and Table 9 shows the performance indexes of several prediction methods, and the result is presented. The MAE and the RMSE of this method are all excellent among the four methods, which show that this method has higher precision than the other methods in the chemical key parameter prediction. From the mean value of the maximum error and the absolute value of the percentile relative error, the method demonstrates more stable than the other methods in predicting the dynamic time series of the chemical industry.Table 9
Performance evaluation of different prediction models.
Method ME MAE RMSE MAPE (%) TIME (s) BP 74.641 9.982 129.060 1.023% 337.414 ELM 70.414 7.916 99.137 0.793% 0.77 Dynamic Recurrent 49.381 5.768 69.117 0.577% 40.341 Proposed Method 31.360 4.879 55.958 0.498% 34.732Figure 10
Prediction error of four prediction models.Among the four methods of this application, ELM has the shortest prediction time but has weaker prediction accuracy and stability than the proposed method in this paper, and the BP neural network is poor in all aspects. After dimensionality reduction, the performance of this method is superior to the single dynamic recurrent model method.


## 4.2. Evaluation and Discussion

In this application, the risk degree of air separation distillation section is assumed to be affected by several key links, such as air separation tower system, air purification system, and air precooling system. Moreover, the risk degree of air separation tower system is considerably influenced by several key links. The evaluation index and evaluation coefficient are presented in Table10.Table 10
Second-level fuzzy evaluation index and evaluation coefficient.
key production links key parameters evaluation labeling evaluation coefficient Air precooling system Circulation water level Q 11 880  Air separation tower system (contain Distillation column) Liquid air level in lower column Q 21 1110 Distillation column pressure Q 22 1050 Liquid oxygen level in the main condenser Q 23 1010 Distillation column top temperature Q 24 1040  Air purification system Adsorption barrel temperature Q 31 890 Adsorption barrel pressure Q 32 870The second-level fuzzy comprehensive evaluation is used to define the evaluation matrixQ2σ2=[Q21,Q22,Q23,Q24] of the fine air separation distillation section, and the key parameter y(t) exceeding the high alarm value G: (25)C2σ2=yt-GG×100%.The following is only for liquid air level in lower column.Multiply the elements corresponding toQ21 and C21. The second-level fuzzy comprehensive evaluation index matrix of liquid air level in lower column:(26)A21=C21·Q21.The 0–25 period of real-time prediction of the excerpt key parameter is presented in Figure11. The yellow line at 1030 mm is LI101 high alarm value and purple line at 1040 mm is High-High alarm value. The current predictive value of time is 24 min. The corresponding value is 960.236 mm.Figure 11
Real-time prediction and display of key parameters.The effect of the fuzzy comprehensive evaluation of historical key parameter becomes smaller with the increase time. This study considers an algorithm for generating an attenuation factor to accurately evaluate the influence of the fuzzy comprehensive evaluation of historical key parameter on MHIs.Assuming that the second-level fuzzy comprehensive evaluation of key parameter isA, allowing A to assign a timeliness weight tp to the numerical adjustment to eliminate the effect of time factors on the overall MHIs rating results, L is the number of abnormal after the last time the real value of the key parameter exceeds the high alarm.(27)∑p=1Ltp=1,1≤p≤L.Time weighttp is called time attenuation factor. In this paper, the attenuation trend of second-level fuzzy comprehensive evaluations with time is presented by decreasing the ratio of equal ratio, where tp-1≠tp≠tp+1.(28)tp-1tp=tptp+1,2≤p≤L-1.Input. L historical key parameter matrix A21 and corresponding time attenuation factor t1,t2,…,tL.Output. Second-level fuzzy comprehensive evaluation value of eliminating time factor is A′(29)A′=∑p=1LA×tp.The interval L is 7 for the last time the segment data exceeds the high alarm value to the last return to the security value, in which the key parameter exceeds the percentage of the normal threshold. The second-level fuzzy comprehensive evaluation and the attenuation factor are shown in Table11.Table 11
Percentage of key parameter over normal threshold, second-level fuzzy comprehensive evaluation, and attenuation factor.
key parameter C 21 A 21 t p 1037.760 0.753% 8.283 0.038 1035.590 0.543% 5.973 0.055 1038.630 0.838% 9.218 0.078 1038.190 0.795% 8.745 0.112 1042.530 1.217% 13.387 0.160 1044.700 1.427% 15.697 0.229 1042.530 1.217% 13.387 0.327The second-level fuzzy comprehensive evaluation value excluding time factor isA′21, which is calculated as follows:(30)A′21=8.2835.9739.2188.74513.38715.69713.3870.0380.0550.0780.1120.1600.2290.327=12.4558.Then, the weight set of the first layer is determined: the effect of each index on the air separation distillation section risk grade, that is, the weight allocation, is shown in Table12.Table 12
Level weight allocation.
Hazard installation Key production number link Weight(D) Air separation distillation section B Air precooling systemB1 0.07 Air separation tower systemB2 0.76 Air purification systemB3 0.17Multiply the second-level fuzzy comprehensive evaluation matrixA2′ and the weight of the first-level fuzzy comprehensive evaluation, so the risk grade index of the air separation distillation section in CIP is as follows:(31)B2=A2σ2′·D2.We assume that the risk grade index of the air precooling system and the air purification system isB1=4.796,B3=6.163, the second-level fuzzy comprehensive evaluation values A′22, A′23 and A′24 of distillation column pressure, liquid oxygen level in the main condenser, and distillation column top temperature are 10.637, 9.58, and 6.37, respectively, which are the key parameters of air separation tower system; therefore, the risk grade indexR = B 1 × 0.07 + B2×0.76 + B3×0.17 = 4.796×0.07 + (10.637+9.58+6.37+12.4558) × 0.76+6.163×0.17 = 31.056 of the air separation distillation section. Combined with the dynamic evaluation level in Table 1, the dynamic grade of the air separation distillation section in CIP is determined to be 31.056, which is a moderate risk grade. The division of the dynamic risk level fully considers the influence of the historical factors of MHIs on the overall risk, and the historical data will have decreasing impact on the overall risk over time.


## 5. Conclusion

This study introduces a dynamic early warning method of MHIs in CIP. Data visual qualitative analysis method and quantitative analysis are conducted to filter the correlation variables. Feature analysis, feature vector acquisition performs the dimensionality reduction of key parameters, effectively improving the performance of dynamic recurrent prediction in the chemical process industry parameters. However, there are still some defects in the current study that we cannot solve. In future research, we need to pay more attention to identification of false alarms and consider human factors and social factors. These research directions will be the key points of our future research.



---

*Source: 6250483-2019-05-20.xml*