CN113254874B - Uncertainty non-stationary industrial process oriented anomaly monitoring method - Google Patents

Uncertainty non-stationary industrial process oriented anomaly monitoring method Download PDF

Info

Publication number
CN113254874B
CN113254874B CN202110397942.8A CN202110397942A CN113254874B CN 113254874 B CN113254874 B CN 113254874B CN 202110397942 A CN202110397942 A CN 202110397942A CN 113254874 B CN113254874 B CN 113254874B
Authority
CN
China
Prior art keywords
stationary
equation
sample
data
uncertainty
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110397942.8A
Other languages
Chinese (zh)
Other versions
CN113254874A (en
Inventor
周东华
吴德浩
陈茂银
纪洪泉
钟麦英
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong University of Science and Technology
Original Assignee
Shandong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong University of Science and Technology filed Critical Shandong University of Science and Technology
Priority to CN202110397942.8A priority Critical patent/CN113254874B/en
Publication of CN113254874A publication Critical patent/CN113254874A/en
Application granted granted Critical
Publication of CN113254874B publication Critical patent/CN113254874B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures

Abstract

The invention discloses an abnormality monitoring method for an uncertain non-stationary industrial process, and particularly relates to the field of industrial process abnormality monitoring. The invention provides a probabilistic stationary subspace analysis method based on the stationary subspace analysis method by considering the process uncertainty. The method models the uncertainty explicitly and effectively separates non-stationary trends from process uncertainty. In consideration of mutual coupling among model parameters, the method utilizes an expectation maximization algorithm to conduct parameter decoupling and deduces a closed-form solution of iterative updating. Based on the model, two detection indexes are provided for anomaly monitoring under a probability framework. Compared with the existing abnormal monitoring method for the non-stationary process, the method provided by the invention eliminates the influence of uncertainty of the process and improves the detection capability of tiny faults in the non-stationary process; and the problem of overfitting of model parameters can be avoided, and a more accurate generative model is established for non-stationary data.

Description

Uncertainty non-stationary industrial process oriented anomaly monitoring method
Technical Field
The invention belongs to the field of industrial process abnormity monitoring, and particularly relates to an abnormity monitoring method for an uncertain non-stable industrial process.
Background
Anomaly monitoring is important for ensuring normal and efficient operation of industrial processes and equipment. Practical industrial processes often exhibit significant non-stationary characteristics, namely: the statistical properties of the process data may change over time. Non-stationary characteristics may be caused by a variety of factors, including raw material variations, load fluctuations, and equipment aging. The non-stationary characteristic affects the application of traditional anomaly monitoring methods such as principal component analysis and the like, and two types of errors often occur in the conventional anomaly monitoring methods. The first type of error is false alarm, that is, the influence of fault is covered by non-stationary trend, thereby resulting in a large amount of false alarm; the second category of errors is false positives, because conventional methods have difficulty tracking changes that are not stationary trends, resulting in higher false positives.
Over the past decade, non-stationary industrial process anomaly monitoring methods have grown in length. To our knowledge, these methods can be mainly divided into four categories, namely: the self-adaptive updating method is based on a method of co-integration analysis, a method of subspace decomposition and a method of trend analysis. Taking recursive principal component analysis as an example, the self-adaptive updating method establishes a data model according to training data, and then updates the parameters of the model by using normal test data. Second, the method based on co-integration analysis wants to find its stationary linear combination among the non-stationary variables. Subspace decomposition-based methods typically separate the stationary subspace from the full space, of which stationary subspace analysis methods are representative. Finally, trend analysis based methods typically extract trend information from non-stationary processes to accomplish the monitoring task.
The above methods are well applicable in anomaly monitoring of non-stationary industrial processes, however they do not take into account the uncertainty issues that exist in non-stationary processes. Actual industrial processes often suffer from various uncertainties, which may be caused by factors such as random noise and unknown disturbances. Process uncertainty presents a challenge to anomaly monitoring of non-stationary processes. On the one hand, the actual non-stationary trend is partially masked by process uncertainties, which may degrade the monitoring performance of the algorithm, especially for minor faults. On the other hand, it is difficult to distinguish true non-stationary trends from changes due to uncertainty, which leads to over-fitting problems of the model parameters. Probabilistic models provide a new view to these challenges, which strikes a good balance between the actual trends of real interest and the changes due to uncertainty. In fact, the measurement data of an industrial process is naturally presented in a statistical form, not in a deterministic manner. Currently, there are some probabilistic latent variable models that have been applied to anomaly monitoring of industrial processes, but they focus on stationary processes, not non-stationary processes.
Disclosure of Invention
In order to solve the problems, the invention provides an abnormity monitoring method for an uncertain non-stationary industrial process, which eliminates the influence of process uncertainty and improves the detection capability of tiny faults in the non-stationary process.
The technical scheme of the invention is as follows:
an abnormity monitoring method facing to uncertain non-stationary industrial process comprises an off-line training stage and an on-line monitoring stage; wherein the content of the first and second substances,
A. the off-line training stage comprises the following steps:
A1. collecting historical data of equipment operation in normal working condition of non-stationary process
Figure BDA0003019269340000021
Wherein N is the number of samples in the historical data set, and m is the number of measurement variables;
A2. the equipment operation historical data is standardized to enable each variable to be in a zero-mean and unit variance standardized form, and the standardization method is shown as a formula (1):
Figure BDA0003019269340000022
where Λ is a diagonal matrix, the diagonal elements consist of the standard deviations of m variables, 1 is a column vector consisting of m 1's,
Figure BDA0003019269340000023
is a sample mean vector of historical data;
A3. carrying out Johansen test on the standardized data matrix X, and determining the number d of stable components, wherein d is more than or equal to 1 and less than or equal to m-1;
A4. the k-th sample X (k) in the normalized data matrix X is decomposed into a probabilistic stationary subspace analysis model shown in equation (2):
Figure BDA0003019269340000024
wherein the content of the first and second substances,
Figure BDA0003019269340000025
is a mixing matrix that is reversible and that,
Figure BDA0003019269340000026
consists of the preceding d columns of A,
Figure BDA0003019269340000027
consists of the rear m-d columns of A,
Figure BDA0003019269340000028
is a stationary component, sn(k) Is a non-stationary component and is,
Figure BDA0003019269340000029
representing process uncertainty, which is independent of s (k), ImIs an identity matrix;
A5. the samples in the normalized data matrix X are divided into n consecutive, non-overlapping data segments, each data segment being denoted as
Figure BDA00030192693400000210
The number of samples of each data segment is recorded as NiAnd is provided with
Figure BDA00030192693400000211
A6. Setting a stationary component ss(k) And non-stationary component sn(k) Is not related, and ss(k) Gaussian distribution obeying a d dimension
Figure BDA00030192693400000212
Is invariant on a time scale, sn(k) Gaussian distribution obeying an m-d dimension
Figure BDA00030192693400000213
It is in each numberIs variable across segments;
A7. combining a4 and a6, the prior distribution of s (k) is:
Figure BDA00030192693400000214
then, the probability density function for x (k) is:
Figure BDA0003019269340000031
A8. estimating parameters of a probabilistic stationary subspace model using an expectation-maximization algorithm
Figure BDA0003019269340000032
A9. For the kth sample x (k), its corresponding local component is calculated by equation (5):
Figure BDA0003019269340000033
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
A10. Calculating the similarity between the sample x (k) and the sample belonging to the ith data segment, as shown in equation (6):
Figure BDA0003019269340000034
wherein the content of the first and second substances,
Figure BDA0003019269340000035
as shown in equation (4), the prior probability
Figure BDA0003019269340000036
Is set to 1/n, and
Figure BDA0003019269340000037
A11. the estimated contribution from samples x (k) is calculated by equation (7):
Figure BDA0003019269340000038
A12. combining equations (2) and (7), the estimated process uncertainty is shown as equation (8):
Figure BDA0003019269340000039
A13. due to the fact that
Figure BDA00030192693400000310
Is a stationary time series, and is therefore designed
Figure BDA00030192693400000311
Statistics monitor changes in process uncertainty, sample x (k) corresponds to
Figure BDA00030192693400000312
The statistic is shown as equation (9):
Figure BDA00030192693400000313
A14. for estimated local components
Figure BDA00030192693400000314
The mean vector and covariance matrix are calculated by equations (10) and (11), respectively:
Figure BDA00030192693400000315
Figure BDA00030192693400000316
wherein R isi＝AΣiAT+σ2Im，i＝1,2,…,n；
A15. Will be provided with
Figure BDA0003019269340000041
Is recorded as
Figure BDA0003019269340000042
Then
Figure BDA0003019269340000043
The mean vector and covariance matrix of the corresponding stationary components are calculated by equations (12) and (13), respectively:
Figure BDA0003019269340000044
Figure BDA0003019269340000045
wherein, W1From an identity matrix ImThe first d columns of (1);
A16. based on
Figure BDA0003019269340000046
Defining a local mahalanobis distance index as shown in equation (14):
Figure BDA0003019269340000047
A17. designing based on the weighted Mahalanobis distance as shown in equation (15)
Figure BDA0003019269340000048
Statistics to monitor the change in stationary components:
Figure BDA0003019269340000049
A18. given a confidence level α, statistics are determined using a kernel density estimation method
Figure BDA00030192693400000410
And
Figure BDA00030192693400000411
are respectively recorded as
Figure BDA00030192693400000412
And
Figure BDA00030192693400000413
B. the on-line monitoring stage comprises the following steps:
B1. acquiring real-time operation data y (t) of a non-stationary process of the equipment, and standardizing y (t):
Figure BDA00030192693400000414
wherein, yiAnd xiThe ith variable of the real-time data before and after normalization respectively,
Figure BDA00030192693400000415
is the mean of the ith variable and is,
Figure BDA00030192693400000416
is the standard deviation of the ith variable;
B2. for the real-time sample x (t) after normalization, its corresponding local component is calculated by equation (17):
Figure BDA00030192693400000417
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
B3. Calculating the similarity between the sample x (t) and the samples belonging to the ith data segment, as shown in equation (18):
Figure BDA00030192693400000418
wherein the content of the first and second substances,
Figure BDA00030192693400000419
prior probability
Figure BDA00030192693400000420
Is set to 1/n, and
Figure BDA00030192693400000421
B4. the estimated contribution from samples x (t) is calculated by equation (19):
Figure BDA0003019269340000051
B5. combining equations (2) and (19), the estimated process uncertainty is shown as equation (20):
Figure BDA0003019269340000052
B6. calculating the real-time sample x (t) corresponds to
Figure BDA0003019269340000053
Statistics:
Figure BDA0003019269340000054
B7. based on
Figure BDA0003019269340000055
Calculating the local mahalanobis distance index as shown in equation (22):
Figure BDA0003019269340000056
B8. calculating the real-time sample x (t) corresponds to
Figure BDA0003019269340000057
Statistics:
Figure BDA0003019269340000058
B9. will make statistics of
Figure BDA0003019269340000059
And
Figure BDA00030192693400000510
respectively in their control limits
Figure BDA00030192693400000511
And
Figure BDA00030192693400000512
by comparison, if
Figure BDA00030192693400000513
And is
Figure BDA00030192693400000514
Judging that the industrial process is in a normal operation condition at present, otherwise, judging that the operation is abnormal.
Preferably, step A8 specifically includes:
A801. initial setting of model parameters
Figure BDA00030192693400000515
Arbitrarily assigning an initial value;
A802. calculating the first moment and the second moment corresponding to s (k) according to the parameter value theta obtained in the last step:
Figure BDA00030192693400000516
〈s(k)s(k)T〉＝HiΣi+〈s(k)〉〈s(k)〉T (25)
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
A803. The model parameters Θ are updated using equations (24) and (25):
Figure BDA00030192693400000517
Figure BDA00030192693400000518
Figure BDA00030192693400000519
Figure BDA0003019269340000061
Figure BDA0003019269340000062
Figure BDA0003019269340000063
wherein, W1From an identity matrix ImThe first d columns of (A) constitute (W)2From an identity matrix ImThe last m-d columns of (A);
A804. if the absolute value of the difference between the norm values of theta obtained by two iterations is less than 10-5Stopping iteration and outputting the estimated value of the model parameter theta, otherwise returning to the step A802 for next iteration.
The invention has the following beneficial technical effects:
the method provided by the invention eliminates the influence of process uncertainty and improves the detection capability of tiny faults in a non-stationary process; and the problem of overfitting of model parameters can be avoided, and a more accurate generative model is established for non-stationary data.
Drawings
FIG. 1 is a flow chart of an anomaly monitoring method for an uncertain non-stationary industrial process according to the present invention;
FIG. 2 is a schematic diagram of a closed-loop controlled continuous stirred tank reactor according to an embodiment of the invention;
FIG. 3 shows a synergistic analysis-T in example 1 of the present invention2A graph of the monitored results of the statistics;
FIG. 4 is a graph showing the results of monitoring the Q statistic of the recursive principal component analysis in example 1 of the present invention;
FIG. 5 shows T of recursive principal component analysis in example 1 of the present invention2A graph of the monitored results of the statistics;
FIG. 6 is a graph showing the results of stationary subspace analysis-Mahalanobis distance monitoring in example 1 of the present invention;
FIG. 7 is T of probability stationary subspace analysis in embodiment 1 of the present inventione 2A graph of the monitored results of the statistics;
FIG. 8 is T of probability stationary subspace analysis in embodiment 1 of the present inventions 2A graph of the monitored results of the statistics;
FIG. 9 shows a co-integration analysis-T in example 2 of the present invention2A graph of the monitored results of the statistics;
FIG. 10 is a graph showing the results of monitoring the Q statistic of the recursive principal component analysis in example 2 of the present invention;
FIG. 11 shows T of recursive principal component analysis in embodiment 2 of the present invention2A graph of the monitored results of the statistics;
FIG. 12 is a graph showing the results of stationary subspace analysis-Mahalanobis distance monitoring in example 2 of the present invention;
FIG. 13 is T of probability stationary subspace analysis in embodiment 2 of the present inventione 2Of statisticsA monitoring result graph;
FIG. 14 is T of probability stationary subspace analysis in embodiment 2 of the present inventions 2And (5) a monitoring result graph of the statistics.
Detailed Description
The invention is described in further detail below with reference to the following figures and detailed description:
as shown in fig. 1, an anomaly monitoring method for uncertain non-stationary industrial processes includes an offline training phase and an online monitoring phase; wherein the content of the first and second substances,
A. the off-line training stage comprises the following steps:
A1. collecting historical data of equipment operation in normal working condition of non-stationary process
Figure BDA0003019269340000071
Wherein N is the number of samples in the historical data set, and m is the number of measurement variables;
A2. the equipment operation historical data is standardized to enable each variable to be in a zero-mean and unit variance standardized form, and the standardization method is shown as a formula (1):
Figure BDA0003019269340000072
where Λ is a diagonal matrix, the diagonal elements consist of the standard deviations of m variables, 1 is a column vector consisting of m 1's,
Figure BDA0003019269340000073
is a sample mean vector of historical data;
A3. carrying out Johansen test on the standardized data matrix X, and determining the number d of stable components, wherein d is more than or equal to 1 and less than or equal to m-1;
A4. it is assumed that the kth sample X (k) in the normalized data matrix X can be decomposed into the form shown in equation (2):
Figure BDA0003019269340000074
this model is called a probabilistic stationary subspace analysis model, where,
Figure BDA0003019269340000075
is a mixing matrix that is reversible and that,
Figure BDA0003019269340000076
consists of the preceding d columns of A,
Figure BDA0003019269340000077
consists of the rear m-d columns of A,
Figure BDA0003019269340000078
is a stationary component, sn(k) Is a non-stationary component and is,
Figure BDA0003019269340000079
representing process uncertainty, which is independent of s (k), ImIs an identity matrix;
A5. the samples in the normalized data matrix X are divided into n consecutive, non-overlapping data segments, each data segment being denoted as
Figure BDA00030192693400000710
The number of samples of each data segment is recorded as NiAnd is provided with
Figure BDA00030192693400000711
A6. Assuming a stationary component ss(k) And non-stationary component sn(k) Is not related, and ss(k) Gaussian distribution obeying a d dimension
Figure BDA00030192693400000712
Is invariant on a time scale, sn(k) Gaussian distribution obeying an m-d dimension
Figure BDA00030192693400000713
It varies across data segments;
A7. the prior distribution of s (k) is known from the joint assumptions in A4 and A6 as:
Figure BDA00030192693400000714
then, the probability density function for x (k) is:
Figure BDA0003019269340000081
A8. estimating parameters of a probabilistic stationary subspace model using an expectation-maximization algorithm
Figure BDA0003019269340000082
A9. For the kth sample x (k), its corresponding local component is calculated by equation (5):
Figure BDA0003019269340000083
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
A10. Calculating the similarity between the sample x (k) and the sample belonging to the ith data segment, as shown in equation (6):
Figure BDA0003019269340000084
wherein the content of the first and second substances,
Figure BDA0003019269340000085
as shown in equation (4), the prior probability
Figure BDA0003019269340000086
Is set to 1/n, and
Figure BDA0003019269340000087
A11. the estimated contribution from samples x (k) can be calculated by equation (7):
Figure BDA0003019269340000088
A12. combining equations (2) and (7), the estimated process uncertainty is shown as equation (8):
Figure BDA0003019269340000089
A13. due to the fact that
Figure BDA00030192693400000810
Is a smooth time sequence and can be designed accordingly
Figure BDA00030192693400000811
Statistics monitor changes in process uncertainty, sample x (k) corresponds to
Figure BDA00030192693400000812
The statistic is shown as equation (9):
Figure BDA00030192693400000813
A14. for estimated local components
Figure BDA00030192693400000814
The mean vector and covariance matrix can be calculated by equations (10) and (11), respectively:
Figure BDA00030192693400000815
Figure BDA00030192693400000816
wherein R isi＝AΣiAT+σ2Im，i＝1,2,…,n；
A15. Will be provided with
Figure BDA00030192693400000817
Is recorded as
Figure BDA00030192693400000818
Then
Figure BDA00030192693400000819
The mean vector and covariance matrix of the corresponding stationary components can be calculated by equations (12) and (13), respectively:
Figure BDA0003019269340000091
Figure BDA0003019269340000092
wherein, W1From an identity matrix ImThe first d columns of (1);
A16. based on
Figure BDA0003019269340000093
Defining a local mahalanobis distance index as shown in equation (14):
Figure BDA0003019269340000094
A17. designing based on the weighted Mahalanobis distance as shown in equation (15)
Figure BDA0003019269340000095
Statistics to monitor the change in stationary components:
Figure BDA0003019269340000096
A18. given a confidence level α, statistics are determined using a kernel density estimation method
Figure BDA0003019269340000097
And
Figure BDA0003019269340000098
are respectively recorded as
Figure BDA0003019269340000099
And
Figure BDA00030192693400000910
B. the on-line monitoring stage comprises the following steps:
B1. acquiring real-time operation data y (t) of a non-stationary process of the equipment, and standardizing y (t):
Figure BDA00030192693400000911
wherein, yiAnd xiThe ith variable of the real-time data before and after normalization respectively,
Figure BDA00030192693400000912
is the mean of the ith variable and is,
Figure BDA00030192693400000913
is the standard deviation of the ith variable;
B2. for the real-time sample x (t) after normalization, its corresponding local component is calculated by equation (17):
Figure BDA00030192693400000914
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
B3. Calculating the similarity between the sample x (t) and the samples belonging to the ith data segment, as shown in equation (18):
Figure BDA00030192693400000915
wherein the content of the first and second substances,
Figure BDA00030192693400000916
prior probability
Figure BDA00030192693400000917
Is set to 1/n, and
Figure BDA00030192693400000918
B4. the estimated contribution from samples x (t) can be calculated by equation (19):
Figure BDA0003019269340000101
B5. combining equations (2) and (19), the estimated process uncertainty is shown as equation (20):
Figure BDA0003019269340000102
B6. calculating the real-time sample x (t) corresponds to
Figure BDA0003019269340000103
Statistics:
Figure BDA0003019269340000104
B7. based on
Figure BDA0003019269340000105
Calculating the local mahalanobis distance index as shown in equation (22):
Figure BDA0003019269340000106
B8. calculating the real-time sample x (t) corresponds to
Figure BDA0003019269340000107
Statistics:
Figure BDA0003019269340000108
B9. will make statistics of
Figure BDA0003019269340000109
And
Figure BDA00030192693400001010
respectively in their control limits
Figure BDA00030192693400001011
And
Figure BDA00030192693400001012
by comparison, if
Figure BDA00030192693400001013
And is
Figure BDA00030192693400001014
Judging that the industrial process is in a normal operation condition at present, otherwise, judging that the operation is abnormal.
The specific content of the step A8 includes:
A801. initial setting of model parameters
Figure BDA00030192693400001015
Arbitrarily assigning an initial value;
A802. calculating the first moment and the second moment corresponding to s (k) according to the parameter value theta obtained in the last step:
Figure BDA00030192693400001016
<s(k)s(k)T>＝HiΣi+<s(k)〉〈s(k)>T (25)
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
A803. The model parameters Θ are updated using equations (24) and (25):
Figure BDA00030192693400001017
Figure BDA00030192693400001018
Figure BDA00030192693400001019
Figure BDA00030192693400001020
Figure BDA0003019269340000111
Figure BDA0003019269340000112
wherein, W1From an identity matrix ImThe first d columns of (A) constitute (W)2From an identity matrix ImThe last m-d columns of (A);
A804. if the absolute value of the difference between the norm values of theta obtained by two iterations is less than 10-5Stopping iteration and outputting the estimated value of the model parameter theta, otherwise returning to the step A802 for the next iterationAnd (4) generation.
In order to help understand the invention and simultaneously visually show the effect of the method for monitoring the abnormal state of the uncertain non-stationary industrial process, the following description is based on two embodiments. The data for both examples was derived from a closed loop controlled continuous stirred tank reactor which provides a library of standard models widely used in the field of monitoring of industrial process anomalies. A schematic of the reactor is shown in fig. 2, where a first order exothermic reaction is performed. Its dynamic model is shown in equation (32):
Figure BDA0003019269340000113
wherein C is the reaction concentration, T is the reaction temperature, TcIs the temperature of cold water, QcIs the cold water flow, viIs measuring noise, Ci、TiAnd TciConstituting the system input.
Since the original continuous stirred tank reactor is a stationary process, system inputs need to be modified in order to simulate non-stationary characteristics. In the modified version, TciIs switched to a steady state, CiAnd TiAdding a non-stationary drift term lambda:
Figure BDA0003019269340000114
wherein the content of the first and second substances,
Figure BDA0003019269340000115
and
Figure BDA0003019269340000116
respectively represent CiAnd TiThe nominal value of (a) is,
Figure BDA0003019269340000117
and
Figure BDA0003019269340000118
representing two gaussian random variables. Thus, the continuous stirred tank reactor exhibits non-stationary characteristics.
In the off-line training phase, 1200 normal samples were collected with a1 minute sampling interval. Each sample consists of 7 measured variables, i.e.: x ═ Ci,Ti,C,T,Qc,Tci,Tc]T. In addition to the method proposed by the present invention, three other monitoring algorithms are also added for comparison, respectively, a method based on co-integration analysis, a recursive principal component analysis method and a method based on stationary subspace analysis. Wherein the method based on the co-integration analysis combines the co-integration analysis with T2The statistics are combined, and the co-integration order is selected as r-5. For the recursive pivot analysis method, the pivot is selected based on an accumulated contribution rate of 85%. For stationary subspace analysis and probabilistic stationary subspace analysis, the number of stationary components is set to d-5, and the number of data segment partitions is set to n-19. In addition, the iteration error of the probabilistic stationary subspace analysis method is assumed to be ∈ 1 × 10-6. Given a confidence level α of 0.99, the control limits for all monitoring statistics are determined by the kernel density estimation method.
TABLE 1 two typical failures in a continuous stirred tank reactor
Figure BDA0003019269340000121
In both embodiments of the invention, two typical failures in a continuous stirred tank reactor were considered separately, including multiplicative and additive failures, as shown in table 1. For both faults, 1200 samples were collected for performance comparison of the algorithm. The fault is introduced starting from the 601 st sample and continuing until the end.
Example 1
The heat transfer coefficient degradation failure was tested in example 1 and the monitoring results of the different algorithms are summarized in fig. 3-8. The fault is a typical tiny fault, and three methods based on the co-integration analysis, the stationary subspace analysis and the probabilistic stationary subspace analysis all show the false negative to different degrees in the early stage of the fault, wherein the three methods are respectively based on the co-integration analysis, the stationary subspace analysis and the probabilistic stationary subspace analysisProbabilistic stationary subspace analysis
Figure BDA0003019269340000122
The statistic has the shortest detection delay as shown in fig. 7.
Figure BDA0003019269340000123
The statistic has the advantage of containing non-stationary information. In addition, as shown in fig. 4 and 5, the recursive principal component analysis method alarms early in the early stage of a fault, but leaves many false alarms as the magnitude of the fault increases. This reflects the failure of the recursive principal component analysis method to discern subtle faults from non-stationary trends.
Example 2
An additive sensor constant offset fault was tested in example 2, which occurred in magnitude 1 in the reaction temperature variation. This is a minor fault of small magnitude compared to the nominal value of the reaction temperature of 430 c. In this example, the monitoring results of the four different methods are shown in fig. 9-14. As can be seen from fig. 9 and 12, the method based on both the co-integration analysis and the stationary subspace analysis can only partially detect the fault, but leaves a certain degree of false negative. Neither of these methods models non-stationary components, thus resulting in unsatisfactory monitoring performance for minor faults. In addition, the recursive principal component analysis method can hardly detect the fault (see fig. 10 and 11) because the fault magnitude is too small to be masked by the non-stationary trend. For probabilistic stationary subspace analysis, it
Figure BDA0003019269340000124
The detection effect of statistics is similar to that of the co-integration analysis and stationary subspace analysis, and
Figure BDA0003019269340000125
the statistics are able to effectively detect the fault at a detection rate of 97.8%. This reflects that the true stationary trend is separated from the process uncertainty, which can effectively improve the algorithm's detection capability for minor faults.
It is to be understood that the above description is not intended to limit the present invention, and the present invention is not limited to the above examples, and those skilled in the art may make modifications, alterations, additions or substitutions within the spirit and scope of the present invention.

