Closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling

Abstract

The present invention discloses a closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling, which belongs to the field of industrial process anomaly monitoring and includes the following steps: Step 1, extracting four types of features, namely, static slow features, dynamic slow features, dynamic independent components and residuals, through a typical-core slow feature analysis method; Step 2, designing a multi-feature distribution approximation strategy based on a Gaussian mixture model, approximating the distribution of the four types of features to a Gaussian distribution; Step 3, constructing an evaluation index based on Mahalanobis distance to achieve anomaly monitoring. The present invention solves the problem of decreased monitoring performance caused by the masking of the closed-loop system regulation effect on abnormal information in actual industrial processes, realizes real-time and accurate monitoring of anomalies, and can lay a foundation for the safe and stable operation of actual industrial processes.

Description

Technical Field

[0001] The present invention belongs to the field of industrial process anomaly monitoring, and in particular relates to a closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling. Background Art

[0002] With the deepening development of economic globalization and the continuous improvement of industrial production efficiency, industrial processes are becoming more complex and large-scale. In order to ensure the reliability and stability of the process, its operating status must be monitored. However, the actual industrial process is affected by external factors (feed fluctuations, external environmental changes, etc.) and internal factors (equipment aging, control system feedback, etc.). The process often has a variety of operating conditions. The change of operating conditions causes the process to show non-stationary characteristics, that is, the statistical characteristics of the data change over time. If this non-stationarity cannot be effectively characterized, it will lead to a high false alarm rate. In addition, a large number of closed-loop controllers in industrial processes are used to maintain the stable operation of the process. The regulation of the controller will cover up abnormal information, which will lead to complex changes in both the process state and dynamic characteristics. If the dynamic information of the closed-loop control process cannot be effectively characterized, it will also lead to a high false alarm rate.

[0003] Therefore, in order to solve the problem of decreased monitoring performance caused by the regulatory effect of closed-loop control systems masking abnormal information of industrial processes, it is of great significance to implement abnormal monitoring of closed-loop industrial processes based on multi-feature Gaussian mixture modeling methods. Summary of the invention

[0004] In order to solve the above problems, the present invention proposes a closed-loop industrial process abnormality monitoring method based on multi-feature Gaussian mixture modeling, designs a typical-kernel slow feature analysis method to extract four types of features that are strongly correlated with process dynamics; adopts a multi-feature distribution approximation strategy based on a Gaussian mixture model to approximate the arbitrary distribution of the extracted features; constructs an evaluation index based on the Mahalanobis distance; solves the problem of decreased monitoring performance caused by the masking of abnormal information in the actual industrial process due to the regulation effect of the closed-loop system, can accurately detect abnormal situations and indicate the impact of abnormalities on the operating status of the system, not only conveys important information that traditional methods cannot represent, but also significantly improves the monitoring performance of the closed-loop process.

[0005] The technical solution of the present invention is as follows:

[0006] A closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling specifically comprises the following steps:

[0007] Step 1: Extract four types of features, namely, static slow features, dynamic slow features, dynamic independent components and residuals, through the typical-kernel slow feature analysis method;

[0008] Step 2: Design a multi-feature distribution approximation strategy based on the Gaussian mixture model to approximate the distribution of the four types of features to a Gaussian distribution;

[0009] Step 3: Construct an evaluation index based on Mahalanobis distance to realize anomaly monitoring.

[0010] Furthermore, the specific process of step 1 is as follows:

[0011] Step 1.1: The Tennessee-Eastman process is a closed-loop industrial process. The control variables and process variables in the Tennessee-Eastman process are obtained as samples to construct a data set. , Indicates the number of samples; Indicates the number of variables; Indicates the first variables; Indicates The first variable samples; is the transpose symbol;

[0012] The controlled variables include feed rate of each feed port, total feed rate, compressor recirculation valve, discharge valve, separator filling flow, stripper liquid product flow, stripper water flow valve, reactor cooling water flow, and condenser cooling water flow;

[0013] Process variables include feed to each feed port, total feed, recycle flow, reactor feed rate, reactor pressure, reactor grade, reactor temperature, discharge rate, product separator temperature, product separator liquid level, product separator pressure, product separator bottom flow, stripper grade, stripper pressure, stripper bottom flow, stripper temperature, stripper flow, compressor power, reactor cooling water outlet temperature, separator cooling water outlet temperature, and each component;

[0014] Step 1.2: At the moment, time series expansion is performed to obtain the past vector and future vector, and the formulas are:

[0015] ;

[0016] ;

[0017] in, for The past vector after the time series expansion; for The future vector after the time series expansion; , , They are The first, second, and third rows of the time-lapse vector Row data; , , They are The first, second, and third rows of the future vector Row data;

[0018] Step 1.3: Construct past data matrix and future data matrix , respectively:

[0019] ;

[0020] ;

[0021] in, for No. Column data; for No. Column data; is the length of the time series expansion; is the maximum order of the matrix;

[0022] Step 1.4, the typical-core slow feature analysis method includes typical variable analysis and core slow feature analysis. The objective functions of the typical variable analysis and core slow feature analysis are as follows:

[0023] ;

[0024] ;

[0025] in, , are different slow feature numbers; For the A slow feature; for The first difference of For the A slow feature; for time The minimum first-order difference of ; for and The cross-covariance matrix of ; for The autocovariance matrix of ; for The autocovariance matrix of ;

[0026] Step 1.5, obtaining static slow features and dynamic slow features based on singular value decomposition;

[0027] Step 1.6: Use dynamic kernel independent element analysis to build a model and obtain dynamic independent components and residuals.

[0028] Furthermore, the specific process of step 1.5 is as follows:

[0029] Step 1.5.1: Perform singular value decomposition on the Hankel matrix. The formula is:

[0030] ;

[0031] in, and They are the orthogonal matrices composed of left singular variables and right singular variables in Hankel matrix decomposition; It is the diagonal matrix composed of singular values ​​when the Hankel matrix is ​​decomposed;

[0032] Step 1.5.2, select observation vectors, calculate the transformation matrix corresponding to the past data matrix and the future data matrix, so as to obtain the state vector with strong time series correlation and the residual vector with weak time series correlation. The specific formula is:

[0033] ;

[0034] ;

[0035] ;

[0036] ;

[0037] in, for The transformation matrix corresponding to the past data matrix of the observation vector; for The transformation matrix corresponding to the future data matrix of the observation vector; is the state vector; is the residual vector; for Before column, front Column corresponds to the front observation vector;

[0038] Step 1.5.3. State vector When the objective function of kernel slow feature analysis is met, features with strong temporal correlation and slow change rate are obtained, which are solved by two-step singular value decomposition. is the state vector No. Column, No. Column corresponds to observation vectors; before performing the two-step singular value decomposition, it is necessary to The state vector kernel matrix is ​​obtained by mapping the kernel function to the high-dimensional data space ,right Centralize the state vector to get the kernel matrix ;

[0039] The first step is to decompose the singular value into The covariance matrix of is decomposed into singular values ​​to obtain the data matrix after the first whitening:

[0040] ;

[0041] ;

[0042] in, For When the covariance matrix of is subjected to singular value decomposition, the orthogonal matrix composed of left singular variables is generated; For The diagonal matrix composed of singular values ​​when the covariance matrix of is subjected to singular value decomposition; It is the data matrix after whitening of the final state vector obtained after the first step of singular value decomposition is completed; For The first whitened data matrix is ​​obtained by performing singular value decomposition on the covariance matrix of ;

[0043] The second step is to calculate the singular value decomposition The first difference of ,right The covariance matrix of is subjected to singular value decomposition to obtain the second whitened data matrix:

[0044] ;

[0045] in, For The second whitened data matrix is ​​obtained by performing singular value decomposition on the covariance matrix of ; is the mapping matrix; For The diagonal matrix composed of singular values ​​when the covariance matrix of is subjected to singular value decomposition; is the identity matrix;

[0046] Step 1.5.4: Calculate static features :

[0047] ;

[0048] Step 1.5.5: The singular values ​​in are arranged in order from small to large, with the first The slow features are divided into static features and dynamic features, so the slow features are divided into static slow features and dynamic slow features, and the fast features are divided into static fast features and dynamic fast features, as follows:

[0049] ;

[0050] ;

[0051] in, It is a dynamic feature; and Represent static slow features and dynamic slow features respectively; It is a static fast feature; It is a dynamic and fast feature; for The first difference of for The first difference of for The first-order difference of .

[0052] Furthermore, the specific process of step 1.6 is as follows:

[0053] Step 1.6.1. Residual vector , is the residual vector No. Column, No. Column corresponds to observation vectors, and compare each observation vector with the previous observation vector The observation vectors are augmented to obtain the augmented residual vector data matrix :

[0054] ;

[0055] in, for The residual vector data matrix augmented at each moment; , , They are The first, second, and third columns of the residual vector at the moment List;

[0056] Step 1.6.2: The residual vector kernel matrix is ​​obtained by mapping the kernel function to the high-dimensional data space ,right Perform centralization to obtain the centralized residual vector kernel matrix ,right The characteristic decomposition is expressed as:

[0057] ;

[0058] in, For the feature vectors; For the The eigenvalues ​​corresponding to the eigenvectors;

[0059] Get the intermediate variable after whitening for:

[0060] ;

[0061] in, is the whitened data matrix of the residual vector; for dimension; is the eigenvalue matrix of the residual vector; is the eigenvector matrix composed of all eigenvectors;

[0062] Step 1.6.3: The unmixing matrix is ​​obtained using the fast independent element analysis algorithm. The unmixing matrix of , for The unit orthogonal matrix of order; is the unmixing matrix Column, No. Column corresponds to Finally, we get the dynamic independent component and residual They are:

[0063] ;

[0064] .

[0065] Furthermore, the specific process of step 2 is as follows:

[0066] Step 2.1: Design static slow features The probability density function of the Gaussian mixture model is:

[0067] ;

[0068] in, is the probability density function of the Gaussian mixture model; is the number of Gaussian components; For the The weights of the Gaussian components; For the The mean of the Gaussian components; For the The variance of the Gaussian components;

