An ammonia synthesis first-stage furnace operation anomaly detection method based on time sequence difference feature extraction

Abstract

The application discloses a kind of based on timing difference feature extraction synthetic ammonia primary furnace operation abnormality detection method, to how adaptive is for synthetic ammonia primary furnace process data real-time extraction timing difference feature, to realize and complete high sensitivity synthetic ammonia primary furnace operation abnormality detection task.Specifically, the method of the present application adaptively establishes timing error generation model by the process data of synthetic ammonia primary furnace.On this basis, using the error of online autoregressive model implements synthetic ammonia primary furnace operation abnormality detection.The biggest technical advantage of the method of the present application is that: the present application adopts a kind of new timing difference feature analysis algorithm technology to carry out feature analysis extraction to the process data of synthetic ammonia primary furnace, can be adaptively updated and obtained from the timing difference feature of maximum timing error.Therefore, from this point of view, the method of the present application can use high sensitivity timing difference feature to implement anomaly detection in theory.

Description

TECHNICAL FIELD

[0001] The application relates to a chemical process anomaly detection method, in particular to a synthetic ammonia primary reformer operation anomaly detection method based on time sequence difference feature extraction. BACKGROUND

[0002] Ammonia is one of important inorganic chemical products and occupies an important position in the national economy. In the synthetic ammonia production process, coke, coal, coke oven gas, natural gas, naphtha, heavy oil and raw material gas (semi-water coal gas) and the like are compressed to a certain pressure after purification and refining (main components are hydrogen and nitrogen) and then enter an ammonia synthesis tower to generate ammonia and other by-products under the action of high temperature and a catalyst. The primary reformer is a key equipment of a synthetic ammonia device of a chemical fertilizer plant, and the synthetic ammonia primary reformer is used for converting natural gas and steam under the action of a catalyst to prepare hydrogen required in the synthetic ammonia production. If the operation level of the synthetic ammonia primary reformer is unstable, the energy consumption of the synthetic ammonia device and the ammonia yield are directly affected.

[0003] It can be said that the operation state monitoring of the synthetic ammonia primary reformer is particularly important, which is not only a high-energy-consumption unit but also determines whether the hydrogen yield and purity can be effectively improved. Since the synthetic ammonia primary reformer is an important link in hydrogen production, according to process experience, it is necessary to maintain the stability of the hydrogen yield by ensuring the stability of the combustion in the furnace. However, the combustion condition of the primary reformer needs a high-priced mass spectrometer, and the synthetic ammonia device considering the design cost generally does not have a matching installation. Fortunately, the relatively low-priced distributed control system (DCS) is an essential software system of the synthetic ammonia complete equipment, which can measure process data such as temperature, pressure, flow, liquid level and the like in real time at a short interval. These easily obtained process data lay a solid data foundation for indirectly implementing data-driven synthetic ammonia primary reformer operation anomaly detection.

[0004] Although monitoring the flow of fuel and air waste can indirectly reflect the combustion condition in the primary reformer, the temperature change caused by the combustion process has time sequence correlation. Therefore, when implementing anomaly detection by using the easily obtained process data, it is necessary to consider the time sequence change feature between the sampling data. The traditional time sequence feature analysis and extraction technology that can be used for anomaly detection usually only focuses on the time sequence relationship feature mining under the normal, i.e. expected, operation state, and it is unknown whether the mined time sequence relationship feature is conducive to anomaly detection. Therefore, from the perspective of time sequence feature analysis, solving the synthetic ammonia primary reformer operation anomaly detection problem must also consider how to extract potential features reflecting time sequence differences. SUMMARY

[0005] The primary technical problem addressed by this invention is how to adaptively extract time-series difference features from the process data of ammonia synthesis primary furnaces in real time, thereby enabling and completing the task of detecting operational anomalies in ammonia synthesis primary furnaces. Specifically, the method adaptively establishes a time-series error generation model based on the process data of ammonia synthesis primary furnaces. Specifically, different online process data corresponds to different time-series error generation models. Based on this, the error of the online autoregressive model is used to detect operational anomalies in ammonia synthesis primary furnaces.