Claims (2)

1. An abnormity monitoring method facing to uncertain non-stationary industrial process is characterized by comprising an off-line training stage and an on-line monitoring stage; wherein the content of the first and second substances,
A. the off-line training stage comprises the following steps:
A1. collecting historical data of equipment operation in normal working condition of non-stationary process
Figure FDA0003019269330000011
Wherein N is the number of samples in the historical data set, and m is the number of measurement variables;
A2. the equipment operation historical data is standardized to enable each variable to be in a zero-mean and unit variance standardized form, and the standardization method is shown as a formula (1):
Figure FDA0003019269330000012
where Λ is a diagonal matrix, the diagonal elements consist of the standard deviations of m variables, 1 is a column vector consisting of m 1's,
Figure FDA0003019269330000013
is a sample mean vector of historical data;
A3. carrying out Johansen test on the standardized data matrix X, and determining the number d of stable components, wherein d is more than or equal to 1 and less than or equal to m-1;
A4. the k-th sample X (k) in the normalized data matrix X is decomposed into a probabilistic stationary subspace analysis model shown in equation (2):
Figure FDA0003019269330000014
wherein the content of the first and second substances,
Figure FDA0003019269330000015
is a mixing matrix that is reversible and that,
Figure FDA0003019269330000016
consists of the preceding d columns of A,
Figure FDA0003019269330000017
consists of the rear m-d columns of A,
Figure FDA0003019269330000018
is a stationary component, sn(k) Is a non-stationary component and is,
Figure FDA0003019269330000019
representing process uncertainty, which is independent of s (k), ImIs an identity matrix;
A5. the samples in the normalized data matrix X are divided into n consecutive, non-overlapping data segments, each data segment being denoted as
Figure FDA00030192693300000110
The number of samples of each data segment is recorded as NiAnd is provided with
Figure FDA00030192693300000111
A6. Setting a stationary component ss(k) And non-stationary component sn(k) Is not related, and ss(k) Gaussian distribution obeying a d dimension
Figure FDA00030192693300000112
Is invariant on a time scale, sn(k) Gaussian distribution obeying an m-d dimension
Figure FDA00030192693300000113
It varies across data segments;
A7. combining a4 and a6, the prior distribution of s (k) is:
Figure FDA00030192693300000114
then, the probability density function for x (k) is:
Figure FDA0003019269330000021
A8. estimating parameters of a probabilistic stationary subspace model using an expectation-maximization algorithm
Figure FDA0003019269330000022
A9. For the kth sample x (k), its corresponding local component is calculated by equation (5):
Figure FDA0003019269330000023
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
A10. Calculating the similarity between the sample x (k) and the sample belonging to the ith data segment, as shown in equation (6):
Figure FDA0003019269330000024
wherein the content of the first and second substances,
Figure FDA0003019269330000025
as shown in equation (4), the prior probability
Figure FDA0003019269330000026
Is set to 1/n, and
Figure FDA0003019269330000027
A11. the estimated contribution from samples x (k) is calculated by equation (7):
Figure FDA0003019269330000028
A12. combining equations (2) and (7), the estimated process uncertainty is shown as equation (8):
Figure FDA0003019269330000029
A13. due to the fact that
Figure FDA00030192693300000210
Is a stationary time series, and is therefore designed
Figure FDA00030192693300000218
Statistics monitor changes in process uncertainty, sample x (k) corresponds to
Figure FDA00030192693300000219
The statistic is shown as equation (9):
Figure FDA00030192693300000211
A14. for estimated local components
Figure FDA00030192693300000212
The mean vector and covariance matrix are calculated by equations (10) and (11), respectively:
Figure FDA00030192693300000213
Figure FDA00030192693300000214
wherein R isi＝AΣiAT+σ2Im，i＝1,2,…,n；
A15. Will be provided with
Figure FDA00030192693300000215
Is recorded as
Figure FDA00030192693300000216
Then
Figure FDA00030192693300000217
The mean vector and covariance matrix of the corresponding stationary components are calculated by equations (12) and (13), respectively:
Figure FDA0003019269330000031
Figure FDA0003019269330000032
wherein, W1From an identity matrix ImThe first d columns of (1);
A16. based on
Figure FDA0003019269330000033
Defining a local mahalanobis distance index as shown in equation (14):
Figure FDA0003019269330000034
A17. based on weighted mahalanobis distanceAs shown in design formula (15)
Figure FDA00030192693300000316
Statistics to monitor the change in stationary components:
Figure FDA0003019269330000035
A18. given a confidence level α, statistics are determined using a kernel density estimation method
Figure FDA00030192693300000317
And
Figure FDA00030192693300000318
are respectively recorded as
Figure FDA0003019269330000036
And
Figure FDA0003019269330000037
B. the on-line monitoring stage comprises the following steps:
B1. acquiring real-time operation data y (t) of a non-stationary process of the equipment, and standardizing y (t):
Figure FDA0003019269330000038
wherein, yiAnd xiThe ith variable of the real-time data before and after normalization respectively,
Figure FDA0003019269330000039
is the mean of the ith variable and is,
Figure FDA00030192693300000310
is the standard deviation of the ith variable;
B2. for the real-time sample x (t) after normalization, its corresponding local component is calculated by equation (17):
Figure FDA00030192693300000311
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
B3. Calculating the similarity between the sample x (t) and the samples belonging to the ith data segment, as shown in equation (18):
Figure FDA00030192693300000312
wherein the content of the first and second substances,
Figure FDA00030192693300000313
prior probability
Figure FDA00030192693300000314
Is set to 1/n, and
Figure FDA00030192693300000315
B4. the estimated contribution from samples x (t) is calculated by equation (19):
Figure FDA0003019269330000041
B5. combining equations (2) and (19), the estimated process uncertainty is shown as equation (20):
Figure FDA0003019269330000042
B6. calculating the real-time sample x (t) corresponds to
Figure FDA00030192693300000416
Statistics:
Figure FDA0003019269330000043
B7. based on
Figure FDA0003019269330000044
Calculating the local mahalanobis distance index as shown in equation (22):
Figure FDA0003019269330000045
B8. calculating the real-time sample x (t) corresponds to
Figure FDA00030192693300000417
Statistics:
Figure FDA0003019269330000046
B9. will make statistics of
Figure FDA00030192693300000418
And
Figure FDA00030192693300000419
respectively in their control limits
Figure FDA0003019269330000047
And
Figure FDA0003019269330000048
by comparison, if
Figure FDA0003019269330000049
And is
Figure FDA00030192693300000410
Judging that the industrial process is in a normal operation condition at present, otherwise, judging that the operation is abnormal.
2. The anomaly monitoring method oriented to the uncertain non-stationary industrial processes as claimed in claim 1, wherein said step A8 specifically comprises:
A801. initial setting of model parameters
Figure FDA00030192693300000411
Arbitrarily assigning an initial value;
A802. calculating the first moment and the second moment corresponding to s (k) according to the parameter value theta obtained in the last step:
Figure FDA00030192693300000412
<s(k)s(k)T>＝HiΣi+<s(k)><s(k)>T (25)
wherein Hi＝I-Σi(σ2I+ATAΣi)-1ATA，i＝1,2,…,n；
A803. The model parameters Θ are updated using equations (24) and (25):
Figure FDA00030192693300000413
Figure FDA00030192693300000414
Figure FDA00030192693300000415
Figure FDA0003019269330000051
Figure FDA0003019269330000052
Figure FDA0003019269330000053
wherein, W1From an identity matrix ImThe first d columns of (A) constitute (W)2From an identity matrix ImThe last m-d columns of (A);
A804. if the absolute value of the difference between the norm values of theta obtained by two iterations is less than 10-5Stopping iteration and outputting the estimated value of the model parameter theta, otherwise returning to the step A802 for next iteration.
CN202110397942.8A 2021-04-14 2021-04-14 Uncertainty non-stationary industrial process oriented anomaly monitoring method Active CN113254874B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110397942.8A CN113254874B (en) 2021-04-14 2021-04-14 Uncertainty non-stationary industrial process oriented anomaly monitoring method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110397942.8A CN113254874B (en) 2021-04-14 2021-04-14 Uncertainty non-stationary industrial process oriented anomaly monitoring method

