Data-driven distributed fault detection method and device

Abstract

The invention discloses a data-driven distributed fault detection method and a data-driven distributed fault detection device, which are used for analyzing and extracting more valuable information from two levels of global variables and local sub-block variables to implement distributed fault detection. The distributed fault detection method disclosed by the invention not only relates to the application of a brand new algorithm of manifold regularization slow feature analysis, but also provides a plurality of improved technical schemes related to variable quantum block division, neighborhood and reconstruction modes thereof, and whitening processing. In order to overcome the defect that feature analysis is only carried out on each sub-block after blocking is completed in the traditional method, the method provided by the invention inherits the modeling thought of mutual staggering from local sub-blocks to global sub-blocks and then to the local sub-blocks, and comprehensively represents the time sequence and spatial distribution features of learning data changes. Besides, the invention further designs a distributed fault detection device based on the same concept, and the central control module can execute the execution process of a manifold regularization slow feature analysis algorithm and also can autonomously divide variable sub-blocks.

Description

Technical Field

[0001] The present invention belongs to the technical field of data-driven fault detection, and particularly relates to a data-driven distributed fault detection method and device. Background Art

[0002] With the development of modern industry, the scale of industrial production is constantly expanding, and the number of operating units and the complexity of each unit are increasing continuously, which poses a great challenge to the establishment of accurate mechanism models. At the same time, with the wide application of various advanced sensors and industrial Internet of Things technologies, a large amount of data collected in the production process has laid a solid data foundation for the implementation of data-driven fault detection. However, considering the trend of large-scale modern industrial processes, traditional centralized data-driven fault detection methods cannot accurately detect faults in local units. To solve this problem, distributed fault detection methods have gradually received extensive attention in the industrial community, and more and more engineers prefer to implement distributed fault monitoring.

[0003] The Chinese invention patent with the application number / patent number 202211396135.5 proposes to establish a fault detection model for each sub-block by using distributed canonical correlation analysis, and then use the average consensus algorithm to implement distributed fault detection fusion. The Chinese invention patent with the application number / patent number 202310230809.2 discloses a plant-level distributed process monitoring method based on knowledge and data dual drive. By implementing variable partitioning through a knowledge causal graph and then introducing standardized mutual information to describe the data correlation within the sub-blocks, fault detection from local to global is achieved. The Chinese invention patent with the application number / patent number 202311381737.8 discloses a distributed data-driven fault detection method and system, which applies an improved average consensus algorithm to achieve distributed fault detection for the finishing mill unit in the rolling process. The methods and techniques disclosed in these patents basically first implement variable partitioning, then establish sub-block models respectively, and finally perform information fusion of the sub-block fault detection results. The relevance between each sub-block cannot be explored during the modeling process.

[0004] In addition, papers on data-driven fault detection also favor distributed implementation methods. For example, a paper titled "Early Fault Detection in the Whole-Plant Process Based on Multi-Subspace Weighted Moving Window Principal Component Analysis" published in the journal "Journal of Zhejiang University (Engineering Science)" in 2024 proposed a two-layer variable block partitioning method that combines process knowledge and data, and then used weighted moving window principal component analysis to establish a fault detection model for each sub-block. A paper titled "Distributed Dynamic Graph Embedding for Quality-Related Monitoring in the Flotation Process" published in the journal "IEEE Transactions on Instrumentation and Measurement" in 2024 divided variable blocks related and unrelated to product quality, and then used a distributed dynamic graph embedding analysis algorithm to establish models for each sub-block to implement fault detection for the industrial flotation process. The methods proposed in these papers also first partition variables into blocks, and then establish sub-block models to mine and analyze the change characteristics of data in each sub-block, and none of them further consider the correlation between the number of sub-blocks.

[0005] Generally speaking, to implement distributed fault detection, first, multiple variable blocks need to be partitioned. Second, it is about how to represent and learn the data change characteristics of each variable block. Finally, it is the fusion of the fault detection results of each sub-block. The partitioning of variable blocks can be achieved through process organizational structure knowledge or process data, and sub-block modeling to extract data change characteristics is the key to fault detection. Considering that data temporal change and spatial distribution change are two main ways to describe data change characteristics, when implementing feature representation learning for each variable sub-block, the coordination from both temporal and spatial distribution perspectives should be considered. In addition, as a whole, for large-scale process objects, how to jointly consider the relationship between global variables of each sub-block as a whole when establishing local models for sub-blocks is very helpful for improving the sensitivity of fault detection. Summary of the Invention

[0006] The main technical problem to be solved by the present invention is: how to comprehensively represent and learn the temporal and spatial distribution characteristics of data change from the two levels of global variables and local sub-block variables, so as to analyze and extract more valuable information for implementing distributed fault detection.

[0007] The technical solution adopted by the present invention to solve the above problems is: a data-driven distributed fault detection method and device; wherein, the implementation process of a data-driven distributed fault detection method disclosed by the present invention includes the following steps 1 to 5.

[0008] Step 1: Obtain the dataset when the industrial process is operating normally and perform standardization processing on it to obtain the training data matrix X. Among them, each column vector in X is composed of the measurement data corresponding to M variables.

[0009] Step 2: Divide the M variables into B non-overlapping variable sub-blocks, and correspondingly divide X into the corresponding B sub-block data matrices X1, X2, …, X B .