[0006] The technical solution adopted by the method of the present invention to solve the above-mentioned problem is: a method for detecting abnormal operation of a synthetic ammonia primary furnace based on time series difference feature extraction, comprising the following steps:

[0007] Step (1): Obtain N sets of process data x1, x2, ..., x2 of the first stage of the synthetic ammonia furnace under normal operation from the DCS historical database. N ; where x i ∈R 13×1 Represents the i-th group of process data, with subscript i∈{1, 2, ..., N}, R 13×1 represents a 13×1 dimensional real number vector, R represents a real number set, and each set of process data includes 13 data, specifically involving three types of data: flow, pressure, and temperature.

[0008] It is worth noting that the arrangement order of m=13 data in each set of process data is as follows: natural gas intake flow rate, exhaust gas intake flow rate, pressure and temperature at the heat exchanger outlet, natural gas temperature at the preheater outlet, gas pressure in the furnace at the outlet of the first furnace, gas temperature at the inlet of the first furnace, gas temperature at the upper left of the first furnace, gas temperature at the upper right of the first furnace, gas temperature in the mixed furnace at the top of the first furnace, conversion gas temperature at the left outlet of the first furnace, conversion gas temperature at the right outlet of the first furnace, and conversion gas temperature at the outlet of the first furnace.

[0009] Step (2): N sets of process data x1, x2, ..., x N Form a training data matrix X = [x1, x2, ..., x N ], the row vectors of each row in X are standardized to obtain the reference data matrix The standardization method is to subtract the mean of the row vector from the row vector and then divide it by the standard deviation of the row vector. 13×N Represents a 13×N-dimensional real matrix.

[0010] Step (3): After setting the time series correlation order equal to D, use the following formula And the following formula ① respectively form the timing matrix X D ∈R 13×(N-D+1) and the time series expansion matrix Z∈R13(D-1)×(N-D+1) :

[0011]

[0012] in, The reference data matrix is ​​represented in turn The column vectors of the 1st, 2nd, to Nth columns in , express The column vector of the D-2th column in , express The column vector of the D-1th column in , express The column vector of the N-2th column in express The column vector of the N-D+1th column in R 13(D-1)×(N-D+1) represents a real matrix of 13(D-1)×(N-D+1), R 13×(N-D+1) Represents a 13×(N-D+1)-dimensional real matrix.

[0013] Step (4): X D The column vectors from the 1st column to the N-D+1th column are taken as online data vectors in turn Implement time series difference feature analysis to calculate the corresponding anomaly detection indicators Q1, Q2, ..., Q N-D+1 Then determine the upper control limit Q of the anomaly detection index lim The specific implementation process is shown in steps (4.1) to (4.6).

[0014] Step (4.1): After initializing i=1, set Equal to X D The column vector of column i in , sets the online time series vector z t is equal to the column vector of the i-th column in Z.

[0015] Step (4.2): X D The remaining ND column vectors except the i-th column vector constitute the reference timing matrix The remaining ND column vectors of Z except the i-th column vector form the reference timing expansion matrix

[0016] Step (4.3): Perform time series difference feature analysis to obtain the online load vector w t and the online coefficient vector β t The specific implementation process is shown in steps (A) to (D).

[0017] Step (A): Initialize the feature vector g∈R m×1 Equal to an arbitrary m×1-dimensional real vector.

[0018] Step (B): according to the formula Calculate the online coefficient vector β t ; wherein, the superscript T represents the transpose of a matrix or a vector, the matrix L w and L D are respectively shown in formula ②:

[0019]

[0020] In the above formula, I D represents a (D-1) × (D-1) dimensional identity matrix, represents the calculation of the Kronecker product between I D and g, and the specific calculation result is as follows:

[0021]

[0022] Step (C): solve the generalized eigenvalue problem G = λL β , and then normalize g using formula ; wherein, the matrix G and the matrix L β are shown in formula ③:

[0023]

[0024] In the above formula, represents the calculation of the Kronecker product between β t and I M , and I M represents a 13 × 13 dimensional identity matrix.

[0025] Step (D): determine whether the characteristic vector g converges, and the convergence criterion is that the elements in g no longer change; if not, return to step (B); if yes, set the online load vector w t = g, and the timing difference feature analysis process is completed, thereby obtaining the online load vector w t and the online coefficient vector β t .

[0026] Step (4.4): calculate the anomaly detection index Q i according to formula ; wherein, online timing error

[0027] Step (4.5): determine whether i is less than N-D+1; if yes, set i = i+1, and then set equal to the column vector of the i-th column in X D , and set z tIf it is equal to the column vector of the i-th column in Z, return to step (4.2); if not, obtain N-D+1 anomaly detection indicators Q1, Q2, ..., Q N-D+1 Then, proceed to step (4.6).

[0028] Step (4.6): For Q1, Q2, ..., Q N-D+1 Kernel Density Estimation (KDE) is performed to obtain the confidence limit Q at a confidence level of α = 99.5%. KDE , and Q1, Q2, ..., Q N-D+1 The average value of the largest anomaly detection index ξ is recorded as Q AVG Then determine the upper control limit Q lim =max{Q KDE , Q AVG}; where ξ represents the largest integer not greater than n×99%, n=N-D+1, max{Q KDE , Q AVG} means taking Q KDE and Q AVG The maximum value in .

[0029] It should be noted that the above steps (A) to (D) are actually a new time series difference feature analysis algorithm technology designed by the method of the present invention. The goal of this algorithm technology is to provide real-time data vector Find the online load vector w t and the online coefficient vector β t , maximize the online timing error While making the reference timing matrix and reference timing expansion matrix The timing error between them is minimized, that is:

[0030]

[0031] In the above formula, the numerator is the square of the online timing error, while the denominator is the square of the length of the timing error vector between the reference timing matrix and the reference timing expansion matrix. The solution to formula (4) can be achieved using the classic Lagrange multiplier method, that is, by constructing the Lagrangian function J as shown below:

[0032]

[0033] Calculate J relative to w respectively t and β t After setting the partial derivative equal to 0, we can get: Gw t =λL β w t And the following equation relationship:

[0034]

[0035] i.e.L w w t = L D β t . Thus, solving w t is equivalent to solving a generalized eigenvalue problem: Gw t = λL β w t , and solving β t is equivalent to a least squares regression problem: When w t converges, β t can be uniquely determined. Therefore, to avoid the problem of using different technical terms to describe the same symbol, the generalized eigenvalue problem in the aforementioned step (C) is equivalently expressed as G = λL β g.

[0036] In addition, it should be noted that the rank of matrix G is equal to 1, so G = λL β g actually has only one eigenvalue and corresponds to one eigenvector. The eigenvalue is the largest eigenvalue. However, matrix L β is not invertible, so solving the generalized eigenvalue in the aforementioned step (C) also needs to be implemented through singular value decomposition, and the specific implementation process is shown in steps (one) to (four).

[0037] Step (one): singular value decomposition L β is implemented, and L β = UΛU T is obtained, only the non-zero singular values greater than 10 -4 are retained, thereby obtaining the unitary matrix U and the singular value diagonal matrix Λ, and then the matrix

[0038] Step (two): after the matrix G β is calculated according to the formula , the eigenvector p is initialized to be an arbitrary 13 × 1-dimensional real vector.

[0039] Step (three): the eigenvector p is updated according to the formula p = G β p and in turn.

[0040] Step (four): it is judged whether the eigenvector p converges or not; if not, return to step (three); if yes, after the final eigenvector p is obtained, step (five) is executed.

[0041] Step (five): the eigenvector g is calculated according to the formula .

[0042] Step (5): Use DCS to obtain the process data x of the first stage of synthetic ammonia furnace at the latest sampling time t ∈R 13×1 Then, the same normalization process as in step (2) is performed to obtain the online data vector corresponding to the latest sampling time.

[0043] Step (6): The online data vector corresponding to the previous sampling moment Online data vectors corresponding to the first two sampling moments Online data vector corresponding to the previous D-1 sampling moments Merge to form an online time series vector

[0044] Step (7): Setup and Then, the online data vector Perform time series difference feature analysis, as shown in steps (A) to (D), to obtain the corresponding online load vector w t and the online coefficient vector β t , and then calculate the anomaly detection index at the latest sampling time

[0045] Step (8): Determine Q t Is it greater than Q lim If not, the first stage of the synthetic ammonia furnace is operating normally, and the process returns to step (5) to continue to use the process data at the latest sampling time to perform abnormality detection; if so, execute step (9) to decide whether to trigger an abnormality alarm.

[0046] Step (9): Return to step (5) and continue to use the process data at the latest sampling moment to perform anomaly detection. If the anomaly detection index of A consecutive sampling moments is greater than Q lim , an abnormal alarm is triggered; otherwise, the synthetic ammonia stage furnace operates normally; where A is equal to the quotient of the shortest time required to trigger an abnormal alarm divided by the sampling interval.

[0047] Through the above-mentioned implementation steps, the main advantages of the method of the present invention are introduced as follows.

[0048] The greatest technical advantage lies in the fact that this invention utilizes a novel time series difference feature analysis algorithm to analyze and extract features from the process data of the synthetic ammonia primary furnace, enabling adaptive updates to obtain time series difference features that maximize the time series error. Therefore, from this perspective, the method of the present invention can theoretically use highly sensitive time series difference features to implement anomaly detection. Compared to traditional methods that focus on analyzing and extracting time series features under normal operating conditions, traditional methods cannot guarantee the sensitivity of the extracted features for implementing anomaly detection. This means that the method of the present invention can better address the real-time detection of whether the synthetic ammonia primary furnace is operating normally, and can always guarantee the sensitivity of the extracted time series difference features to process data under abnormal operating conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 Schematic diagram of the implementation process of the method of the present invention. DETAILED DESCRIPTION

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

[0051] The present invention discloses a method for detecting abnormal operation of ammonia synthesis furnace based on time series difference feature extraction. Figure 1 The implementation flow diagram shown and sampling data of a domestic ammonia synthesis process are used to illustrate the specific implementation method of the method of the present invention.

[0052] Step (1): Obtain N sets of process data x1, x2, ..., x2 of the first stage of the synthetic ammonia furnace under normal operation from the DCS historical database. N .

[0053] Step (2): N sets of process data x1, x2, ..., x N The training data matrix X = [x1, x2, ..., x N ], the row vectors of each row in X are standardized to obtain the reference data matrix

[0054] Step (3): After setting the time series correlation order equal to D, use the following formula And the above formula ① respectively form the time series matrix X D ∈R 13×(N-D+1) and the time series expansion matrix Z∈R 13(D-1)×(N-D+1) .

[0055] Step (4): X D The column vectors from the 1st column to the N-D+1th column are taken as online data vectors in turn Implement time series difference feature analysis to calculate the corresponding anomaly detection indicators Q1, Q2, ..., Q N-D+1Afterwards, the control upper limit Q of the anomaly detection index is determined lim The specific implementation process is shown in steps (4.1) to (4.6).

[0056] Step (5): Obtain the process data x of the synthetic ammonia primary reformer at the latest sampling time point by using the DCS t ∈R 13×1 Afterwards, the same standardization processing as in step (2) is performed on it, so as to obtain the online data vector x

[0057] Step (6): Merge the online data vector x corresponding to the previous sampling time point, the online data vector x corresponding to the previous two sampling time points, and the online data vector x corresponding to the previous D-1 sampling time points to form an online time sequence vector x

[0058] Step (7): Set and Afterwards, the time sequence difference feature analysis is performed on the online data vector x , and the specific implementation process is shown in steps (A) to (D), so as to obtain the corresponding online load vector w t and the online coefficient vector β t , and further calculate the anomaly detection index Q

[0059] Step (8): Determine whether Q t is greater than Q lim ; if not, the synthetic ammonia primary reformer is in normal operation, and the step (5) is returned to continue performing anomaly detection by using the process data at the latest sampling time point; if yes, the step (9) is executed to determine whether to trigger an abnormal alarm.

[0060] Step (9): The step (5) is returned to continue performing anomaly detection by using the process data at the latest sampling time point, and if the anomaly detection indexes at the continuous A sampling time points are all greater than Q lim , an abnormal alarm is triggered; otherwise, the synthetic ammonia primary reformer is in normal operation.

Claims

1. A method for detecting abnormal operation of a synthetic ammonia primary furnace based on time series difference feature extraction, characterized in that: The specific steps include the following: Step (1): Obtain N sets of process data x1, x2, ..., x2 of the first stage of the synthetic ammonia furnace under normal operation from the DCS historical database. N ; where x i ∈R 13×1 Represents the i-th group of process data, with subscript i∈{1, 2, ..., N}, R 13×1 Represents a 13×1 dimensional real number vector, where each set of process data includes 13 data; Step (2): N sets of process data x1, x2, ..., x N The training data matrix X = [x1, x2, ..., x N ], the row vectors of each row in X are standardized to obtain the reference data matrix The standardization method is to subtract the mean of the row vector from the row vector and then divide it by the standard deviation of the row vector. 13×N represents a 13×N-dimensional real matrix, and R represents a real number set; Step (3): After setting the time series correlation order equal to D, use the following formula And the following formula ① respectively form the timing matrix X D ∈R 13×(N-D+1) and the time series expansion matrix Z∈R 13(D-1)×(N-D+1) : in, The reference data matrix is ​​represented in turn The column vectors of the 1st, 2nd, to Nth columns in , express The column vector of the D-2th column in , express The column vector of the D-1th column in , express The column vector of the N-2th column in express The column vector of the N-D+1th column in R 13(D-1)×(N-D+1) represents a real matrix of 13(D-1)×(N-D+1), R 13 ×(N-D+1) Represents a real matrix of 13×(N-D+1) dimensions; Step (4): X D The column vectors from the 1st column to the N-D+1th column are taken as online data vectors in turn Implement time series difference feature analysis to calculate the corresponding anomaly detection indicators Q1, Q2, ..., Q N-D+1 Then determine the upper control limit Q of the anomaly detection index lim , the specific implementation process is shown in steps (4.1) to (4.6); Step (4.1): After initializing i=1, set Equal to X D The column vector of column i in , sets the online time series vector z t is equal to the column vector of the i-th column in Z; Step (4.2): X D The remaining ND column vectors except the i-th column vector constitute the reference timing matrix The remaining ND column vectors of Z except the i-th column vector form the reference timing expansion matrix Step (4.3): Perform time series difference feature analysis to obtain the online load vector w t and the online coefficient vector β t ; Step (4.4): According to the formula Calculate the anomaly detection indicator Q i ;in, Online Timing Error Step (4.5): Determine whether i is less than N-D+1; if so, set i=i+1 and then set Equal to X D The column vector of column i in , and set z t If it is equal to the column vector of the i-th column in Z, return to step (4.2); if not, obtain N-D+1 anomaly detection indicators Q1, Q2, ..., Q N-D+1 Then, execute step (4.6); Step (4.6): For Q1, Q2, ..., Q N-D+1 Implement kernel density estimation to obtain the confidence limit Q at confidence level α = 99.5% KDE , and Q1, Q2, ..., Q N-D+1 The average value of the largest anomaly detection index ξ is recorded as Q AVG Then determine the upper control limit Q lim =max{Q KDE , Q AVG }; where ξ represents the largest integer not greater than n×99%, n=N-D+1, max{Q KDE , Q AVG } means taking Q KDE and Q AVG The maximum value in ; Step (5): Use DCS to obtain the process data x of the first stage of synthetic ammonia furnace at the latest sampling time t ∈R 13×1 Then, the same normalization process as in step (2) is performed to obtain the online data vector corresponding to the latest sampling time. Step (6): The online data vector corresponding to the previous sampling moment Online data vectors corresponding to the first two sampling moments Online data vector corresponding to the previous D-1 sampling moments Merge to form an online time series vector Step (7): Setup and Then, the online data vector Perform time series difference feature analysis to obtain the corresponding online load vector w t and the online coefficient vector β t , and then according to the formula Calculate the anomaly detection index Q at the latest sampling time t ; Step (8): Determine Q t Is it greater than Q lim If not, the ammonia synthesis furnace is operating normally, and the process returns to step (5) to continue to use the process data at the latest sampling time to perform abnormality detection; if so, step (9) is executed to decide whether to trigger an abnormality alarm; Step (9): Return to step (5) and continue to use the process data at the latest sampling moment to perform anomaly detection. If the anomaly detection index of A consecutive sampling moments is greater than Q lim , an abnormal alarm is triggered; otherwise, the synthetic ammonia stage furnace operates normally; where A is equal to the quotient of the shortest time required to trigger an abnormal alarm divided by the sampling interval.