Publications (2)

Publication Number Publication Date
CN113254874A CN113254874A (en) 2021-08-13
CN113254874B true CN113254874B (en) 2022-04-15

Family

ID=77220674

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110397942.8A Active CN113254874B (en) 2021-04-14 2021-04-14 Uncertainty non-stationary industrial process oriented anomaly monitoring method

Country Status (1)

Country Link
CN (1) CN113254874B (en)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110942258A (en) * 2019-12-10 2020-03-31 山东科技大学 Performance-driven industrial process anomaly monitoring method
CN112040858A (en) * 2017-10-19 2020-12-04 脸谱科技有限责任公司 System and method for identifying biological structures associated with neuromuscular source signals
CN112069457A (en) * 2020-08-13 2020-12-11 山东科技大学 Non-stationary dynamic process abnormity monitoring method based on dynamic stationary subspace analysis

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10157319B2 (en) * 2017-02-22 2018-12-18 Sas Institute Inc. Monitoring, detection, and surveillance system using principal component analysis with machine and sensor data

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112040858A (en) * 2017-10-19 2020-12-04 脸谱科技有限责任公司 System and method for identifying biological structures associated with neuromuscular source signals
CN110942258A (en) * 2019-12-10 2020-03-31 山东科技大学 Performance-driven industrial process anomaly monitoring method
CN112069457A (en) * 2020-08-13 2020-12-11 山东科技大学 Non-stationary dynamic process abnormity monitoring method based on dynamic stationary subspace analysis

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Reduced Stationary Subspace Analysis for Anomaly Detection in Nonstationary Industrial Processes;Dehao Wu 等;《2020 Chinese Automation Congress》;20201108;第6612-6616页 *
风电机组智能状态评估与故障预测研究进展;李重桂等;《电站系统工程》;20200715(第04期);全文 *