[0010] Step 3: Use the manifold regularization slow feature analysis algorithm to perform feature analysis on X1, X2, …, X B to obtain the corresponding sub-block decomposition matrices W1, W2, …, W B and sub-block loading matrices P1, P2, …, P B .

[0011] Step 4: Calculate the corresponding fault detection index for each column vector in X, and then determine the upper limit of the change of the fault detection index

[0012] Step 5: Obtain a set of data at the latest sampling time t of the industrial process and perform standardization processing on it to obtain the column vector x t After that, calculate the corresponding fault detection index Then, according to whether it exceeds to judge whether the industrial process is malfunctioning.

[0013] In the above implementation steps, the core technology of the method of the present invention lies in the process of using the manifold regularization slow feature analysis algorithm to perform feature analysis on the B sub-block data matrices in Step 3. This process includes the following Steps 3-1 to 3-5.

[0014] Step 3-1: Determine the corresponding neighborhood for each column vector in X b to obtain the corresponding neighborhood reconstruction matrix Φ b of X b , and then calculate the difference matrix b of X where, X b represents the b-th sub-block data matrix, and the number b = 1, 2, …, B.

[0015] Step 3-2: Perform whitening processing on X b to obtain the whitening matrix D b After that, calculate the symmetric matrix G through the formula and perform singular value decomposition on G, and denote the decomposition vector corresponding to the smallest singular value as w b; where \(I\) represents the identity matrix, the superscript \(T\) represents the transpose of a matrix or vector, and \(\eta\) is the regularization coefficient.

[0016] Step 3-3: After combining the \(B\) decomposed vectors \(w_1, w_2, \ldots, w\) B into the global transformation vector \(w\), and combining the \(B\) sub-block data matrices \(X_1, X_2, \ldots, X\) B into the global data matrix \(X_0\), then calculate the global load vector \(p\).

[0017] Step 3-4: Corresponding to \(p\), divide it into \(B\) load vectors \(p_1, p_2, \ldots, p\) B , and then through update each sub-block data matrix, and return to Step 3-2 until each sub-block data matrix corresponds to \(K\) decomposed vectors and \(K\) load vectors; where \(w\) b and \(p\) b respectively represent the decomposed vector and load vector corresponding to the \(b\)-th sub-block data matrix.

[0018] Step 3-5: Combine the \(K\) decomposed vectors and \(K\) load vectors corresponding to each sub-block data matrix into sub-block decomposition matrices \(W_1, W_2, \ldots, W\) B and sub-block load matrices \(P_1, P_2, \ldots, P\) B .

[0019] The calculation process of the fault detection index involved in Step 3 and Step 6 includes the following two steps.

[0020] Step A1: According to the \(B\) variable sub-blocks in Step 2, divide any given column vector \(x\) into the corresponding \(B\) sub-block column vectors \(x_1, x_2, \ldots, x\) B and then through and calculate the sub-block eigenvector \(s\) b and sub-block error vector \(\varepsilon\) b in sequence; where \(W\) b and \(P\) b respectively represent the \(b\)-th sub-block decomposition matrix and the \(b\)-th sub-block load matrix.

[0021] Step A2: Calculate the \(b\)-th sub-block detection index through , and then set the fault monitoring index corresponding to \(x\) where represents the upper limit of the change of the \(b\)-th sub-block detection index.

[0022] It should be further noted that the implementation process of the above Step 3-2 is actually to solve the solution of the following minimization problem:

[0023] In the above formula ①, the purpose of minimizing is to achieve slow feature extraction, because the column vectors in the difference matrix reflect the temporal variation of the column vectors in the atomic block data matrix X b ; in addition, adding a regularization term on the basis of minimizing is to further take into account the neighborhood relationship of the column vectors in X b ; that is, after the transformation of v, the corresponding feature components should not only change slowly, but also ensure that the data points close in the original space distribution are also close in the feature space.

[0024] In summary, the method of the present invention analyzes and extracts the temporal and spatial neighborhood distribution features in each sub-block data matrix through step 3-2. For the solution of the above formula ①, the Lagrange multiplier method can be used to achieve it, and then it is transformed into two processes in step 3-2, that is: first perform whitening processing, and then perform singular value decomposition.

[0025] In addition, in steps 3-3 to 3-4, the method of the present invention selects to combine the obtained B decomposition vectors w1, w2,..., w B After that, through the participation of the global data matrix X0, the corresponding global load vector is calculated, and then it is divided into the load vectors corresponding to each sub-block data matrix. It can be seen that when determining the decomposition vectors of multiple sub-block data matrices, the present invention adheres to the mutually staggered solution idea from local sub-blocks to the global and then to local sub-blocks, not only considering the independence of each variable sub-block, but also considering its cooperation as a complete whole. Therefore, it can comprehensively represent the temporal and spatial distribution features of the learning data change from two levels of global variables and local sub-block variables.

[0026] Considering that the implementation of distributed fault detection first requires dividing the M variables of the industrial process into B variable sub-blocks, the division method directly restricts - which sub-block data matrices corresponding to the variable sub-blocks will be subjected to feature analysis subsequently. Although the most direct way is to divide the variables belonging to each production unit of the industrial process into a variable sub-block, this method requires the participation of human subjective factors. Especially, there are feedback loops between different production units of the industrial process, resulting in ambiguous division and attribution of some variables.

[0027] As an improvement of the method of the present invention, step 2 divides the M variables into non-overlapping B variable sub-blocks according to steps 2-1 to 2-3 shown below.

[0028] Step 2-1: According to the formula C = XX TAfter calculating matrix C and then calculating the eigenvector β corresponding to the maximum eigenvalue of C, update β through ; where |β| represents taking the absolute value of all elements in β.

[0029] Step 2-2: Denote the average value of all elements in β as β0, and then select the corresponding variables from the M variables according to the positions of the elements in β that are greater than β0 to form a variable sub-block.

[0030] Step 2-3: Delete the row vectors at the same positions in X according to the positions of the elements in β that are greater than β0, and then return to Step 2-1 to continue dividing to obtain the next variable sub-block.

[0031] It should be noted that calculating the eigenvector corresponding to the maximum eigenvalue of C = XX T in the above Step 2-1 is actually to achieve maxβ T XX T β, that is: achieving β T XX T β to maximize, so as to know which variables in X have the greatest correlation. Therefore, this improvement scheme realizes the division of variable sub-blocks by the way of the greatest correlation between variables, that is, dividing the variables with greater correlation into one variable sub-block, which is an optimal scheme for dividing variable sub-blocks that does not depend on process knowledge and does not require subjective human participation.

[0032] In addition, in Step 3-1, determining the corresponding near-neighborhoods for the column vectors in X b to obtain the corresponding near-neighbor reconstruction matrix is for subsequent feature extraction of the spatial distribution relationship of the column vectors in the sub-block data matrix. Therefore, how to determine the near-neighborhood is directly related to the representation learning of the characteristics of the data spatial distribution relationship.

[0033] To better describe the spatio-temporal distribution relationship in each sub-block data matrix, an improvement scheme of the method of the present invention is to determine the near-neighborhoods of each column vector according to the following Steps 3-1.A to 3-1.B.

[0034] Step 3-1.A: Calculate the distance between any two column vectors in X b , and then determine the c column vectors with the smallest distance to the b i-th column vector in X as the corresponding near-neighborhood; where c is an integer not less than 5 and not greater than 10.

[0035] Step 3-1.B: Find the d column vectors with the sampling time closest to b from the column vectors of X , and supplement these d column vectors into in the corresponding neighborhood; where d is the smallest even number not less than c.

[0036] In the technical field, determining the neighborhood through step 3-1.A is the most common method, that is, the c-order neighborhood rule well-known to those skilled in the art. However, this conventional method of determining the neighborhood only describes the spatial distribution relationship between data points from the perspective of spatial distance. In contrast, on the basis of determining the neighborhood by distance, the above step 3-1.B further supplements the column vectors close in sampling time into the neighborhood, so as to describe the spatio-temporal distribution relationship of data points from both spatial and temporal perspectives, which is beneficial to subsequent feature analysis and extraction.

[0037] As a further improvement of the method of the present invention, it is preferably to determine the neighborhood reconstruction matrix Φ in step 3-1 according to steps B1 to B3 as follows b .

[0038] Step B1: Initialize Φ b to a zero matrix, and merge the column vectors in the corresponding neighborhood into the input matrix and record the position numbers of the column vectors in the corresponding neighborhood in X b .

[0039] Step B2: After calculating the neighborhood reconstruction coefficient vector through , then update where , denotes the sum of all elements in the calculation of .

[0040] Step B3: According to the position numbers recorded in step B1, assign the elements in to the elements in the same positions of the i-th column vector of Φ b .

[0041] The traditional method of determining the neighborhood reconstruction matrix according to the neighborhood is to directly calculate the average value or the distance-weighted average value of all column vectors in , while the above step B2 calculates the reconstruction coefficient vector through the least squares algorithm, which can describe the relationship between and its neighborhood from the perspective of minimizing the reconstruction error.

[0042] As a further improvement of the method of the present invention, step 3-1 preferably calculates the difference matrix in the following manner

[0043] First, let X bThe column vectors from the 1st column to the second last column in form the input matrix Z1, and X b The column vectors from the 2nd column to the last column in form the output matrix Z2.

[0044] Secondly, calculate the time series regression matrix C according to the formula ; finally, calculate the difference matrix of X through b For those skilled in the art, the first-order difference is the most commonly used way to calculate the difference, that is, the difference matrix can be directly obtained according to Z2 - Z1

[0045] In contrast, the above improved technical solution uses the time series regression matrix to calculate the difference matrix, and can take into account the cross-correlation between variables, because the calculation of the time series regression matrix C satisfies the requirement of minimizing the sum of all elements in Z2 - CZ1. Physically speaking, the first-order difference reflects the speed of change, and the second-order difference reflects the acceleration of change. Therefore, the difference matrix calculated using the second-order difference can better reflect the essence of the change, and minimizing the second-order difference can analyze and extract the time series features with slow changes at a deeper level.

[0046] As a further improvement of the method of the present invention, when calculating the difference matrix for X

[0047] b it is preferred to calculate the difference matrix according to the second-order difference, that is: calculate the column vectors in according to ; where is obtained Among them, represents the difference vector corresponding to x0 in , x -1 , x0, x +1 represents any three column vectors with continuous sampling times in X b .

[0048] As a further improvement of the method of the present invention, the present invention provides a simple and direct technical solution for determining the regularization coefficient η, specifically: First, record the sum of the diagonal elements in as α1, and record the sum of the diagonal elements in as α2; secondly, set η = α2 / (α1 + α2).

[0049] Since the sum of the diagonal elements in the matrix is equal to the sum of its singular values, reflecting the energy information of the data change in the matrix, setting the regularization coefficient η in this way is very objective and direct, and is easy for engineers to implement.

[0050] As a further improvement of the method of the present invention, it is preferably to perform the whitening process in step 3-2 according to the following steps C1 to C3.

[0051] Step C1: Set Y = X b After that, according to the formula Y τ = [y τ , y τ , …, y τ+n-d T Form d sub-block time series matrices Y1, Y2, …, Y d ; where n is equal to the total number of column vectors in Y, and y τ , y τ , …, y τ+n-d respectively represent the column vectors of the τ-th column, the τ+1-th column, to the τ+n-d-th column in Y, and the time series number τ = 1, 2, …, d

[0052] Step C2: Set the whitening vector equal to any non-zero vector of M b ×1 dimension, and then repeat the following steps C2-1 to C2-2 until converges; where M b is equal to the total number of row vectors in X b .

[0053] Step C2-1: Calculate the matrix G1 according to the formula , and then solve the eigenvector α corresponding to the largest eigenvalue of G1.

[0054] Step C2-2: Calculate the matrix G2 according to the formula , and then solve the eigenvector ξ corresponding to the largest eigenvalue of G2, and then update through where α τ represents the τ-th element in α.

[0055] Step C3: Calculate the transformation vector g through , and then update Y according to , and then form d sub-block time series matrices Y1, Y2, …, Y d again, and then return to step C2 until M b whitening vectors and M b transformation vectors are obtained.

[0056] Step C4: Combine the obtained M b transformation vectors into a transformation matrix G, and combine the obtained M b whitening vectors into a whitening matrix D b After that, then through D b = D b (G T D​b ) -1 Update the whitening matrix D b 。

[0057] The purpose of whitening is to perform orthogonal decomposition on X b That is, through The score matrix S obtained by calculation b The row vectors in are orthogonal to each other. The most direct way to perform traditional whitening is through singular value decomposition to obtain orthogonal vectors. However, this technical solution obtains the corresponding whitening matrix through the eigenanalysis of the sub-block time series matrix, which can take into account the time series change characteristics of the data during whitening, so as to achieve more comprehensive feature analysis and extraction.

[0058] The present invention also provides a data-driven distributed fault detection device, including a data acquisition module, a central control module, a display module, and a fault alarm; the functions of each component are as follows.

[0059] The data acquisition module can obtain a set of data of the industrial process at the latest sampling moment t according to the sampling time interval and immediately send it to the central control module.

[0060] After receiving a set of measurement data, the central control module first performs standardization processing on it to obtain the column vector x t , and then calculates its corresponding fault detection index, and sends the fault detection index to the display module, and at the same time judges whether the fault detection index is greater than the threshold If so, trigger the fault alarm; if not, keep the fault alarm silent.

[0061] After receiving the fault detection index, the display module draws a monitoring curve graph with the latest sampling moment t as the abscissa and the fault detection index and the threshold as the ordinate.

[0062] The core key of this distributed fault detection device is that: the central control module stores a program for executing steps 3-1 to 3-5, so that after receiving the initialization instruction, it executes steps 1 to 4 to obtain the sub-block decomposition matrices W1, W2,... W B and the sub-block loading matrices P1, P2,... P B , as well as the threshold In addition, the central control module calculates the fault detection index according to steps A1 to A2.

[0063] Therefore, both the fault detection device disclosed in the present invention and the previously disclosed fault detection method perform eigenanalysis on the B sub-block data matrices using the manifold regularization slow feature analysis algorithm, that is, based on the same inventive concept.

[0064] In addition, the present invention also provides an improved distributed fault detection device, and a program for implementing Steps 2-1 to 2-3 is further stored in the central control module of the device, so as to divide M variables in the industrial process into B non-overlapping variable sub-blocks.

[0065] This improved distributed fault detection device can divide variables without the need for process knowledge or human participation, making it easier to be applied in practice. Description of the Drawings

[0066] Figure 1 is a schematic diagram of a fluid catalytic cracking process.

[0067] Figure 2 is a schematic flow diagram of the implementation feature analysis related to the present invention.

[0068] Figure 3 The fault detection curve graph corresponding to the first embodiment.

[0069] Figure 4 The fault detection curve graph corresponding to the second embodiment.

[0070] Figure 5 is an implementation flow chart of the whitening process improvement scheme.

[0071] Figure 6 is the fault detection curve graph corresponding to the eighth embodiment.

[0072] Figure 7 is a schematic composition structure diagram of the distributed fault detection device related to the present invention. Detailed Description of the Invention

[0073] The present invention will be described in detail below with reference to the drawings and specific embodiments.

[0074] In order to clearly and specifically describe the specific embodiments of the fault detection method related to the present invention, it is applied to the fault detection task of the fluid catalytic cracking process as shown in Figure 1 in. Figure 1 The process shown in consists of three main parts: a reactor, a regenerator, and a fractionator, which have important application values in the oil refining industry. The preheated heavy oil feed is sent to the bottom of the riser and then contacts the hot catalyst from the regenerator. The hydrocarbon vapor undergoes a catalytic cracking reaction during the process of passing through the riser, and then the cracked product vapor is sent to the fractionator. In the fractionator, the product vapor is separated into lighter fractions according to their different boiling points, such as diesel and gasoline. At each sampling moment, data of a total of M = 92 variables including temperature, pressure, and flow rate can be collected.

[0075] ForFigure 1 The process object shown applies a data-driven distributed fault detection method disclosed by the present invention, including the following steps 1 to 5 shown below.

[0076] Step 1: Obtain a data set when the industrial process is operating normally, and perform standardization processing on it to obtain a training data matrix X; among them, each column vector in X is composed of data corresponding to M = 92 variables.

[0077] Step 2: According to the Figure 1 composition structure of the fluid catalytic cracking process shown in, divide these M = 92 variables into B = 5 non-overlapping variable sub-blocks, and then correspondingly divide X into corresponding B = 5 sub-block data matrices X1, X2,..., X B .

[0078] Step 3: Use the manifold regularization slow feature analysis algorithm to perform feature analysis on X1, X2,..., X B to obtain corresponding sub-block decomposition matrices W1, W2,..., W B and sub-block load matrices P1, P2,..., P B ; the corresponding implementation process is as Figure 2 shown, including the following steps 3-1 to 3-5 shown below.

[0079] Step 3-1: Determine the corresponding neighborhood for each column vector in X b to obtain the corresponding neighborhood reconstruction matrix Φ b of X b , and then calculate the difference matrix b of X where X b represents the b-th sub-block data matrix, and the number b = 1, 2,..., 5

[0080] Step 3-2: Perform whitening processing on X b to obtain the whitening matrix D b , and then calculate the symmetric matrix G through the formula . Perform singular value decomposition on G, and denote the decomposition vector corresponding to the smallest singular value as w b ; among them, the regularization coefficient η = 1.

[0081] Step 3-3: After merging the B decomposition vectors w1, w2,..., w B into a global decomposition vector w, merge the B sub-block data matrices X1, X2,..., X B into a global data matrix X0, and then calculate the global load vector p through .

[0082] Step 3-4: Divide p into B load vectors p1, p2, …, p correspondingly B , and then through After updating each sub-block data matrix, return to Step 3-2 until each sub-block data matrix corresponds to K = 10 decomposition vectors and K = 10 load vectors.

[0083] Step 3-5: Combine the K decomposition vectors and K load vectors corresponding to each sub-block data matrix into sub-block decomposition matrices W1, W2, …, W B and sub-block load matrices P1, P2, …, P B .

[0084] Step 4: Calculate the corresponding fault detection index for each column vector in X, and then determine the upper limit of the change of the fault detection index Among them, the calculation process of the fault detection index includes the following Steps A1 to A2, and the upper limit of the change is determined using the kernel density estimation method.

[0085] Step A1: Divide any given column vector x into corresponding B sub-block column vectors x1, x2, …, x according to the B variable sub-blocks in Step 2 B , and then through and Calculate the sub-block eigenvector s b and sub-block error vector ε b in sequence.

[0086] Step A2: Calculate the b-th sub-block detection index through , and then set the fault detection index corresponding to x Among them, represents the upper limit of the change of the b-th sub-block detection index.

[0087] It should be noted that after having the sub-block detection indexes corresponding to each column vector in X, the maximum value of them can be directly determined as the corresponding upper limit of the change, which is the simplest and most direct way to determine the upper limit of the change; of course, the average value of the largest multiple of the fault detection indexes can also be taken as the upper limit of the change. In this embodiment, using the kernel density estimation method to determine the upper limit of the change is only a specific description of the implementation manner involved in the present invention, rather than a limitation on the technical solution of the present invention.

[0088] Step 5: Obtain a set of data of the industrial process at the latest sampling moment t, and perform standardization processing on it to obtain the column vector x t , and then calculate its corresponding fault detection index and determine whether the industrial process has a fault at the latest sampling moment t according to whether it exceeds .

[0089] Another set of fault detection metrics corresponding to the data at 1000 latest sampling moments is plotted in Figure 3 Since the upper change limit in Step 4 remains unchanged once determined, therefore Figure 3 the horizontal line in

[0090] From Figure 3 it can be found that starting from the 400th latest sampling moment, the fault detection metrics at the corresponding sampling moments are located above the horizontal line in Figure 3 for the vast majority of the time, that is, the fault detection metrics are greater than the threshold, and a fault alarm can be triggered.

[0091] In the second embodiment of the present invention, when dividing variable sub - blocks, it no longer depends on the composition structure of the process, but divides M = 92 variables into B = 9 non - overlapping variable sub - blocks according to the following Steps 2 - 1 to 2 - 3.

[0092] Step 2 - 1: Calculate matrix C according to the formula C = XX T After calculating the maximum eigenvalue of C and its corresponding eigenvector β, update β through where |β| represents taking the absolute value of all elements in β.

[0093] Step 2 - 2: Denote the average value of all elements in β as β0, and then select the corresponding variables from the M variables according to the positions of the elements in β that are greater than β0 to form a variable sub - block.

[0094] Step 2 - 3: Delete the row vectors at the same positions in X according to the positions of the elements in β that are greater than β0, and then return to Step 2 - 1 to continue dividing to obtain the next variable sub - block.

[0095] Based on this improved variable sub - block division implementation plan, fault detection is performed on the process data at the same 1000 latest sampling moments, and the obtained fault detection curve graph is as shown in Figure 4 From Figure 4 and Figure 3 it can be found that after the fault occurs, that is, starting from the 401st sampling moment, Figure 4 the fault detection curve shown in

[0096] The third embodiment of the present invention is based on the above - mentioned first or second embodiment, and determines the neighborhood in Step 3 - 1 according to the following Steps 3 - 1.A to 3 - 1.B.

[0097] Step 3 - 1.A: Calculate the distance between any two column vectors in X b and then compare with Xb the $i$-th column vector in The $c$ column vectors with the smallest distances among them are determined as the corresponding neighborhood; where $c = 5$, and the number $i = 1, 2, \ldots, N$ b , and $N$ b is equal to the total number of column vectors in $X$ b .

[0098] Step 3-1.B. From the column vectors of $X$ b , find the $d$ column vectors whose sampling times are closest to, and supplement these $d$ column vectors into the corresponding neighborhood; where $d = 6$.

[0099] The fourth embodiment of the present invention is based on any one of the above first to third embodiments, and determines the neighborhood reconstruction matrix $\varPhi$ in Step 3-1 according to the following Steps B1 to B3 b .

[0100] Step B1. Initialize $\varPhi$ b as an $N$ b × $N$ b dimensional zero matrix, merge the column vectors in the corresponding neighborhood into the input matrix and record the position numbers of the column vectors in the corresponding neighborhood in $X$ b ; where $N$ b is equal to the total number of column vectors in $X$ b , denotes the $i$-th column vector in $X$ b , and the number $i = 1, 2, \ldots, N$ b .

[0101] Step B2. After calculating the neighborhood reconstruction coefficient vector through, then update by the formula where , denotes the sum of all elements in the calculation of .

[0102] Step B3. According to the position numbers recorded in Step B1, assign the elements in to the elements in the same positions of the $i$-th column vector of $\varPhi$ b .

[0103] The fifth and sixth embodiments of the present invention respectively relate to improvements in calculating the difference matrix; among them, the fifth embodiment calculates the difference matrix in the following manner

[0104] First, form the input matrix Z1 from the column vectors of the first column to the penultimate column in X b , and form the output matrix Z2 from the column vectors of the second column to the last column in X b . Secondly, calculate the time series regression matrix C according to the formula . Finally, calculate the difference matrix of X through b .

[0105] The sixth embodiment calculates the difference matrix according to the calculation method of the second-order difference . Specifically: calculate to obtain each difference vector in ; among them, represents the difference vector corresponding to x0 in -1 , x +1 , x0, x b represent any three column vectors with continuous sampling times in X

[0106] In the above embodiments, the regularization coefficient is set to η = 1. In the seventh embodiment of the present invention, on the basis of the above six embodiments, the regularization coefficient is particularly set to η = α2 / (α1 + α2); where α1 is the sum of the elements on the diagonal of , and α2 is the sum of the elements on the diagonal of .

[0107] The eighth embodiment of the present invention relates to an improvement in the whitening process, and its implementation flowchart is as shown in Figure 5 , and the corresponding whitening process includes the following steps C1 to C4

[0108] Step C1: After setting Y = X b , form d sub-block time series matrices Y1, Y2,..., Y τ = [y τ , y τ ,..., y τ+n-d T according to the formula d .

[0109] Step C2: After initializing to be equal to any non-zero vector of dimension M b ×1, repeat the following steps C2-1 to C2-2 until converges

[0110] Step C2-1: Calculate the matrix G1 according to the formula , and then solve the eigenvector α corresponding to the maximum eigenvalue of G1​

[0111] Step C2-2: Calculate matrix G2 according to the formula , then solve the eigenvector ξ corresponding to the maximum eigenvalue of G2, and update where α represents the τ-th element in α. τ

[0112] Step C3: Calculate the transformation vector g through , then update Y according to , and form d sub-block time series matrices Y1, Y2, …, Y d , and then return to Step C2 until obtaining M b whitening vectors and M b transformation vectors.

[0113] Step C4: Combine the obtained M b transformation vectors into a transformation matrix G, and combine the obtained M b whitening vectors into a whitening matrix D b , and then update the whitening matrix D through D b = D b (G T D b ) -1 . b

[0114] All of the above seven embodiments can perform whitening processing on X according to the whitening processing implementation process in Figure 5 to obtain the whitening matrix D b . Taking the second embodiment as an example and combining it with the whitening processing implementation process in b Figure 5 , the corresponding fault detection curve is as shown in Figure 6 .

[0115] It can be found from Figure 6 that compared with Figure 3 and Figure 4 , using the above Steps C1 to C4 to perform whitening processing achieves a more excellent fault detection effect.

[0116] The present invention also provides a data-driven distributed fault detection device, and its structural composition schematic diagram is as shown in Figure 7 , including a data acquisition module, a central control module, a display module, and a fault alarm.

[0117] Figure 7The data acquisition module therein can obtain a set of data of the industrial process at the latest sampling moment t according to the sampling time interval and immediately send it to the central control module; after receiving a set of measurement data, the central control module first performs standardization processing on it to obtain the column vector x t , then calculates its corresponding fault detection index, and sends the fault detection index to the display module, and at the same time judges whether the fault detection index is greater than the threshold If so, trigger the fault alarm; if not, keep the fault alarm silent; after receiving the fault detection index, the display module draws a monitoring curve graph with the latest sampling moment t as the abscissa and the fault detection index and the threshold as the ordinate.

[0118] Figure 7 The program for implementing steps 3-1 to step 3-5 is stored in the central control module shown in, so that after receiving the initialization instruction, steps 1 to 4 are implemented to obtain the sub-block decomposition matrices W1, W2,...W for calculating the fault detection index B and the sub-block load matrices P1, P2,...P B , as well as the threshold Moreover, the central control module calculates the fault detection index according to steps A1 to A2.

[0119] Finally, as an improvement of this fault detection device, Figure 7 The program for implementing steps 2-1 to step 2-3 is also stored in the central control module in, so as to divide the M variables into B variable sub-blocks that do not overlap with each other.

[0120] The explanatory description of the technical solution disclosed in the present invention in combination with the accompanying drawings is illustrative rather than restrictive. The technical solution disclosed in the present invention is not limited to the above-mentioned embodiments. Modifications made without departing from the purpose of the technical solution disclosed in the present invention and within the scope of protection of the claims all fall within the protection scope of the present invention.

Claims

1. A data-driven distributed fault detection method, comprising the following steps 1 to step 5: Step 1, obtain a data set when the industrial process is operating normally, and perform standardization processing on it to obtain a training data matrix X; Among them, Each column vector of X is composed of data corresponding to M variables; Step 2: Divide the M variables into B non-overlapping variable sub-blocks, and correspondingly divide X into B sub-block data matrices X1, X2, …, X B ; Step 3: Use the manifold regularization slow feature analysis algorithm to perform feature analysis on X1, X2, …, X B to obtain the corresponding sub-block decomposition matrices W1, W2, …, W B and sub-block loading matrices P1, P2, …, P B ; Step 4: Calculate the corresponding fault detection indexes for each column vector in X, and then determine the upper limit of the change of the fault detection indexes Step 5: Obtain a set of data of the industrial process at the latest sampling moment t, and perform standardization processing on it to obtain a column vector x t After that, calculate its corresponding fault detection index Then, based on whether it exceeds to determine whether a fault occurs in the operation of the industrial process; Its feature is that the implementation process of step 3 includes the following steps 3-1 to step 3-5; Step 3-1: For X b Determine the corresponding neighborhood for each column vector in b to obtain the neighborhood reconstruction matrix Φ of X b Then, calculate the difference matrix of X b wherein, X represents the b-th sub-block data matrix, and the number b = 1, 2, …, B; b ​ Step 3-2: Perform whitening processing on X b to obtain the whitening matrix D b After that, through the formula calculate the symmetric matrix G, perform singular value decomposition on G, and denote the decomposition vector corresponding to the smallest singular value as w b ; where I represents the identity matrix, the superscript T represents the transpose of the matrix or vector, and the value range of the regularization coefficient η is 0 < η ≤ 1; Step 3-3: Combine B decomposed vectors w1, w2, …, w B into a global decomposed vector w, and combine B sub-block data matrices X1, X2, …, X B into a global data matrix X0, and then calculate the global load vector p through the following operations; Step 3-4: Correspondingly divide p into B load vectors p1, p2, …, p B , and then through After updating each sub-block data matrix, return to Step 3-2 until each sub-block data matrix corresponds to K decomposition vectors and K load vectors; where, w b and p b respectively represent the decomposition vector and the load vector corresponding to the b-th sub-block data matrix; Step 3-5: Combine the K decomposition vectors and K loading vectors corresponding to each sub-block data matrix into sub-block decomposition matrices W1, W2, …, W B and sub-block loading matrices P1, P2, …, P B ; The calculation process of the fault detection index includes the following steps A1 to step A2: Step A1. According to the B variable sub-blocks in Step 2, any given column vector x is divided into corresponding B sub-block column vectors x1, x2, …, x B and then, through and successively calculate the sub-block eigenvector s b and the sub-block error vector ε b ; Step A2, by calculating the detection index of the b-th sub-block and then setting the fault detection index corresponding to x where represents the upper limit of the change in the detection index of the b-th sub-block.

2. The distributed fault detection method according to claim 1, wherein Divide the M variables into B non-overlapping variable sub-blocks according to the following steps 2-1 to step 2-3: Step 2-1: Calculate matrix C according to the formula C = XX T After calculating matrix C and then calculating the eigenvector β corresponding to the maximum eigenvalue of C, update β through where |β| represents taking the absolute value of all elements in β; Step 2-2, record the average value of all elements in β as β0, and then select the corresponding variables from the M variables according to the positions of the elements in β that are greater than β0 to form a variable sub-block; Step 2-3, delete the row vectors at the same positions in X according to the positions of the elements in β that are greater than β0, and then return to step 2-1 to continue dividing to obtain the next variable sub-block.

3. A distributed fault detection method according to claim 1 or 2, characterized in that The determination process of the neighborhood in step 3-1 includes the following steps 3-1.A to step 3-1.B: Step 3-1.A. Calculate X b the distance between any two column vectors in, and then the c column vectors with the smallest distance from the b i-th column vector in are determined as the corresponding neighborhood; where c is an integer not less than 5 and not greater than 10, and the number i = 1, 2,..., N b , and N b is equal to the total number of column vectors in X b ; Step 3-1.B. From the column vectors of X b , find d column vectors whose sampling times are closest to , and supplement these d column vectors into the corresponding neighborhood; where d is the smallest even number not less than c.

4. A distributed fault detection method according to claim 1 or 2, characterized in that, Determine the nearest neighbor reconstruction matrix Φ in step 3-1 according to steps B1 to B3 shown below b : Step B1. Initialize Φ b as an N b ×N b zero matrix, and merge the column vectors in the corresponding neighborhood of into the input matrix and record the position numbers of the column vectors in the corresponding neighborhood of in X b ; where N b is equal to the total number of column vectors in X b , represents the i-th column vector in X b , and the number i = 1, 2, …, N b ; Step B2: After calculating the nearest neighbor reconstruction coefficient vector , update it through the formula where denotes the sum of all elements in the calculation of ; ​ Step B3. According to the position numbers recorded in Step B1, assign the elements in to the elements at the same positions in the i-th column vector of Φ b .

5. A distributed fault detection method according to claim 1 or 2, characterized in that, In step 3-1, the differential matrix is calculated in the following way: First, the column vectors from the first column to the second last column in X b form the input matrix Z1, and the column vectors from the second column to the last column in X b form the output matrix Z2; Second, the time series regression matrix C is calculated according to the formula ; Finally, the differential matrix of X b is calculated through 6. A distributed fault detection method according to claim 1 or 2, characterized in that In step 3-1, the differential matrix is calculated in the following way: According to calculate to obtain each differential vector in; where represents the differential vector corresponding to x0 in, x -1 , x0, x +1 represents any three column vectors with continuous sampling times in X b .

7. A distributed fault detection method according to claim 1 or 2, characterized in that The determination method of the regular coefficient η is as follows: First, denote the sum of the elements on the diagonal in as α1, and denote the sum of the elements on the diagonal in as α2; Second, set η = α2 / (α1 + α2).

8. A distributed fault detection method according to claim 1 or 2, characterized in that, In step 3-2, for X b The process of performing whitening processing includes steps C1 to C4 shown below: Step C1: Set Y = X b After that, according to the formula Y τ = [y τ , y τ+1 , …, y τ+n-d T Construct d sub-block timing matrices Y1, Y2, …, Y d ; where n is equal to the total number of column vectors in Y, and y τ , y τ+1 , …, y τ+n-d represent the column vectors of the τ-th column, the (τ + 1)-th column, up to the (τ + n - d)-th column in Y respectively, and the timing number τ = 1, 2, …, d;​ Step C2: Set the whitening vector θ to be equal to any non - zero vector of dimension M×1, and then repeatedly execute the following steps C2 - 1 to C2 - 2 until θ converges; where M b is equal to the total number of row vectors in X b ×1, and then repeat the following steps C2 - 1 to C2 - 2 until θ converges; where M b equals the total number of row vectors of X Step C2-1: After calculating the matrix G1 according to the formula G1 = [Y1θ, Y2θ, …, Y d θ] T [Y1θ, Y2θ, …, Y d θ], solve the eigenvector α corresponding to the maximum eigenvalue of G1; Step C2-2: According to the formula calculate matrix G2. After solving for the eigenvector ξ corresponding to the maximum eigenvalue of G2, update θ through ; where α τ represents the τ-th element in α. Step C3: Calculate the transformation vector g through g = YY T θ, and then update Y according to Y = Y - gθ T After updating Y, reorganize the d sub-block time series matrices Y1, Y2, …, Y d again, and then return to Step C2 until M b whitening vectors and M b transformation vectors are obtained; Step C4. Merge the obtained M b transform vectors into a transformation matrix G, and merge the obtained M b whitening vectors into a whitening matrix D b . Then, update the whitening matrix D b through D b = D T (G b D -1 ) b .

9. A data-driven distributed fault detection device, comprising a data acquisition module, a central control module, a display module, and a fault alarm; wherein, The data acquisition module can obtain a set of data of the industrial process at the latest sampling moment t in real time according to the sampling time interval, and immediately send it to the central control module; After the central control module receives a set of measurement data, it first performs standardization processing on it to obtain the column vector x t , and then calculates its corresponding fault detection index. After that, it sends the fault detection index to the display module and simultaneously determines whether the fault detection index is greater than the threshold If so, trigger the fault alarm; if not, keep the fault alarm silent; After receiving the fault detection index, the display module plots a monitoring curve graph with the latest sampling moment t as the abscissa and the fault detection index and the change upper limit as the ordinate; It is characterized in that a program for implementing steps 3-1 to 3-5 is stored in the central control module, so that after receiving an initialization instruction, steps 1 to 4 are implemented to obtain sub-block decomposition matrices W1, W2, … W for calculating fault detection metrics B and sub-block load matrices P1, P2, … P B , as well as the upper limit of change The central control module calculates the fault detection metric according to steps A1 to A2.

10. A distributed fault detection device according to claim 9, wherein, The central control module also stores a program for implementing steps 2-1 to step 2-3, so as to divide the M variables of the industrial process into B non-overlapping variable sub-blocks.

Images