[0069] Step 2.2: Design the log-likelihood function as:

[0070] ;

[0071] in, is the weight; is the mean; is the variance; Represents the total number of static slow features; Indicates A static slow feature;

[0072] Step 2.3, use the expectation maximization algorithm to solve the parameters of the log-likelihood function, including weight, mean, and variance;

[0073] Step 2.4: After the solution is completed, the static slow characteristics of the Gaussian distribution are obtained , dynamic slow characteristics after Gaussian distribution , dynamic independent components after Gaussian distribution And the residual after Gaussian distribution .

[0074] Furthermore, the specific process of step 2.3 is as follows:

[0075] Step 2.3.1: Given the current parameters, each sample belongs to The probability of a Gaussian component :

[0076] ;

[0077] in, For the A latent variable with 1 Gaussian component;

[0078] Step 2.3.2, maximize the log-likelihood function and update the parameters according to the following three formulas:

[0079] ;

[0080] ;

[0081] ;

[0082] Step 2.3.3, iteratively execute steps 2.3.1 to 2.3.2 until the model converges, after which the required parameters are obtained and the log-likelihood function is calculated;

[0083] Step 2.3.4: Use the Bayesian Information Criterion to introduce a regularization term based on the log-likelihood function to penalize the complexity of the model, as shown in the following formula:

[0084] ;

[0085] in, is the value of the Bayesian Information Criterion.

[0086] Furthermore, the specific process of step 3 is as follows:

[0087] Step 3.1: Based on the idea of ​​Bayesian inference probability, The local probability density under each Gaussian component is:

[0088] ;

[0089] in, For the The local probability density under Gaussian components; for To The Mahalanobis distance of the Gaussian components is greater than To The probability of the Mahalanobis distance of the Gaussian components; For the Gaussian Components , , The set representation of three parameters; is the Mahalanobis distance; For the The static slow speed characteristics after Gaussian distribution corresponding to the Gaussian components;

[0090] Step 3.2 Belong to The posterior probability of the Gaussian components is , combining the local probability density with the posterior probability to obtain the four required evaluation indicators:

[0091] ;

[0092] ;

[0093] ;

[0094] ;

[0095] in, , , , They are the evaluation indicators corresponding to static slow features, dynamic slow features, dynamic independent components and residuals;

[0096] Step 3.3: Use kernel density estimation to calculate the control limits of the four evaluation indicators, which are the control limits of the static slow feature evaluation indicator , control limits of dynamic slow characteristic evaluation index , control limits of dynamic independent component evaluation indicators and the control limits of the residual evaluation index ; Compare the evaluation index with the control limit. When it exceeds the control limit, it is an abnormal situation.

[0097] The beneficial technical effects brought about by the present invention are as follows.

[0098] 1. The present invention extracts features with strong dynamic correlation and slow changes and features with weak dynamic correlation and slow changes through the typical-kernel slow feature analysis method, and extracts nonlinear independent components from features with weak dynamic correlation through dynamic kernel independent element analysis to obtain the most comprehensive dynamic information; for each extracted feature, an approximate approximation strategy based on Gaussian mixture model is used to approximate the arbitrary distribution of the multi-features of the extracted features; and anomaly monitoring is achieved using an evaluation index based on Mahalanobis distance;

[0099] 2. The present invention adopts an abnormal monitoring method based on multi-feature Gaussian mixture modeling to realize abnormal monitoring of closed-loop industrial processes, reduce the impact of the closed-loop controller's regulation on the dynamic characteristics and operating status of the process, realize effective characterization of process dynamic information and accurately identify changes in process operating status, effectively improve abnormal monitoring performance, and ensure the smooth operation of actual industrial processes. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 This is a closed-loop industrial process anomaly monitoring framework diagram based on multi-feature Gaussian mixture modeling of the present invention.

[0101] Figure 2 It is a monitoring result diagram of the static slow characteristic evaluation index for the random variable fault of the Tennessee-Eastman process in the experiment of the present invention.

[0102] Figure 3 It is a monitoring result diagram of dynamic slow characteristic evaluation index for random variable fault of Tennessee-Eastman process in the experiment of the present invention.

[0103] Figure 4 It is a monitoring result diagram of dynamic independent component evaluation index for random variable fault of Tennessee-Eastman process in the experiment of the present invention.

[0104] Figure 5It is a graph of residual evaluation index monitoring results for random variable faults in the Tennessee-Eastman process in the experiment of the present invention. DETAILED DESCRIPTION

[0105] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments:

[0106] The present invention proposes a closed-loop industrial process abnormality monitoring method based on multi-feature Gaussian mixture modeling. The method extracts the features with strong dynamic correlation and slow change speed in the process through the typical-kernel slow feature analysis method, so as to effectively avoid the influence of static information on the process dynamics; then, a multi-feature approximation strategy based on the Gaussian mixture model is adopted to approximate the arbitrary distribution of features, and an evaluation index based on the Mahalanobis distance is established; the problem that the nonlinear relationship between variables in the actual industrial process and the abnormal information are masked under the regulation of the closed-loop control system, resulting in the decline of abnormal monitoring performance, is solved. The proposed method can detect abnormal situations in a timely and accurate manner, and indicate the influence of abnormalities on the operating status of the system. It not only transmits important information that traditional methods cannot cover, but also significantly improves the monitoring performance of the closed-loop process.