Also Published As

Publication number Publication date
CN113254874A (en) 2021-08-13

Similar Documents

Publication Publication Date Title
Zhao et al. Fault-relevant principal component analysis (FPCA) method for multivariate statistical modeling and process monitoring
Chen et al. Dynamic process fault monitoring based on neural network and PCA
Lee et al. Statistical monitoring of dynamic processes based on dynamic independent component analysis
Yu A particle filter driven dynamic Gaussian mixture model approach for complex process monitoring and fault diagnosis
Jiang et al. Weighted kernel principal component analysis based on probability density estimation and moving window and its application in nonlinear chemical process monitoring
US7233884B2 (en) Methodology for temporal fault event isolation and identification
Chen et al. Cognitive fault diagnosis in tennessee eastman process using learning in the model space
Kariwala et al. A branch and bound method for isolation of faulty variables through missing variable analysis
Zhao et al. Online fault prognosis with relative deviation analysis and vector autoregressive modeling
CN108549908B (en) Chemical process fault detection method based on multi-sampling probability kernel principal component model
Harrou et al. A statistical fault detection strategy using PCA based EWMA control schemes
Wu et al. Locality preserving randomized canonical correlation analysis for real-time nonlinear process monitoring
Bakdi et al. An improved plant‐wide fault detection scheme based on PCA and adaptive threshold for reliable process monitoring: Application on the new revised model of Tennessee Eastman process
Jia et al. Dynamic higher-order cumulants analysis for state monitoring based on a novel lag selection
Zhang et al. Fault detection and diagnosis strategy based on a weighted and combined index in the residual subspace associated with PCA
Lan et al. Dynamic statistical process monitoring based on generalized canonical variate analysis
Shang et al. Efficient recursive canonical variate analysis approach for monitoring time‐varying processes
CN110751217A (en) Equipment energy consumption ratio early warning analysis method based on principal component analysis
Jiang et al. Probabilistic monitoring of chemical processes using adaptively weighted factor analysis and its application
CN116627116B (en) Process industry fault positioning method and system and electronic equipment
CN113254874B (en) Uncertainty non-stationary industrial process oriented anomaly monitoring method
Venkateswaran et al. Design of functional observers for fault detection and isolation in nonlinear systems in the presence of noises
Fei et al. Online process monitoring for complex systems with dynamic weighted principal component analysis
Zhang et al. A comparison of different statistics for detecting multiplicative faults in multivariate statistics-based fault detection approaches
Peng et al. Phase partition and fault diagnosis of batch process based on keca angular similarity

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant