CN109917777B - Fault detection method based on mixed multi-sampling rate probability principal component analysis model - Google Patents

Fault detection method based on mixed multi-sampling rate probability principal component analysis model Download PDF

Info

Publication number
CN109917777B
CN109917777B CN201910304064.3A CN201910304064A CN109917777B CN 109917777 B CN109917777 B CN 109917777B CN 201910304064 A CN201910304064 A CN 201910304064A CN 109917777 B CN109917777 B CN 109917777B
Authority
CN
China
Prior art keywords
sampling rate
mode
principal component
component analysis
training sample
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
CN201910304064.3A
Other languages
Chinese (zh)
Other versions
CN109917777A (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.)
Zhejiang Lover Health Science and Technology Development Co Ltd
Original Assignee
Zhejiang Lover Health Science and Technology Development Co Ltd
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 Zhejiang Lover Health Science and Technology Development Co Ltd filed Critical Zhejiang Lover Health Science and Technology Development Co Ltd
Priority to CN201910304064.3A priority Critical patent/CN109917777B/en
Publication of CN109917777A publication Critical patent/CN109917777A/en
Application granted granted Critical
Publication of CN109917777B publication Critical patent/CN109917777B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Abstract

The invention discloses a fault detection method based on a mixed multi-sampling rate probability principal component analysis model, which comprises the following steps: (1) collecting data of different sampling rates under multiple modes of normal operation of a chemical process to be monitored to form a training sample set for modeling; (2) preprocessing a training sample set; (3) constructing a mixed multi-sampling rate probability principal component analysis model by utilizing the preprocessed training sample set; (4) collecting monitoring samples of a sampling rate process of a chemical process to be monitored on line; (5) based on the constructed mixed multi-sampling rate probability principal component analysis model, calculating monitoring statistics of monitoring samples; (6) and judging whether the fault occurs according to the monitoring statistics. Compared with the traditional multi-mode modeling method, the multi-mode modeling method has the advantages that the multi-sampling rate data characteristics and the multi-mode process characteristics can be considered simultaneously, so that the applicability and the fault detection precision in the multi-mode multi-sampling rate process fault detection problem are improved.

Description