[0107] like Figure 1 As shown in the figure, a closed-loop industrial process abnormality monitoring method based on multi-feature Gaussian mixture modeling is shown, such as the process related to the condenser in the Tennessee-Eastman process, in which the controller, actuators such as valves, controlled objects such as reactors and condensers, and detection devices such as temperature, flow and pressure sensors together constitute a closed-loop industrial process. For similar closed-loop industrial processes of this type, abnormality monitoring is performed through a strategy based on typical-kernel slow feature analysis combined with Gaussian mixture model, which specifically includes the following steps:

[0108] Step 1: Extract four types of features: static slow features, dynamic slow features, dynamic independent components and residuals through the typical-kernel slow feature analysis method; the specific process is as follows:

[0109] Step 1.1: Obtain the control variables and process variables of the Tennessee-Eastman process as samples and construct a data set , Indicates the number of samples; Indicates the number of variables; Indicates the first variables; Indicates The first variable samples; is the transpose symbol;

[0110] The controlled variables include feed rate at port A, feed rate at port D, feed rate at port E, total feed rate, compressor recirculation valve, discharge valve, separator filling flow, stripper liquid product flow, stripper water flow valve, reactor cooling water flow, and condenser cooling water flow;

[0111] Process variables include feed at port A, feed at port D, feed at port E, total feed, recycle flow, reactor feed rate, reactor pressure, reactor grade, reactor temperature, discharge rate, product separator temperature, product separator liquid level, product separator pressure, product separator bottom flow, stripper grade, stripper pressure, stripper bottom flow, stripper temperature, stripper flow, compressor power, reactor cooling water outlet temperature, separator cooling water outlet temperature, and individual components.

[0112] Step 1.2: At the moment, time series expansion is performed to obtain the past vector and future vector, and the formulas are:

[0113] ;

[0114] ;

[0115] in, for The past vector after the time series expansion; for The future vector after the time series expansion; , , They are The first, second, and third rows of the time-lapse vector Row data; , , They are The first, second, and third rows of the future vector Row data;

[0116] Step 1.3 and The mean values ​​of and ,Will and After removing the mean, they are and , construct the past data matrix and future data matrix , respectively:

[0117] ;

[0118] ;

[0119] in, for No. Column data; for No. Column data; is the length of the time series expansion, which can be determined based on the autocorrelation of the variable ; is the maximum order of the matrix; in the present invention, ;

[0120] Step 1.4, the typical-core slow feature analysis method includes typical variable analysis and core slow feature analysis. The objective functions of the typical variable analysis and core slow feature analysis are as follows:

[0121] ;

[0122] ;

[0123] in, , are different slow feature numbers; For the A slow feature; for The first difference of For the A slow feature; for time The minimum first-order difference of ; for and The cross-covariance matrix of ; for The autocovariance matrix of ; for The autocovariance matrix of ;

[0124] Step 1.5: Obtain static slow features and dynamic slow features based on singular value decomposition; the specific process is as follows:

[0125] Step 1.5.1: Solve the objective function of canonical variable analysis to obtain the characteristics of strong time series correlation. Perform singular value decomposition on the Hankel matrix to maximize the correlation between the two data matrices. The formula is:

[0126] ;

[0127] in, and They are the orthogonal matrices composed of left singular variables and right singular variables in Hankel matrix decomposition; It is the diagonal matrix composed of singular values ​​when the Hankel matrix is ​​decomposed;

[0128] Step 1.5.2, select observation vectors, calculate the transformation matrix corresponding to the past data matrix and the future data matrix, so as to obtain the state vector with strong time series correlation and the residual vector with weak time series correlation. The specific formula is:

[0129] ;

[0130] ;

[0131] ;

[0132] ;

[0133] in, for The transformation matrix corresponding to the past data matrix of the observation vector; for The transformation matrix corresponding to the future data matrix of the observation vector; is the state vector; is the residual vector; for Before column, front Column corresponds to the front observation vector, The number of The size of the singular value in is determined, and the corresponding larger singular value is retained. In the embodiment of the present invention, ;

[0134] Step 1.5.3. State vector When the objective function of kernel slow feature analysis is met, features with strong temporal correlation and slow change rate are obtained, which can be solved by two-step singular value decomposition. is the state vector No. Column, No. Column corresponds to observation vectors; before performing the two-step singular value decomposition, it is necessary to The state vector kernel matrix is ​​obtained by mapping the kernel function to the high-dimensional data space ,right Centralize the state vector to get the kernel matrix ;

[0135] The first step is to decompose the singular value into The covariance matrix of is decomposed into singular values ​​to obtain the data matrix after the first whitening:

[0136] ;

[0137] ;

[0138] in, For When the covariance matrix of is subjected to singular value decomposition, the orthogonal matrix composed of left singular variables is generated; For The diagonal matrix composed of singular values ​​when the covariance matrix of is subjected to singular value decomposition; It is the data matrix after whitening of the final state vector obtained after the first step of singular value decomposition is completed; For The first whitened data matrix is ​​obtained by performing singular value decomposition on the covariance matrix of ;

[0139] The second step is to calculate the singular value decomposition The first difference of ,right The covariance matrix of is subjected to singular value decomposition to obtain the second whitened data matrix:

[0140] ;

[0141] in, For The second whitened data matrix is ​​obtained by performing singular value decomposition on the covariance matrix of ; is the mapping matrix; For The diagonal matrix composed of singular values ​​when the covariance matrix of is subjected to singular value decomposition; is the identity matrix;

[0142] Step 1.5.4: Calculate static features :

[0143] ;

[0144] Step 1.5.5: The singular values ​​in the are arranged in order from small to large. The smaller the singular value, the slower the change speed of the corresponding feature. According to the change speed of the feature, the slow feature and the fast feature can be arranged in front. The slow feature is shown in the example of the present invention. ; Features are divided into static features and dynamic features, so slow features are divided into static slow features and dynamic slow features, and fast features are divided into static fast features and dynamic fast features, as follows:

[0145] ;

[0146] ;

[0147] in, It is a dynamic feature; and Represent static slow features and dynamic slow features respectively; It is a static fast feature; It is a dynamic and fast feature; for The first difference of for The first difference of for The first difference of

[0148] Step 1.6: Traditional canonical variable analysis believes that the time series correlation of the residual vector is weak, but there is still effective dynamic information of the process in this part of the vector. Therefore, dynamic kernel independent element analysis is used to further model this part of the vector to obtain dynamic independent components and residuals; the specific process is:

[0149] Step 1.6.1. Residual vector , is the residual vector No. Column, No. Column corresponds to observation vectors, first compare each observation vector with the previous observation vector The observation vectors are augmented to obtain the augmented residual vector data matrix :

[0150] ;

[0151] in, for The residual vector data matrix augmented at each moment; The number can be selected according to the dynamic nature of the process, and 1 or 2 is generally selected based on empirical values; , , They are The first, second, and third columns of the residual vector at the moment List;

[0152] Step 1.6.2: The residual vector kernel matrix is ​​obtained by mapping the kernel function to the high-dimensional data space ,right Perform centralization to obtain the centralized residual vector kernel matrix ,right The eigendecomposition can be expressed as:

[0153] ;

[0154] in, For the feature vectors; For the The eigenvalues ​​corresponding to the eigenvectors;

[0155] Get the intermediate variable after whitening for:

[0156] ;

[0157] in, Whiten the data matrix for the residual vector; for dimension; is the eigenvalue matrix of the residual vector; is the eigenvector matrix consisting of all eigenvectors, which is a matrix with dimension The eigenvector matrix of the corresponding residual vector ; represents a diagonal matrix;

[0158] Step 1.6.3: The unmixing matrix is ​​obtained by using the fast independent element analysis algorithm, and its order is ,therefore, The unmixing matrix of , for The unit orthogonal matrix of order; is the unmixing matrix Column, No. Column corresponds to Finally, we get the dynamic independent component and residual They are:

[0159] ;

[0160] ;

[0161] In summary, through the method based on typical-kernel slow feature analysis, four types of features, namely static slow features, dynamic slow features, dynamic independent components and residuals, are successfully extracted; these features effectively capture the essential information of the process, and avoid the interference of static information on the dynamic characteristics of the process through low-dimensional representation, thereby maximizing the accuracy of dynamic feature extraction.

[0162] Step 2: Design a multi-feature distribution approximation strategy based on a Gaussian mixture model to approximate the distribution of the four types of features to a Gaussian distribution. In actual industrial processes, the features extracted by typical-kernel slow feature analysis and dynamic kernel independent element analysis usually do not obey the Gaussian distribution. In order to reduce the impact of non-Gaussian characteristics on monitoring results, a Gaussian mixture model is used to approximate features with arbitrary distribution. The specific process is as follows:

[0163] Step 2.1: Design static slow features The probability density function of the Gaussian mixture model is:

[0164] ;

[0165] in, is the probability density function of the Gaussian mixture model; is the number of Gaussian components; For the The weights of the Gaussian components need to satisfy and ; For the The mean of the Gaussian components; For the The variance of the Gaussian components;

[0166] Step 2.2: In order to estimate the parameters in the model, use maximum likelihood estimation to assist in solving the problem. Introduce a latent variable ,definition For the The latent variable with Gaussian components is The probability density function of , therefore, the conditional probability at the sample point is Probability model with Gaussian components The log-likelihood function is expressed as:

[0167] ;

[0168] in, is the weight; is the mean; is the variance; Represents the total number of static slow features; Indicates A static slow feature;

[0169] Step 2.3: Use the expectation maximization algorithm to solve the parameters of the log-likelihood function, including weight, mean, and variance. The specific process is as follows:

[0170] Step 2.3.1: Given the current parameters, each sample belongs to The probability of a Gaussian component :

[0171] ;

[0172] Step 2.3.2, maximize the log-likelihood function and update the parameters according to the following three formulas:

[0173] ;

[0174] ;

[0175] ;

[0176] Step 2.3.3, iteratively execute steps 2.3.1 to 2.3.2 until the model converges. After convergence, the required parameters can be obtained and the log-likelihood function can be calculated;

[0177] Step 2.3.4: To prevent overfitting caused by too many Gaussian components, the Bayesian Information Criterion is used. The core idea is to introduce a regularization term based on the log-likelihood function to penalize the complexity of the model, as shown in the following formula:

[0178] ;

[0179] in, is the value of the Bayesian Information Criterion.

[0180] Step 2.4: Use the multi-feature distribution approximation strategy based on Gaussian mixture model to transform the static slow features , dynamic slow features , dynamic independent component and residual The arbitrary distribution of the four types of features is approximately Gaussian distribution, and the static slow features after Gaussian distribution are obtained , dynamic slow characteristics after Gaussian distribution , dynamic independent components after Gaussian distribution And the residual after Gaussian distribution ;

[0181] Step 3: Construct an evaluation index based on Mahalanobis distance to realize abnormality monitoring; the specific process is as follows:

[0182] Step 3.1: Based on the idea of ​​Bayesian inference probability, Belong to The posterior probability of the Gaussian components is , the local probability density of the sample under each Gaussian component is:

[0183] ;

[0184] in, For the The local probability density under Gaussian components; for To The Mahalanobis distance of the Gaussian components is greater than To The probability of the Mahalanobis distance of the Gaussian components; For the Gaussian Components , , The set representation of three parameters; For samples to The Mahalanobis distance of the Gaussian components; is the Mahalanobis distance; For the The static slow speed characteristics after Gaussian distribution corresponding to the Gaussian components; is the covariance matrix;

[0185] Step 3.2: Combine the local probability density with the posterior probability to obtain the four required evaluation indicators:

[0186] ;

[0187] ;

[0188] ;

[0189] ;

[0190] in, , , , They are the evaluation indicators corresponding to static slow features, dynamic slow features, dynamic independent components and residuals;

[0191] Step 3.3: Use kernel density estimation to calculate the control limits of the four evaluation indicators, which are the control limits of the static slow feature evaluation indicator , control limits of dynamic slow characteristic evaluation index , control limits of dynamic independent component evaluation indicators and the control limits of the residual evaluation index ; Compare the evaluation index with the control limit. When it exceeds the control limit, it is an abnormal situation.

[0192] Aiming at the problem that abnormal information in actual industrial processes is masked by the closed-loop system regulation, resulting in the degradation of monitoring performance, an abnormal monitoring method based on multi-feature Gaussian mixture modeling is proposed. The specific process is as follows: Figure 1This method first extracts static slow features, dynamic slow features, dynamic independent components and residuals through the transformation matrix, mapping matrix and unmixing matrix in the typical-kernel slow feature analysis; secondly, a multi-feature distribution approximation strategy based on the Gaussian mixture model is designed to obtain the static slow features, dynamic slow features, dynamic independent components and residuals after Gaussian distribution; then, an evaluation index based on the Mahalanobis distance is constructed. , , , and use the kernel density estimation method to calculate the control limits , , and ; Ultimately, real-time and accurate monitoring of anomalies is achieved.

[0193] In order to prove the feasibility and superiority of the present invention, the Tennessee-Eastman process data under closed-loop control was collected, including 11 control variables and 41 process variables. These variables and their descriptions are shown in Table 1:

[0194] Table 1 Variables and their descriptions

[0195] .

[0196] The training data uses normal data of the Tennessee-Eastman process, which contains 500 samples; the test data set contains 960 samples, and the fault begins to appear from the 161st sample. The abnormal monitoring results of the random variable fault of the test data set are as follows: Figure 2-Figure 5 shown.

[0197] Figure 2 Demonstrated static slow feature evaluation metrics for random variable faults Abnormal monitoring results show that the evaluation indicators The control limit is exceeded first, and its fault detection delay is 23 samples. Figure 3 Demonstrated dynamic slow feature evaluation metrics for random variable faults Abnormal monitoring results show that the evaluation indicators Out of control limit, significantly deviates from preset state, fault monitoring delay is 8 samples. Figure 4 Demonstrated dynamic independent component evaluation metrics for random variable faults Abnormal monitoring results show that The fault detection delay is 9 samples. Figure 5 Demonstrates residual evaluation metrics for random variable faults Abnormal monitoring results show that The fault detection delay is 28 samples.

[0198] from Figure 2-Figure 5The results show that the adjustment of the closed-loop controller cannot restore the process to normal or stabilize it to another state, so that the essential state of the process is affected. The proposed closed-loop industrial process abnormality monitoring method based on multi-feature Gaussian mixture modeling can clearly indicate the changes in the process operating state under various faults and the regulation effect of the closed-loop controller on the process operating state, and can effectively improve the closed-loop process monitoring performance.

[0199] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling, characterized in that: The specific steps include: Step 1: Extract four types of features, namely, static slow features, dynamic slow features, dynamic independent components and residuals, through the typical-kernel slow feature analysis method; Step 2: Design a multi-feature distribution approximation strategy based on the Gaussian mixture model to approximate the distribution of the four types of features to a Gaussian distribution; Step 3: Construct an evaluation index based on Mahalanobis distance to realize abnormality monitoring; the specific process is as follows: Step 3.1: Based on the Bayesian inference probability, s′ slow The local probability density under each Gaussian component is: in, is the local probability density under the h-th Gaussian component; s′ slow The Mahalanobis distance to the hth Gaussian component is greater than s′ slow,h The probability of the Mahalanobis distance to the h-th Gaussian component; C h is the hth Gaussian component π h , μ h ,∑ h The set of three parameters represents: D(·) is the Mahalanobis distance; s′ slow,h is the static slow feature after Gaussian distribution corresponding to the hth Gaussian component; Step 3.2, s′ slow The posterior probability belonging to the hth Gaussian component is P post (C h |s′ slow ), combining the local probability density with the posterior probability to obtain the four required evaluation indicators: Among them, BIP ss 、BIP ds , BIP SPE They are the evaluation indicators corresponding to static slow features, dynamic slow features, dynamic independent components and residuals; Step 3.3: Use kernel density estimation to calculate the control limits of the four evaluation indicators, namely, the control limit of the static slow characteristic evaluation indicator BIP ctr,ss , control limit BIP of dynamic slow characteristic evaluation index ctr,ds , control limits of dynamic independent component evaluation indicators and the control limits of the residual evaluation index BIP ctr,SPE ; Compare the evaluation index with the control limit. When it exceeds the control limit, it is an abnormal situation.

2. The closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling according to claim 1 is characterized in that: The specific process of step 1 is as follows: Step 1.1: The Tennessee-Eastman process is a closed-loop industrial process. The control variables and process variables in the Tennessee-Eastman process are obtained as samples, and the data set Y(N×J)=[y1,y2,…,y J ], N represents the number of samples; J represents the number of variables; y J =[y 1,J ,y 2,J ,…,y N,J ] T represents the Jth variable in the data set; y N,J represents the Nth sample in the Jth variable; T is the transposition symbol; The controlled variables include feed rate of each feed port, total feed rate, compressor recirculation valve, discharge valve, separator filling flow, stripper liquid product flow, stripper water flow valve, reactor cooling water flow, and condenser cooling water flow; Process variables include feed to each feed port, total feed, recycle flow, reactor feed rate, reactor pressure, reactor grade, reactor temperature, discharge rate, product separator temperature, product separator liquid level, product separator pressure, product separator bottom flow, stripper grade, stripper pressure, stripper bottom flow, stripper temperature, stripper flow, compressor power, reactor cooling water outlet temperature, separator cooling water outlet temperature, and each component; Step 1.2: At time t, perform time series expansion to obtain the past vector and future vector. The formulas are: and past,t =[and t-1 ,and t-2 ,…,and t-q ] T ; and future,t =[and t ,and t+1 ,…,and t+q-1 ] T ; Among them, y past,t is the past vector after time series expansion at time t; yf uture,t is the future vector after time series expansion at time t; y t-1 ,y t-2 ,y t-q are the first, second, and qth row of data of the vector at time t; t ,y t+1 ,y t+q-1 They are the first, second and qth row of data of the future vector at time t respectively; Step 1.3: Construct past data matrix Y past and the future data matrix Y future , respectively: in, Y past The q+Qth column data; Y future The q+Qth column data; q is the length of the time series expansion; Q is the maximum order of the matrix; Step 1.4, the typical-core slow feature analysis method includes typical variable analysis and core slow feature analysis. The objective functions of the typical variable analysis and core slow feature analysis are as follows: Among them, i and j are different slow feature numbers; s j is the jth slow feature; For j The first difference of i is the ith slow feature; minΔs j is time t j The minimum first-order difference of ; Y past and Y future The cross-covariance matrix of ; Y past The autocovariance matrix of ; Y future The autocovariance matrix of ; Step 1.5, obtaining static slow features and dynamic slow features based on singular value decomposition; Step 1.6: Use dynamic kernel independent element analysis to build a model and obtain dynamic independent components and residuals.

3. The closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling according to claim 2 is characterized in that: The specific process of step 1.5 is as follows: Step 1.5.1: Perform singular value decomposition on the Hankel matrix. The formula is: Among them, U Hankel and V Hankel are orthogonal matrices composed of left singular variables and right singular variables when the Hankel matrix is ​​decomposed; Λ Hankel It is the diagonal matrix composed of singular values ​​when the Hankel matrix is ​​decomposed; Step 1.5.2, select r observation vectors, calculate the transformation matrix corresponding to the past data matrix and the future data matrix, so as to obtain the state vector with strong time series correlation and the residual vector with weak time series correlation. The specific formula is: E=L r Y past ; in, is the transformation matrix corresponding to the past data matrix of r observation vectors; L r is the transformation matrix corresponding to the future data matrix of r observation vectors; Z is the state vector; E is the residual vector; U Hankel,r For U Hankel The first r columns of , the first r columns correspond to the first r observation vectors; Step 1.5.

3. State vector When the objective function of kernel slow feature analysis is met, features with strong time series correlation and slow change speed are obtained, which are solved by two-step singular value decomposition. r is the rth column of the state vector Z, and the rth column corresponds to the rth observation vector; before performing the two-step singular value decomposition, it is necessary to map Z to the high-dimensional data space through the kernel function to obtain the state vector kernel matrix K Z , for K Z Centralize the state vector to get the kernel matrix The first step is to decompose the singular value into The covariance matrix of is decomposed into singular values ​​to obtain the data matrix after the first whitening: in, For When the covariance matrix of is subjected to singular value decomposition, the orthogonal matrix composed of left singular variables is generated; For The diagonal matrix composed of singular values ​​when the covariance matrix of is subjected to singular value decomposition; Γ is the data matrix after the whitening of the final state vector obtained after the first step of singular value decomposition is completed; For The first whitened data matrix is ​​obtained by performing singular value decomposition on the covariance matrix of ; The second step is to calculate the first-order difference of Γ right The covariance matrix of is subjected to singular value decomposition to obtain the second whitened data matrix: in, For The second whitened data matrix is ​​obtained by performing singular value decomposition on the covariance matrix of ; W is the mapping matrix; For The diagonal matrix composed of singular values ​​when the covariance matrix of is subjected to singular value decomposition; I is the unit matrix; Step 1.5.4, calculate the static characteristic s: Step 1.5.5: The singular values ​​in are arranged in ascending order, and the first w are slow features; features are divided into static features and dynamic features, so slow features are divided into static slow features and dynamic slow features, and fast features are divided into static fast features and dynamic fast features, as follows: s=[s slow ,s fast ] T ; Among them, d is the dynamic feature; s slow and d slow Represent static slow features and dynamic slow features respectively; s fast is a static fast feature; d fast It is a dynamic and fast feature; is the first-order difference of s; For slow The first difference of For fast The first-order difference of .

4. The closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling according to claim 3 is characterized in that: The specific process of step 1.6 is as follows: Step 1.6.1, residual vector E = [e1, e2, ..., e r ],e r is the rth column of the residual vector E, the rth column corresponds to the rth observation vector, and each observation vector is augmented with the b observation vectors before each observation vector to obtain the augmented residual vector data matrix E aug : E aug (t)=[e1(t),e2(t-1),…,e r (t-b)]; Among them, E aug (t) is the augmented residual vector data matrix at time t; e1(t), e2(t-1), e r (tb) are the first, second, and rth columns of the residual vector at time t; Step 1.6.2: E aug The residual vector kernel matrix K is obtained by mapping the kernel function to the high-dimensional data space. E , for K E Perform centralization to obtain the centralized residual vector kernel matrix right The characteristic decomposition is expressed as: Among them, β l is the lth eigenvector; l is the eigenvalue corresponding to the lth eigenvector; The intermediate variable u after whitening is: Among them, Φ E is the whitened data matrix of the residual vector; L is the dimension; is the eigenvalue matrix of the residual vector; is the eigenvector matrix composed of all eigenvectors; Step 1.6.3, use the fast independent element analysis algorithm to obtain the unmixing matrix for u, The unmixing matrix of for The unit orthogonal matrix of order; is the unmixing matrix Column, No. Column corresponds to Finally, the dynamic independent component S and residual R are obtained as follows:

5. The closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling according to claim 4 is characterized in that: The specific process of step 2 is: Step 2.1: Design static slow features slow The probability density function of the Gaussian mixture model is: Where p(·) is the probability density function of the Gaussian mixture model; H is the number of Gaussian components; π h is the weight of the hth Gaussian component; μ h is the mean of the hth Gaussian component; ∑ h is the variance of the hth Gaussian component; Step 2.2: Design the log-likelihood function as: Among them, π is the weight; μ is the mean; ∑ is the variance; Represents the total number of static slow features; Indicates A static slow feature; Step 2.3, use the expectation maximization algorithm to solve the parameters of the log-likelihood function, including weight, mean, and variance; Step 2.4: After the solution is completed, the static slow feature s′ after Gaussian distribution is obtained slow , dynamic slow feature d′ after Gaussian distribution slow , dynamic independent component S′ after Gaussian distribution and residual R′ after Gaussian distribution.

6. The closed-loop industrial process anomaly monitoring method based on multi-feature Gaussian mixture modeling according to claim 5 is characterized in that: The specific process of step 2.3 is as follows: Step 2.3.1, given the current parameters, the probability γ(·) that each sample belongs to the hth Gaussian component: Among them, v h is the latent variable of the h-th Gaussian component; Step 2.3.2, maximize the log-likelihood function and update the parameters according to the following three formulas: Step 2.3.3, iteratively execute steps 2.3.1 to 2.3.2 until the model converges, after which the required parameters are obtained and the log-likelihood function is calculated; Step 2.3.4: Use the Bayesian Information Criterion to introduce a regularization term based on the log-likelihood function to penalize the complexity of the model, as shown in the following formula: BIC=-2lnp(s slow |π,μ,∑)+3ln(N); Among them, BIC is the value of the Bayesian Information Criterion.

Images