2. The method for detecting abnormal operation of a synthetic ammonia primary furnace based on time series difference feature extraction according to claim 1, characterized in that: Each set of process data obtained using DCS is specifically a 13×1-dimensional real number vector consisting of m=13 data. The arrangement order of these 13 data is as follows: natural gas intake flow rate, exhaust gas intake flow rate, pressure and temperature at the heat exchanger outlet, natural gas temperature at the preheater outlet, furnace gas pressure at the outlet of a first-stage furnace, gas temperature at the inlet of a first-stage furnace, gas temperature at the upper left of a first-stage furnace, gas temperature at the upper right of a first-stage furnace, mixed furnace gas temperature at the top of a first-stage furnace, conversion gas temperature at the left outlet of a first-stage furnace, conversion gas temperature at the right outlet of a first-stage furnace, and conversion gas temperature at the outlet of a first-stage furnace.

3. The method for detecting abnormal operation of a synthetic ammonia primary furnace based on time series difference feature extraction according to claim 1, characterized in that: The specific implementation process of implementing the time series difference feature analysis in step (4.4) and step (7) is as follows: Step (A): Initialize the feature vector g∈R m×1 Equal to an arbitrary m×1-dimensional real vector; Step (B): According to the formula Calculate the online coefficient vector β t ; Where the superscript T represents the transpose of a matrix or vector, and the matrix L w and L D As shown in formula ②: In the above formula, Indicates calculation I D The Kronecker product between I and g, D represents the (D-1)×(D-1)-dimensional identity matrix; Step (C): Solve the generalized eigenvalue problem G = λL β After obtaining the eigenvector g corresponding to the maximum eigenvalue λ in g, use the formula Normalize g; where matrix G and matrix L β As shown in formula ③: In the above formula, Indicates the calculation of β t and I M The Kronecker product between M Represents the 13×13 dimensional identity matrix; Step (D): Determine whether the eigenvector g converges; if not, return to step (B); if so, set the online load vector w t =g, the time series difference characteristic analysis process is completed, and the online load vector w is obtained t and the online coefficient vector β t .

4. The method for detecting abnormal operation of a synthetic ammonia primary furnace based on time series difference feature extraction according to claim 3 is characterized in that: In step (C), the generalized eigenvalue problem G=λL is solved. β The specific implementation process of the eigenvector g corresponding to the maximum eigenvalue λ in g is as follows: Step 1: Matrix L β Implement singular value decomposition L β =UΛU T , only keep all values ​​greater than 10 -4 The non-zero singular values ​​of , thus obtaining the unitary matrix U and the singular value diagonal matrix Λ, and then calculating the matrix Step 2: According to the formula Calculate the matrix G β After that, initialize the eigenvector p to an arbitrary 13×1 dimensional real vector; Step (3): According to the formula p=G β p and Update to obtain the feature vector p; Step (4): Determine whether the eigenvector p has converged; if not, return to step (3); if so, obtain the final eigenvector p and execute step (5); Step 5: According to the formula Calculate the eigenvector g.

Images