Fault detection method based on mixed multi-sampling rate probability principal component analysis model
Technical Field
The invention belongs to the technical field of fault detection, and particularly relates to a fault detection method based on a mixed multi-sampling rate probability principal component analysis model.
Background
With the development of modern industry, process safety and product quality are widely regarded. With the wide application of Distributed Control Systems (DCS) in the industrial field, a large number of process variables can be collected and stored by various sensors with high sampling rates, while key quality variables related to production safety and product quality need to be collected in a low sampling rate manner and acquired through testing, so that the multi-sampling rate characteristics of data and the difficult acquirement of important variable data are caused, which is a challenge to the management of actual industrial engineering. Meanwhile, with the continuous progress of process monitoring (MSPM) and soft measurement technology based on multivariate statistical analysis, the dimension reduction, reconstruction and visualization of mass data variables are realized, and the method is widely applied to the fields of pharmacy, chemical industry, pollution control and the like. Traditional static Principal Component Analysis (PCA) and least squares estimation (PLS) models can effectively extract cross-correlation of variables, but are not effective in the face of multi-modal problems. While techniques based on Gaussian Mixture Models (GMMs) and mixed probability models can effectively extract and analyze multi-modal data characteristics of the process, they cannot utilize complete multi-sample rate data. The method based on multi-sampling probability principal component analysis can completely utilize multi-sampling rate data information and effectively estimate model parameters by utilizing an expectation-maximization (EM) algorithm, but has poor processing effect on multi-modal data. Therefore, there is a need to provide a soft measurement technique that can not only fully utilize multi-sampling rate data information, but also fully consider the characteristics of multi-modal data in industrial processes.
Disclosure of Invention
The invention aims to provide a multi-mode multi-sampling rate industrial process fault detection method based on a mixed multi-sampling rate probability principal component analysis model aiming at the defects of the prior art.
Aiming at the problem of multi-modal working condition fault detection in the chemical process, the multi-sampling rate data under a plurality of normal working modes are collected by a discrete control system, and a mixed multi-sampling rate probability principal component analysis model is established. The model structure may be estimated by an expectation-maximization algorithm. On the basis, sampling is carried out on the data of the on-line pre-decarburization unit to obtain a multi-sampling rate test sample, latent variables of the test sample are estimated by utilizing the existing model structure, and statistical monitoring indexes at the moment are calculated to realize the fault diagnosis result of the pre-decarburization unit.
A fault detection method based on a mixed multi-sampling rate probability principal component analysis model comprises the following steps:
(1) collecting data of different sampling rates under M modes of normal operation of a chemical process to be monitored, wherein M is more than or equal to 1, and forming a training sample set X for modeling;
(2) preprocessing a training sample set X;
(3) constructing a mixed multi-sampling rate probability principal component analysis model by utilizing the preprocessed training sample set X;
(4) on-line collection of monitoring sample x of sampling rate process of chemical process to be monitorednew；
(5) Based on the constructed mixed multi-sampling rate probability principal component analysis model, calculating monitoring statistics of monitoring samples;
(6) and judging whether the fault occurs according to the monitoring statistics.
Preferably, in step (1) or/and step (4), the training sample or the monitoring sample is collected by using a discrete control system in the invention. The discrete control system is used as a mature data acquisition system, and the stability and accuracy of data can be integrally guaranteed.
The multi-mode process means that for a certain chemical unit, multiple products may respectively correspond to one working condition according to the requirements of different products, and the variable number and the sampling rate under each working condition may change; taking the pre-decarburization unit in the synthetic ammonia process as an example, multiple working conditions may exist according to different component contents of production raw materials or different requirements of final products, and the working conditions may have differences in parameters such as stable working condition operating points, reaction internal parameters, final product quality indexes and the like. Meanwhile, the difference of sampling rates in different working conditions brings obstacles to modeling and fault detection of the process. The invention provides a technical scheme aiming at the technical problem in the prior art.
Of course, the method of the present invention can also be applied to simple single-mode chemical processes. Preferably, in the present invention, the chemical process is a multi-modal chemical process, such as a pre-decarbonization process in an ammonia synthesis process or a synthesis unit of an organic compound-like intermediate. Preferably, the chemical process to be monitored is a pre-decarburization process of an ammonia synthesis process.
In the present invention, the training sample set X can be represented as:
X＝[X(1)X(2)… X(i)… X(S)]∈RN×d
wherein: n is the number of samples, S is the number of sampling rates,
Figure BDA0002029210180000031
representing a matrix of data at the ith sampling rate, having a number of samples NiAnd sample dimension di，
Figure BDA0002029210180000032
For training sample set XThe total dimension.
Preferably, the training sample set X is preprocessed so that all elements in each data sample in the training sample set X fluctuate around 0, where greater than 0 indicates higher than average level, less than 0 indicates lower than average level, and there is a linear correlation relationship with latent variables, so that the following static model (in m-mode) can be obtained:
Figure BDA0002029210180000033
wherein, therein
Figure BDA0002029210180000034
For training samples at the ith sampling rate, Wm(i) And mum(i) The divergence matrix and mean value corresponding to the sample at the ith sampling rate in the m mode ηm(i) η for the noise corresponding to the sample at the ith sampling rate in the m-modem(i) Obey a 0-mean Gaussian distribution, i.e. satisfy
Figure BDA0002029210180000035
βm(i) Is the variance of the noise at that sampling rate,
Figure BDA0002029210180000036
is di×diThe unit diagonal matrix of (1), wherein βm(i)＝σm(i)2，σm(i) Is the standard deviation of the noise at that sampling rate. t is tn,m∈RqThe implicit variable of the training sample under the m mode is subject to the standard normal distribution, namely t is satisfiedn,m～N(0,Iq) Q is the dimension of an implicit variable, IqIs an identity diagonal matrix of q × q.
Thus for the global model, the divergence matrix WmMean value of μmAnd noise ηmOf the variance matrix sigmamCan be expressed as:
Wm＝[Wm(1)；Wm(2)；...；Wm(i)；...；Wm(S)]
μm＝[μm(1)；μm(2)；...；μm(r)；...；μm(S)]
Figure BDA0002029210180000037
Wm(i) a divergence matrix corresponding to a training sample at the ith sampling rate in the mth mode; mu.sm(i) The mean value corresponding to the training sample under the ith sampling rate in the mth mode;
Figure BDA0002029210180000038
corresponding noise variance of the training sample under the mth mode;
for the nth training sample xnThe probability density function can be expressed as:
Figure BDA0002029210180000039
πmis the probability that the sample belongs to the mth mode; p (x)n| m) is the probability that the training sample occurs in the mth modality; t is tn,mThe implicit variable of the sample in the m mode is taken; p (x)n|tn,mM) is the data sample with respect to the hidden variable t in the mth modalityn,mThe conditional probability of (a); p (t)n,mIm) is the probability that the hidden variable of the data sample occurs in the m-th modality.
Thus, the model parameter set for mrmppc a (i.e., the mixed multisampling rate probabilistic principal component analysis model) is:
{πm,Wm(i),μm(i),βm(i)},(i＝1,2,...,S；m＝1,2,...,M)
πmthe probability that the current training sample belongs to the m mode is obtained; wm(i) A divergence matrix of the training sample under the ith sampling rate under the m mode; mu.sm(i) β is the mean value of the training sample at the ith sampling rate in the m modem(i) The noise variance of the training sample at the ith sampling rate in M modes is shown, and M is the total mode number.
Updating model parameters by adopting an expectation maximization algorithm during the constructed mixed multi-sampling rate probability principal component analysis model, and estimating the posterior probability of latent variables by using the current model parameters in the step E; in the M step, updating the parameters of the mixed multi-sampling rate model in a mode of a maximum likelihood function; and repeating the step E and the step M until reaching the model convergence condition.
First, model parameters { π are matchedm,Wm(i),μm(i),βm(i) Initializing (i ═ 1, 2.., S, and M ═ 1, 2.., M) randomly; at the same time, we define a certain training sample xnComprising SnA different sampling rate, then xnCan be expressed as
Figure BDA0002029210180000041
Wherein
Figure BDA0002029210180000042
Is a sequence number of a sampling rate and satisfies
Figure BDA0002029210180000043
Figure BDA0002029210180000044
Is as follows
Figure BDA0002029210180000045
Training samples at a sampling rate. Then it can be written accordingly:
Figure BDA0002029210180000046
is the sample xnCorresponding mean value in the m-th mode:
Figure BDA0002029210180000047
Figure BDA0002029210180000048
is in the m-th mode
Figure BDA0002029210180000049
Sample x at one sampling ratenA corresponding mean value;
Figure BDA00020292101800000410
is the sample x in the m-th modenThe corresponding divergence matrix in the mth mode:
Figure BDA00020292101800000411
Figure BDA00020292101800000412
is in the m-th mode
Figure BDA00020292101800000413
The sample x at a sampling ratenA corresponding divergence matrix;
Figure BDA0002029210180000051
is the sample xnCorresponding noise variance in the mth mode:
Figure BDA0002029210180000052
Figure BDA0002029210180000053
is as follows
Figure BDA0002029210180000054
Sample x at one sampling ratenVariance of the corresponding noise;
Figure BDA0002029210180000055
is di×diThe unit diagonal matrix of (2);
diag { } is a diagonal matrix
And then in the step E of model parameter estimation, obtaining an updated value of the model latent variable estimation according to the initial value of the current model parameter, wherein the main formula is as follows:
Figure BDA0002029210180000056
Figure BDA0002029210180000057
Figure BDA0002029210180000058
to simplify the above formula, we define:
Figure BDA0002029210180000059
<zn,m>for the training sample xn(ii) a posterior probability expectation belonging to the mth modality;<tn,m>for the training sample xnThe posterior probability expectation of the hidden variable in the m-th mode.
Figure BDA00020292101800000510
For the training sample xnA posterior probability covariance matrix of the hidden variables in the m-th mode.
Comparing the maximum likelihood value theta corresponding to the new model parameternewMaximum likelihood value Θ o corresponding to its original model parameterldIf | | | Θnew-Θold||2If yes, entering the fourth step, otherwise, continuing to iterate the EM algorithm, wherein the complete log-maximum likelihood estimation formula of the model is as follows for a threshold value of model convergence:
Figure BDA0002029210180000061
where Θ represents a maximum likelihood function value, const represents an arbitrary constant, and trace () represents a trace of a matrix.
In step M, obtaining model parameters { pi ] according to the result of step Em,Wm(i),μm(i),βm(i) The update values of (i ═ 1, 2.., S; M ═ 1, 2.., M) are as follows:
Figure BDA0002029210180000062
Figure BDA0002029210180000063
Figure BDA0002029210180000064
Figure BDA0002029210180000065
xn(i) is a sample xnA subvector consisting of variables at the ith sampling rate; Σ represents the sum over all acquired samples at that sampling rate; trace () represents the traces of the matrix.
In step (4), collecting new multiple sampling rate process monitoring sample x of chemical process on linenewThe monitoring sample contains SnewA different sampling rate, then xnewCan be expressed as:
Figure BDA0002029210180000066
wherein
Figure BDA0002029210180000067
A sequence number representing a sampling rate and satisfying
Figure BDA0002029210180000068
Figure BDA0002029210180000069
Is as follows
Figure BDA00020292101800000610
Monitoring samples at individual sampling rates.
In the step (5) and the step (6):
(5-1) solving monitoring sample x based on constructed mixed multi-sampling rate probability principal component analysis modelnewT in the m-th mode2Statistics:
Figure BDA0002029210180000071
and SPE statistics:
Figure BDA0002029210180000072
obtaining M T2Statistics:
Figure BDA0002029210180000073
and SPE statistics:
Figure BDA0002029210180000074
after the monitoring sample is obtained, the same pretreatment and standardization are carried out on the monitoring sample, and the data sample x can be obtained according to the constructed mixed multi-sampling rate probability principal component analysis modelnewThe corresponding score belongs to the mean vector under M modes
Figure BDA0002029210180000075
Divergence matrix
Figure BDA0002029210180000076
Sum noise covariance matrix
Figure BDA0002029210180000077
Respectively as follows:
Figure BDA0002029210180000078
Figure BDA0002029210180000079
Figure BDA00020292101800000710
wherein
Figure BDA00020292101800000711
For monitoring a sample xnewIn the m-th mode
Figure BDA00020292101800000712
A divergence matrix at one sampling rate.
Figure BDA00020292101800000713
For monitoring a sample xnewIn the m-th mode
Figure BDA00020292101800000714
Mean vector at each sampling rate.
Figure BDA00020292101800000715
For monitoring a sample xnewIn the m-th mode
Figure BDA00020292101800000716
The noise variance at each sampling rate.
Firstly, the monitoring sample x is obtainednewPosterior probability expectation belonging to the m-th mode<znew,m>Namely:
Figure BDA00020292101800000717
the monitoring sample x can be obtained under the model of the m-th modenewIs a hidden variable tnew,mExpected value of posterior distribution of (1):
Figure BDA00020292101800000718
where we define for simplicity the formula:
Figure BDA00020292101800000719
the monitoring sample xnewIn the m-th modeIs/are as follows
Figure BDA00020292101800000720
The statistics may be obtained by:
Figure BDA00020292101800000721
further, we can obtain a monitoring sample xnewConditional probability distribution in the m-th mode
Figure BDA00020292101800000722
Wherein:
Figure BDA00020292101800000723
then the observation sample xnewThe residual at the mth mode is:
Figure BDA0002029210180000081
the sample xnewIn the m-th mode
Figure BDA0002029210180000082
The statistics may be expressed as:
Figure BDA0002029210180000083
thus obtaining M total T2Statistics and SPE statistics.
(5-2) calculating the monitor sample xnewIn the m-th mode
Figure BDA0002029210180000084
And
Figure BDA0002029210180000085
counting the probability of failure;
for the obtained M T2Statistics and SPE statistics, calculating monitorSample xnewIn the m-th mode
Figure BDA0002029210180000086
And
Figure BDA0002029210180000087
the probability of the fault occurrence of the statistic is respectively as follows:
Figure BDA0002029210180000088
Figure BDA0002029210180000089
wherein
Figure BDA00020292101800000810
And
Figure BDA00020292101800000811
is assumed to be:
Figure BDA00020292101800000812
Figure BDA00020292101800000813
where (1- α) is the confidence level, it may be set to 0.99, i.e., α ═ 0.01.
Sample xnewIn that
Figure BDA00020292101800000814
And
Figure BDA00020292101800000815
the conditional probabilities of normal samples (N) and fault samples (F) in the statistics are:
Figure BDA00020292101800000816
Figure BDA00020292101800000817
Figure BDA00020292101800000818
Figure BDA00020292101800000819
wherein
Figure BDA00020292101800000820
And
Figure BDA00020292101800000821
respectively as statistical confidence limits in the mth mode
Figure BDA00020292101800000822
The distribution, g and h, can be approximated by:
Figure BDA0002029210180000091
Figure BDA0002029210180000092
wherein
Figure BDA0002029210180000093
And
Figure BDA0002029210180000094
statistics calculated for samples belonging to the mth modality in the modeling data, respectively. mean represents mean and Var represents variance.
(5-3) binding-monitoring sample xnewPosterior probability of each mode to obtain fused
Figure BDA0002029210180000095
And SPEnewStatistics are obtained.
Combining on-line samples xnewThe posterior probability of each mode is<znew,m>Then after fusion
Figure BDA0002029210180000096
And SPEnewThe statistics are:
Figure BDA0002029210180000097
Figure BDA0002029210180000098
in step (6), the online sample x is samplednewStatistic of (2)
Figure BDA0002029210180000099
And SPEnewAre compared to the value of the confidence level α to determine whether the sample is a fault.
The invention relates to a fault detection method based on a mixed multi-sampling rate probability principal component analysis model, which comprises the steps of establishing the multi-sampling rate probability principal component analysis model under each mode, fusing a plurality of sub-mode models by a mixed model method, extracting mode information and variable autocorrelation relation of a process, diagnosing faults by using the autocorrelation relation and providing a corresponding online fault detection statistic construction method. Compared with the traditional multi-mode modeling method, the MrMPPCA model (namely a mixed multi-sampling rate probability principal component analysis model) provided by the invention can simultaneously consider the multi-sampling rate data characteristics and the multi-mode process characteristics, so that the applicability and the fault detection precision on the fault detection problem of the multi-mode multi-sampling rate process are improved.
Detailed Description
The invention is further explained by taking a pre-decarbonization unit in the synthetic ammonia process as an example:
a fault detection method based on a mixed multi-sampling rate probability principal component analysis model is disclosed. Aiming at the problem of multi-mode working condition fault detection of a pre-decarburization unit in a synthetic ammonia process, the method firstly utilizes a discrete control system to collect multi-sampling rate data under a plurality of normal working modes and establishes a mixed multi-sampling rate probability principal component analysis model. The model structure is estimated by an expectation-maximization algorithm. On the basis, sampling is carried out on the data of the on-line pre-decarburization unit to obtain a multi-sampling rate test sample, latent variables of the test sample are estimated by utilizing the existing model structure, and statistical monitoring indexes at the moment are calculated to realize the fault diagnosis result of the pre-decarburization unit.
The multi-mode process means that for a certain chemical unit, multiple products may respectively correspond to one working condition according to the requirements of different products, and the variable number and the sampling rate under each working condition may change; taking the pre-decarburization unit in the synthetic ammonia process as an example, multiple working conditions may exist according to different component contents of production raw materials or different requirements of final products, and the working conditions may have differences in parameters such as stable working condition operating points, reaction internal parameters, final product quality indexes and the like. Meanwhile, the difference of sampling rates in different working conditions brings obstacles to modeling and fault detection of the process.
The invention relates to a multi-mode multi-sampling rate fault detection method based on a mixed multi-sampling rate probability principal component analysis model and a synthetic ammonia process pre-decarburization process, which comprises the following steps:
the first step is as follows: collecting data of different sampling rates under M (M is more than or equal to 0) modes in normal operation in the pre-decarburization process of the synthetic ammonia process by using a distributed control system, and forming a training sample set X for modeling to be expressed as:
X＝[X(1)X(2)... X(i)... X(S)]∈RN×d
wherein: n is the number of samples, S is the number of sampling rates,
Figure BDA0002029210180000101
representing a matrix of data at the ith sampling rate, having a number of samples NiAnd sample dimension di，
Figure BDA0002029210180000102
Is the total dimension of the training sample set X.
The second step is that: preprocessing and normalizing the data set X, i.e. for a training sample, averaging all elements in the sample, and then subtracting the average value from each element, so that the respective normalized variable value (or element) fluctuates around 0, more than 0 indicates higher than average level, less than 0 indicates lower than average level, and there is a linear correlation with latent variables, the following static model can be obtained (m.di.m. 1,2 … M in the mth modality):
Figure BDA0002029210180000103
wherein, therein
Figure BDA0002029210180000104
For training samples at the ith sampling rate, Wm(i) And mum(i) The divergence matrix and mean value corresponding to the sample at the ith sampling rate in the mth mode, ηm(i) η for the noise corresponding to the sample at the ith sampling rate in the mth modem(i) Obey a 0-mean Gaussian distribution, i.e. satisfy
Figure BDA0002029210180000111
βm(i) Is the variance of the noise at that sampling rate,
Figure BDA0002029210180000112
is di×diThe unit diagonal matrix of (1), wherein βm(i)＝σm(i)2，σm(i) Is the standard deviation of the noise at that sampling rate. t is tn,m∈RqThe implicit variable of the training sample under the m mode is subject to the standard normal distribution, namely t is satisfiedn,m～N(0,Iq) Q is the dimension of an implicit variable, IqIs the unit diagonal matrix of q × q for the global model, the divergence matrix W is thereforemMean value of μmAnd noise ηmOf the variance matrix sigmamCan be expressed as:
Wm＝[Wm(1)；Wm(2)；...；Wm(i)；...；Wm(S)]
μm＝[μm(1)；μm(2)；...；μm(i)；...；μm(S)]
Figure BDA0002029210180000113
Wm(i) a divergence matrix corresponding to a training sample at the ith sampling rate in the mth mode; mu.sm(i) The mean value corresponding to the training sample under the ith sampling rate in the mth mode;
Figure BDA0002029210180000114
corresponding noise variance of the training sample under the mth mode;
for the nth training sample xnThe probability density function can be expressed as:
Figure BDA0002029210180000115
πmthe probability that the training sample belongs to the m mode is taken as the probability; p (x)n| m) is the probability that the training sample occurs in the mth modality; t is tn,mThe hidden variable of the training sample under the m mode; p (x)n|tn,mM) is the data sample with respect to the hidden variable t in the mth modalityn,mThe conditional probability of (a); p (t)n,mIm) is the probability that the hidden variable of the data sample occurs in the m-th modality.
Thus, the set of model parameters for MrMPPCA (i.e., the mixed multisampling rate probabilistic principal component analysis model) is
{πm,Wm(i),μm(i),βm(i)},(i＝1,2,...,S；m＝1,2,...,M)
Wm(i) A divergence matrix corresponding to a training sample at the ith sampling rate in the mth mode; mu.sm(i) Is the ith sampling rate in the mth modeMean of training samples βm(i) And the noise variance corresponding to the training sample at the ith sampling rate in the mth mode.
The third step: updating the model parameters by using an Expectation Maximization (EM) algorithm, and estimating the posterior probability of the latent variable by using the current model parameters in the step E; in M, updating the parameters of the mixed multi-sampling rate model by means of a maximum likelihood function. And finally, repeating the step E and the step M until reaching the model convergence condition.
First, model parameters { π are matchedm,Wm(i),μm(i),βm(i) Initializing (i ═ 1, 2.., S, and M ═ 1, 2.., M) randomly; at the same time, we define a certain training sample xnComprising SnA different sampling rate, then xnCan be expressed as
Figure BDA0002029210180000121
Wherein
Figure BDA0002029210180000122
Is a sequence number of the sampling rate and satisfies
Figure BDA0002029210180000123
Figure BDA0002029210180000124
Is as follows
Figure BDA0002029210180000125
Training samples at a sampling rate. Then it can be written accordingly:
Figure BDA0002029210180000126
is the sample xnCorresponding mean value at the mth modality:
Figure BDA0002029210180000127
Figure BDA0002029210180000128
is the m < th > modality
Figure BDA0002029210180000129
Sample x at one sampling ratenA corresponding mean value;
Figure BDA00020292101800001210
is the sample xnIn the m-th mode
Figure BDA00020292101800001211
Corresponding divergence matrix at each sampling rate:
Figure BDA00020292101800001212
Figure BDA00020292101800001213
is the sample xnIn the m-th mode
Figure BDA00020292101800001214
The sample x at a sampling ratenA corresponding divergence matrix;
Figure BDA00020292101800001215
is the sample xnCorresponding noise variance in the mth mode:
Figure BDA00020292101800001216
Figure BDA00020292101800001217
for the m-th mode
Figure BDA00020292101800001218
Sample x at one sampling ratenVariance of the corresponding noise;
Figure BDA00020292101800001219
is di×diThe unit diagonal matrix of (2);
diag { } is a diagonal matrix.
And then in the step E of model parameter estimation, obtaining an updated value of the model latent variable estimation according to the initial value of the current model parameter, wherein the main formula is as follows:
Figure BDA00020292101800001220
Figure BDA00020292101800001221
Figure BDA00020292101800001222
to simplify the above formula, we define:
Figure BDA0002029210180000131
<zn,m>for the training sample xn(ii) a posterior probability expectation belonging to the mth modality;<tn,m>for the training sample xnThe posterior probability expectation of the hidden variable in the m-th mode.
Figure BDA0002029210180000132
For the training sample xnA posterior probability covariance matrix of the hidden variables in the m-th mode.
Comparing the maximum likelihood value theta corresponding to the new model parameternewMaximum likelihood value theta corresponding to original model parameteroldIf | | | Θnew-Θold||2If yes, entering the fourth step, otherwise, continuing to iterate the EM algorithm, wherein the complete log-maximum likelihood estimation formula of the model is as follows for a threshold value of model convergence:
Figure BDA0002029210180000133
where Θ represents a maximum likelihood function value, const represents an arbitrary constant, and trace () represents a trace of a matrix.
In step M, obtaining model parameters { pi ] according to the result of step Em,Wm(i),μm(i),βm(i) The update values of (i ═ 1, 2.., S; M ═ 1, 2.., M) are as follows:
Figure BDA0002029210180000134
Figure BDA0002029210180000135
Figure BDA0002029210180000136
Figure BDA0002029210180000141
xn(i) is a sample xnA subvector consisting of variables at the ith sampling rate; Σ represents the sum over all acquired samples at that sampling rate; trace () represents the traces of the matrix.
The fourth step: on-line collection of multiple sampling rate process monitoring samples x for a new synthetic ammonia pre-decarbonization processnewThe monitoring sample contains SnewA different sampling rate, then xnewCan be expressed as:
Figure BDA0002029210180000142
wherein
Figure BDA0002029210180000143
A sequence number representing a sampling rate and satisfying
Figure BDA0002029210180000144
Figure BDA0002029210180000145
Is as follows
Figure BDA0002029210180000146
Monitoring samples at individual sampling rates.
For the monitoring sample xnewThe same pretreatment and normalization as in the second step were performed. The monitoring sample x can be obtained according to the constructed mixed multi-sampling rate probability principal component analysis modelnewThe corresponding score belongs to the mean vector under M modes
Figure BDA0002029210180000147
Divergence matrix
Figure BDA0002029210180000148
Sum noise covariance matrix
Figure BDA0002029210180000149
Respectively as follows:
Figure BDA00020292101800001410
Figure BDA00020292101800001411
Figure BDA00020292101800001412
wherein
Figure BDA00020292101800001413
For monitoring a sample xnewIn the m-th mode
Figure BDA00020292101800001414
A divergence matrix at one sampling rate.
Figure BDA00020292101800001415
For monitoring a sample xnewIn the m-th mode
Figure BDA00020292101800001416
Mean vector at each sampling rate.
Figure BDA00020292101800001417
For monitoring a sample xnewIn the m-th mode
Figure BDA00020292101800001418
The noise variance at each sampling rate.
Firstly, the monitoring sample x is obtainednewPosterior probability expectation z belonging to the m-th modenew,m>Namely:
Figure BDA00020292101800001419
the monitoring sample x can be obtained under the model of the m-th modenewIs a hidden variable tnew,mExpected value of posterior distribution of (1):
Figure BDA0002029210180000151
where we define for simplicity the formula:
Figure BDA0002029210180000152
the monitoring sample xnewIn the m-th mode
Figure BDA0002029210180000153
The statistics may be obtained by:
Figure BDA0002029210180000154
further, we can obtain a monitoring sample xnewConditional probability distribution in the m-th mode
Figure BDA0002029210180000155
Wherein:
Figure BDA0002029210180000156
then the observation sample xnewThe residual at the mth mode is:
Figure BDA0002029210180000157
the sample xnewIn the m-th mode
Figure BDA0002029210180000158
The statistics may be expressed as:
Figure BDA0002029210180000159
thus obtaining M total T2Statistics and SPE statistics.
Then the observation sample x is calculatednewIn the m-th mode
Figure BDA00020292101800001510
And
Figure BDA00020292101800001511
the probability of the fault occurrence of the statistic is respectively as follows:
Figure BDA00020292101800001512
Figure BDA00020292101800001513
wherein
Figure BDA00020292101800001514
And
Figure BDA00020292101800001515
is assumed to be:
Figure BDA00020292101800001516
Figure BDA00020292101800001517
where (1- α) is the confidence level, it may be set to 0.99, i.e., α ═ 0.01.
Sample xnewIn that
Figure BDA00020292101800001518
And
Figure BDA00020292101800001519
the conditional probabilities of normal samples (N) and fault samples (F) in the statistics are:
Figure BDA0002029210180000161
Figure BDA0002029210180000162
Figure BDA0002029210180000163
Figure BDA0002029210180000164
wherein
Figure BDA0002029210180000165
And
Figure BDA0002029210180000166
respectively as statistical confidence limits in the mth mode
Figure BDA0002029210180000167
The distribution, g and h, can be approximated by:
Figure BDA0002029210180000168
Figure BDA0002029210180000169
wherein
Figure BDA00020292101800001610
And
Figure BDA00020292101800001611
statistics calculated for samples belonging to the mth modality in the modeling data, respectively. mean represents mean and Var represents variance.
Combining on-line samples xnewThe posterior probability of each mode is<znew,m>Then after fusion
Figure BDA00020292101800001612
And SPEnewThe statistics are:
Figure BDA00020292101800001613
Figure BDA00020292101800001614
and a sixth step: to sample x onlinenewStatistic of (2)
Figure BDA00020292101800001615
And SPEnewMay be compared to the confidence level α to determine whether the sample is a fault.

Claims (5)

1. A fault detection method based on a mixed multi-sampling rate probability principal component analysis model is characterized by comprising the following steps:
(1) collecting data of different sampling rates under M modes of normal operation of a chemical process to be monitored, wherein M is more than or equal to 1, and forming a training sample set X for modeling;
(2) preprocessing a training sample set X;
(3) constructing a mixed multi-sampling rate probability principal component analysis model by utilizing the preprocessed training sample set X;
(4) on-line collection of monitoring samples x of a multi-sampling-rate process of a chemical process to be monitorednew；
(5) Based on the constructed mixed multi-sampling rate probability principal component analysis model, calculating monitoring statistics of monitoring samples;
(6) judging whether the fault occurs according to the monitoring statistics;
preprocessing a training sample set X to enable all elements in each data sample in the training sample set X to fluctuate around 0;
a linear correlation relationship exists between the preprocessed training sample set X and latent variables;
when constructing the mixed multi-sampling rate probability principal component analysis model, the model parameter set is as follows:
{πm,Wm(i),μm(i),βm(i)},(i＝1,2,...,S；m＝1,2,...,M)
πmthe probability that the current training sample belongs to the m mode is obtained; wm(i) A divergence matrix of the training sample under the ith sampling rate under the m mode; mu.sm(i) β is the mean value of the training sample at the ith sampling rate in the m modem(i) Noise variance of a training sample under the ith sampling rate under M modes is obtained, and M is the total mode number;
in the step (5) and the step (6):
(5-1) solving monitoring sample x based on constructed mixed multi-sampling rate probability principal component analysis modelnewT in the m-th mode2Statistics:
Figure FDA0002549832500000011
and SPE statistics:
Figure FDA0002549832500000012
obtaining M T2Statistics:
Figure FDA0002549832500000013
and SPE statistics:
Figure FDA0002549832500000014
(5-2) calculating the monitor sample xnewIn the m-th mode
Figure FDA0002549832500000015
And
Figure FDA0002549832500000016
counting the probability of failure;
(5-3) binding-monitoring sample xnewPosterior probability of each mode to obtain fused
Figure FDA0002549832500000021
And SPEnewStatistics are obtained.
2. The fault detection method based on the mixed multisampling rate probabilistic principal component analysis model according to claim 1, wherein in step (1) or step (4), a discrete control system is used to collect data.
3. The method of claim 1, wherein the chemical process to be monitored is a pre-decarbonization process of a synthetic ammonia process.
4. The fault detection method based on the hybrid multisampling rate probabilistic principal component analysis model according to claim 1, wherein the model parameters are updated by adopting an expectation-maximization algorithm when the hybrid multisampling rate probabilistic principal component analysis model is constructed, and the posterior probability of latent variables is estimated by using the current model parameters in step E of the expectation-maximization algorithm; in M steps in the expectation-maximization algorithm, updating the parameters of the mixed multi-sampling rate model in a mode of a maximization likelihood function; and repeating the step E and the step M until reaching the model convergence condition.
5. The fault detection method based on the mixed multisampling rate probabilistic principal component analysis model as claimed in claim 1, wherein in the step (6), the monitoring sample x isnewFused statistics of
Figure FDA0002549832500000022
And SPEnewAre compared to the value of the confidence level α to determine whether the sample is a fault.
CN201910304064.3A 2019-04-16 2019-04-16 Fault detection method based on mixed multi-sampling rate probability principal component analysis model Active CN109917777B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910304064.3A CN109917777B (en) 2019-04-16 2019-04-16 Fault detection method based on mixed multi-sampling rate probability principal component analysis model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910304064.3A CN109917777B (en) 2019-04-16 2019-04-16 Fault detection method based on mixed multi-sampling rate probability principal component analysis model

Publications (2)

Publication Number Publication Date
CN109917777A CN109917777A (en) 2019-06-21
CN109917777B true CN109917777B (en) 2020-08-25

Family

ID=66977421

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910304064.3A Active CN109917777B (en) 2019-04-16 2019-04-16 Fault detection method based on mixed multi-sampling rate probability principal component analysis model

Country Status (1)

Country Link
CN (1) CN109917777B (en)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110988674A (en) * 2019-11-19 2020-04-10 中南大学 Health state monitoring method and system of permanent magnet synchronous motor and mobile terminal
CN111913446B (en) * 2020-06-13 2022-03-18 宁波大学 Fault detection method suitable for multi-sampling-rate chemical process
CN112286169B (en) * 2020-10-20 2022-02-01 浙江钱江机器人有限公司 Industrial robot fault detection method
CN112489841A (en) * 2020-11-23 2021-03-12 浙江大学 Water level fault-tolerant control method for steam generator of nuclear power unit
CN114326486B (en) * 2021-12-13 2024-01-19 山东科技大学 Process monitoring method based on probability slow feature analysis and elastic weight consolidation
CN114358332A (en) * 2021-12-23 2022-04-15 国网江苏省电力有限公司南通供电分公司 Power grid upgrading work order feature portrait and analysis method based on probability

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777830A (en) * 2015-04-01 2015-07-15 浙江大学 Multi-work-condition process monitoring method based on KPCA (kernel principal component analysis) mixture model
CN105469101A (en) * 2015-12-31 2016-04-06 北京工业大学 Mixed two-dimensional probabilistic principal component analysis method
CN107292323A (en) * 2016-03-31 2017-10-24 日本电气株式会社 Method and apparatus for training mixed model
CN108549908A (en) * 2018-04-13 2018-09-18 浙江科技学院 Chemical process fault detection method based on more sampled probability core principle component models
CN109085805A (en) * 2018-07-24 2018-12-25 浙江科技学院 A kind of industrial process fault detection method based on multi-sampling rate Factor Analysis Model

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6643569B2 (en) * 2001-03-30 2003-11-04 The Regents Of The University Of Michigan Method and system for detecting a failure or performance degradation in a dynamic system such as a flight vehicle

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777830A (en) * 2015-04-01 2015-07-15 浙江大学 Multi-work-condition process monitoring method based on KPCA (kernel principal component analysis) mixture model
CN105469101A (en) * 2015-12-31 2016-04-06 北京工业大学 Mixed two-dimensional probabilistic principal component analysis method
CN107292323A (en) * 2016-03-31 2017-10-24 日本电气株式会社 Method and apparatus for training mixed model
CN108549908A (en) * 2018-04-13 2018-09-18 浙江科技学院 Chemical process fault detection method based on more sampled probability core principle component models
CN109085805A (en) * 2018-07-24 2018-12-25 浙江科技学院 A kind of industrial process fault detection method based on multi-sampling rate Factor Analysis Model

Also Published As

Publication number Publication date
CN109917777A (en) 2019-06-21

Similar Documents

Publication Publication Date Title
CN109917777B (en) Fault detection method based on mixed multi-sampling rate probability principal component analysis model
CN112202736B (en) Communication network anomaly classification method based on statistical learning and deep learning
CN107153874B (en) Water quality prediction method and system
CN109000940B (en) Abnormal axle temperature diagnosis method and system for rolling stock
CN102789545B (en) Based on the Forecasting Methodology of the turbine engine residual life of degradation model coupling
CN108960303B (en) Unmanned aerial vehicle flight data anomaly detection method based on LSTM
US8725456B1 (en) Decomposition technique for remaining useful life prediction
Coble et al. Applying the general path model to estimation of remaining useful life
CN109085805B (en) Industrial process fault detection method based on multi-sampling-rate factor analysis model
CN108549908B (en) Chemical process fault detection method based on multi-sampling probability kernel principal component model
CN111222290A (en) Large-scale equipment residual service life prediction method based on multi-parameter feature fusion
CN102208028A (en) Fault predicting and diagnosing method suitable for dynamic complex system
CN109739214B (en) Method for detecting intermittent faults in industrial process
CN109298633A (en) Chemical production process fault monitoring method based on adaptive piecemeal Non-negative Matrix Factorization
Dong A tutorial on nonlinear time-series data mining in engineering asset health and reliability prediction: concepts, models, and algorithms
CN112800616A (en) Equipment residual life self-adaptive prediction method based on proportional acceleration degradation modeling
Wang et al. A hidden semi-markov model with duration-dependent state transition probabilities for prognostics
CN114943179A (en) Reliability evaluation and residual life prediction method based on multi-source degradation data fusion
JP3718765B2 (en) Plant diagnostic equipment
CN109325065B (en) Multi-sampling-rate soft measurement method based on dynamic hidden variable model
CN110084301B (en) Hidden Markov model-based multi-working-condition process working condition identification method
Grebenyuk Monitoring and identification of structural shifts in processes with a unit root
CN111126477A (en) Learning and reasoning method of hybrid Bayesian network
CN112016208B (en) Hidden fault diagnosis method and system considering disturbance
CN113269327A (en) Flow anomaly prediction method based on machine learning

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