Fault Detection Method Based on Dynamic Non-Stationary Projection Structure

Abstract

The present invention discloses a fault detection method based on a dynamic non-stationary projection structure. The method consists of forming a training sample set from the measured variables collected during the normal operation of a chemical process and performing normalization processing; combining the normalized training sample set, using the expectation-maximization algorithm and the forward-backward algorithm to construct and train a dynamic non-stationary projection structure; solving the characteristic statistic through the dynamic non-stationary projection structure and determining the control limit of the characteristic statistic; collecting the measured variables during the chemical process to be measured online and performing normalization processing; combining the normalized variables to be measured and using the dynamic non-stationary projection structure to solve the characteristic statistic during the chemical process to be measured, and then judging whether there is a fault during the chemical process to be measured according to these characteristic statistics. The present invention provides a complete monitoring framework for industrial process fault detection, can effectively extract the dynamic relationship of the measured variables, and is more suitable for monitoring non-stationary industrial processes with dynamic characteristics.

Description

Technical Field

[0001] The present invention relates to a fault detection method in the field of industrial process monitoring, and more particularly to a fault detection method based on a dynamic non-stationary projection structure. Background Art

[0002] In order to ensure the safe and stable operation of industrial processes and improve production efficiency, it is necessary to detect faults existing in the process operation at an early stage. The advent of the big data era has promoted the widespread application of data-driven process monitoring technologies. Among them, the multivariate statistical variable analysis method, as an important branch of data-driven methods, provides a strong guarantee for enhancing process reliability and product quality. In most multivariate analysis methods, the stationarity of data is a necessary assumption. In other words, the non-stationary characteristics of the process are not considered. However, non-stationary characteristics not only exist widely but also are important factors dominating industrial operation processes, such as the blast furnace ironmaking process and the penicillin fermentation process. The non-stationary characteristics of industrial processes are mainly manifested as adjustments in production plans, switches in operation stages, equipment wear, and drifts in process parameters. The measurement variables collected from industrial processes dominated by non-stationary characteristics inevitably have non-stationary trends, resulting in unreliable monitoring results of traditional multivariate analysis methods.

[0003] To solve the problems brought by non-stationary characteristics in the field of process monitoring, scholars have proposed a series of methods. First is the adaptive modeling method, whose main idea is to adapt to the switching of process conditions by continuously updating model parameters. Its advantage is that the modeling is simple and effective, but the problem is that frequent updates may cause the monitoring model to adapt to some early fault trends, resulting in missed alarms. Then there is the cointegration analysis method. By extracting the cointegration relationship between non-stationary variables, stationary cointegration variables are obtained, and then the cointegration variables are monitored. Essentially, the method based on cointegration analysis monitors the cointegration relationship between non-stationary variables, so it requires all variables to have the same integration order. When there is no cointegration relationship between variables, cointegration analysis cannot be implemented. Finally, the method based on subspace decomposition attempts to separate the steady-state subspace from the original data space and study it, such as stationary subspace analysis.

[0004] Although the above non-stationary methods can be applied to non-stationary industrial process monitoring. However, most of the existing methods do not consider and extract the dynamic characteristics of non-stationary process data. For actual industrial processes, feedback control is widespread and makes the measurement variables have autocorrelation. In order to model the dynamics of process data, the dynamic stationary subspace analysis method was proposed and successfully applied to non-stationary process monitoring. However, due to factors such as unmeasurable disturbances and sensor errors, the measurement variables often have uncertainties. To solve this problem, a feasible method is to establish a dynamic non-stationary probability model in a generative framework. Summary of the Invention

[0005] To overcome the defects of the prior art, the present invention provides a fault detection method based on a dynamic non-stationary projection structure. By adopting a specific model structure to model the uncertainty and dynamics of an industrial process, while extracting the stationary features and non-stationary features of the industrial process, and then using the extracted features to develop three statistics, namely, residual statistics, stationary feature statistics, and non-stationary feature statistics, for detecting faults in non-stationary industrial processes and realizing the monitoring of the operating state of non-stationary industrial processes.

[0006] To achieve the object of the present invention, the technical solution adopted by the present invention is as follows:

[0007] A fault detection method based on a dynamic non-stationary projection structure, the method comprising the following steps:

[0008] Step 1: A training sample set is composed of measurement variables X collected during a chemical operation process without faults and related to the faults.

[0009] Step 2: Normalize the measurement variables X of all training samples in the training sample set in Step 1 to obtain a normalized training sample set X * ;

[0010] Step 3: Based on the normalized training sample set X * , construct a dynamic non-stationary projection structure and train the projection structure.

[0011] Step 4: Combine the normalized training sample set to determine the control limits of the respective multiple characteristic statistics of the projection structure by solving the multiple characteristic statistics of the projection structure.

[0012] Step 5: Online collect the measured variables x(c) during the chemical operation process to be measured and normalize them to obtain a normalized set of measured variables to be measured.

[0013] Step 6: Combine the normalized set of measured variables to be measured and use the dynamic non-stationary projection structure to determine the respective characteristic statistics during the chemical operation process to be measured, and then determine whether there are faults during the chemical operation process to be measured according to the magnitude relationship between the respective characteristic statistics during the chemical operation process to be measured and the control limits of the multiple characteristic statistics of the dynamic non-stationary projection structure.

[0014] Preferably, the training sample set X * is expressed as: X * = [x * (1), x * (2), x * (t)...x * (N)], t ∈ [1, N], where x* (t) represents the measurement variable of the t-th training sample after normalization, and N represents the number of training samples in the training sample set; each sample corresponds to a moment.

[0015] Preferably, the method for constructing the dynamic non-stationary projection structure is as follows:

[0016]

[0017] Among them, B ∈ R m×p represents a linear superposition matrix, and B s ∈ R m×a is composed of the first a columns of B, and B n ∈ R m×(p-a) is composed of the last (p - a) columns of B, s s (t) ∈ R a×1 represents the stationary feature of the t-th training sample, and s n (t) ∈ R (p-a)×1 represents the non-stationary feature of the t-th training sample, and represents the full feature matrix of the t-th training sample, where R represents the set of real numbers, m represents the dimension of the measurement variable, p is the dimension of the full feature s(t) of the t-th training sample, and a is the dimension of the stationary feature s s (t) of the t-th training sample; e(t) represents the noise of the t-th training sample after normalization.

[0018] Preferably, the specific method for training the projection structure is as follows: First, input the training sample set into the projection structure, and calculate the posterior hidden state and posterior joint hidden state of the hidden Markov chain through the forward-backward algorithm and in combination with the current model parameters where A represents the state transition matrix, B represents the emission matrix, and π i represents the probability that the hidden state at the initial moment is i, and μ s represents the mean of the stationary features of the training samples, and Σ s represents the covariance of the stationary features of the training samples, represents the mean of the non-stationary features of the training samples when the hidden state is i, represents the covariance of the non-stationary features of the training samples when the hidden state is i, and σ 2 represents the noise intensity of the measurement variable; then calculate the first moment of the local full feature of the training samples and the second moment of the local full feature; then, in combination with the posterior hidden state, posterior hidden state transition variable, first moment of the local full feature, and second moment of the local full feature, use the likelihood function to perform iterative calculations repeatedly until the likelihood function converges, and update the current model parameters to obtain the optimal parameter set of the dynamic non-stationary projection structure Complete the training of the dynamic non-stationary projection structure.

[0019] Preferably, the 8 parameters in the optimal parameter set correspond one-to-one with the 8 parameters in the model parameter Θ, and are specifically determined by the following formula:

[0020]

[0021]

[0022]

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] where γ i (t) represents the posterior hidden state distribution at time t, and ξ ij (t) represents the posterior joint hidden state distribution at time t, x * (t) represents the sample at time t, i.e., the t-th sample, <s i (t)> and respectively represent the first moment and the second moment of the local full feature s i (t) of the t-th training sample. <·> represents the operation of taking the expectation, and s i (t) represents the local full feature of the t-th training sample when the hidden state at time t is i; I represents the number of hidden states, and W s ∈R p×a is the first a columns of the identity matrix I p , and W n ∈R p×(p-a) is the last (p - a) columns of the identity matrix I p . Tr(*) represents the trace of the matrix.

[0029] Preferably, the multiple feature statistics of the projection structure include stationary feature statistics non-stationary feature statistics and the residual statistic SPE; the control limits of the multiple feature statistics of the projection structure respectively include the control limit of the stationary feature statistics the control limit of the non-stationary feature statistics and the control limit of the residual statistic

[0030] The calculation formula of the residual statistic SPE is as follows:

[0031] SPE = [SPE(1), SPE(2), SPE(t)... SPE(N)], t ∈ [1, N]

[0032] Among them, SPE(t) represents the residual of the t-th training sample, and is calculated by the following formula:

[0033]

[0034] In the formula: e(t) represents the residual of the t-th training sample;

[0035] The stationary feature statistic The calculation formula is:

[0036]

[0037] Among them, represents the stationary feature statistic of the t-th training sample, and is calculated by the following formula:

[0038]

[0039] In the formula: represents the posterior probability that the hidden state corresponding to the t-th training sample is i, represents the new local full feature of the t-th training sample when the hidden state at time t is i The Mahalanobis distance between the corresponding stationary feature and its own first moment;

[0040] The non-stationary feature statistic The calculation formula is:

[0041]

[0042] Among them, represents the non-stationary feature statistic of the t-th training sample, and the specific calculation formula is as follows:

[0043]

[0044] In the formula: represents the new local full feature of the t-th training sample when the hidden state at time t is i The corresponding non-stationary feature The Mahalanobis distance between and its own local mean.

[0045] Preferably, the control limits of the multiple feature statistics of the projection structure are determined by combining the confidence level and using the kernel density estimation method.

[0046] Preferably, each characteristic statistic in the chemical process to be measured includes a stationary characteristic statistic to be measured a non-stationary characteristic statistic to be measured a residual statistic SPE(c) to be measured, and its calculation method is the same as that of the stationary characteristic statistic in step 4 the non-stationary characteristic statistic and the residual statistic SPE is consistent

[0047] Preferably, determining whether there is a fault in the chemical process to be measured in step 6 specifically includes

[0048] If each characteristic statistic in the chemical process to be measured is respectively less than the control limit of the corresponding characteristic statistic determined in step 4, that is, the stationary characteristic statistic in the chemical process to be measured is less than the control limit of the stationary characteristic statistic and the non-stationary characteristic statistic in the chemical process to be measured is less than the control limit of the non-stationary characteristic statistic and the residual statistic SPE(c) in the chemical process to be measured is less than the control limit of the residual statistic then there is no fault in the chemical process to be measured; otherwise, there is a fault in the chemical process to be measured, thereby completing the fault detection of the chemical process to be measured

[0049] The beneficial effects of the present invention are

[0050] This method models the uncertainty and dynamic characteristics of process data, projects the data into a stationary subspace and a non-stationary subspace, and establishes specific monitoring statistics in all subspaces, providing a complete monitoring scheme for industrial process fault detection. By introducing a hidden Markov model to model the non-stationary trend, the present invention can not only extract the current non-stationary characteristics, but also predict the changes of non-stationary characteristics, thereby realizing the dynamic characteristic monitoring of non-stationary processes. The present invention can effectively extract the dynamic information of non-stationary processes, so it is more suitable for monitoring non-stationary industrial processes with dynamic characteristics BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 is a flowchart of a fault detection method based on a dynamic non-stationary projection structure according to the present invention

[0052] Figure 2 is a schematic diagram of the penicillin fermentation process in an embodiment of the present invention

[0053] Figure 3 is a monitoring result diagram of stationary subspace analysis-Mahalanobis distance in Embodiment 1 of the present invention

[0054] Figure 4 Monitoring result graph of probability stationary subspace analysis-residual statistic in Embodiment 1 of the present invention ;

[0055] Figure 5 Monitoring result graph of probability stationary subspace analysis-stationary feature statistic in Embodiment 1 of the present invention ;

[0056] Figure 6 Monitoring result graph of dynamic non-stationary projection structure-residual statistic SPE in Embodiment 1 of the present invention

[0057] Figure 7 Monitoring result graph of dynamic non-stationary projection structure-stationary feature statistic in Embodiment 1 of the present invention ;

[0058] Figure 8 Monitoring result graph of dynamic non-stationary projection structure-non-stationary feature statistic in Embodiment 1 of the present invention ;

[0059] Figure 9 Monitoring result graph of stationary subspace analysis-Mahalanobis distance in Embodiment 2 of the present invention

[0060] Figure 10 Monitoring result graph of probability stationary subspace analysis-residual statistic in Embodiment 2 of the present invention ;

[0061] Figure 11 Monitoring result graph of probability stationary subspace analysis-stationary feature statistic in Embodiment 2 of the present invention ;

[0062] Figure 12 Monitoring result graph of dynamic non-stationary projection structure-residual statistic SPE in Embodiment 2 of the present invention

[0063] Figure 13 Monitoring result graph of dynamic non-stationary projection structure-stationary feature statistic in Embodiment 2 of the present invention ;

[0064] Figure 14 Monitoring result graph of dynamic non-stationary projection structure-non-stationary feature statistic in Embodiment 2 of the present invention ; Detailed implementation manners

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

[0066] As Figure 1As shown, in a preferred embodiment of the present invention, a fault detection method based on a dynamic non-stationary projection structure is provided. The method specifically includes the following steps:

[0067] Step 1: A training sample set is composed of the measured variables X collected during the chemical operation without faults; the process variable X needs to be a measured variable related to the existence of faults. For example, for the penicillin fermentation process in the subsequent embodiments, it may include variables such as air flow rate, stirring power, and bottom feed rate.

[0068] Step 2: The measured variables X in Step 1 are normalized to obtain the normalized training sample set X * .

[0069] Step 3: Based on the normalized training sample set X * , a dynamic non-stationary projection structure is constructed using a specific algorithm, and the projection structure is trained using a specific algorithm.

[0070] Step 4: Combining the normalized training sample set, multiple characteristic statistics of the projection structure are solved through the projection structure, and then the control limits of the multiple characteristic statistics of the projection structure are determined respectively.

[0071] Step 5: During online operation, the measured variables x(c) of the chemical operation to be measured are collected and normalized to obtain the normalized variable set to be measured.

[0072] Step 6: Combining the normalized variable set to be measured, the characteristic statistics in the chemical operation to be measured are determined using the dynamic non-stationary projection structure. Then, according to the magnitude relationship between the characteristic statistics in the chemical operation to be measured and the control limits of the multiple characteristic statistics of the dynamic non-stationary projection structure, it is judged whether there is a fault in the chemical operation to be measured.

[0073] The specific implementation methods and their effects of Steps 1 to 6 in this embodiment are described in detail below.

[0074] In Step 1 of this embodiment, the training sample set X * is expressed as:

[0075] X * = [x * (1), x * (2), x * (t)... x * (N)], t ∈ [1, N]

[0076] where x *(t) represents the measured variable of the t-th training sample after normalization, and N represents the number of training samples in the training sample set; each sample corresponds to a moment.

[0077] The multiple characteristic statistics of the projection structure include stationary characteristic statistics Non-stationary characteristic statistics Residual statistic SPE; the control limits of the multiple characteristic statistics of the projection structure respectively include the control limits of the stationary characteristic statistics Control limits of non-stationary characteristic statistics Control limits of residual statistics

[0078] In step 2 of this embodiment above, the normalization process is specifically determined by the following formula:

[0079]

[0080] Among them, X is the training sample set, R represents the set of real numbers, m is the number of measured variables; μ is the mean of all measured variables in the training sample set; δ is the standard deviation of all measured variables in the training sample set; X * is the normalized training sample set.

[0081] In step 3 of this embodiment above, the dynamic non-stationary projection structure constructed by a specific algorithm is as follows:

[0082]

[0083] Among them, B ∈ R m×p represents a linear superposition matrix, B s ∈ R m×a is composed of the first a columns of B, B n ∈ R m×(p-a) is composed of the last (p - a) columns of B, s s (t) ∈ R a×1 represents the stationary characteristic of the t-th training sample, s n (t) ∈ R (p-a)×1 represents the non-stationary characteristic of the t-th training sample, and represents the full characteristic matrix of the t-th training sample, where, R represents the set of real numbers, m represents the dimension of the measured variable, p is the dimension of the stationary characteristic s(t) of the t-th training sample, a is the dimension of the stationary characteristic s s (t) of the t-th training sample; e(t) represents the noise of the t-th training sample after normalization; e(t) follows a Gaussian distribution with a mean of 0 and a covariance of σ 2 I p of, I p represents the identity matrix of dimension p, σ2 Represents the noise intensity of the measurement variable.

[0084] Among them, the stationary feature s of the t-th training sample s The prior distribution of (t) is specifically as follows:

[0085]

[0086] In the formula, Represents a Gaussian distribution; p(s s (t)) represents the prior distribution of the stationary feature s of the t-th training sample s (t), μ s And Σ s Respectively represent the mean and covariance of the stationary feature s of the t-th training sample s (t).

[0087] Among them, the conditional prior distribution of the non-stationary feature s of the t-th training sample n (t) is as follows:

[0088]

[0089] In the formula, z(t) represents the hidden state of the hidden Markov chain, I represents the number of hidden states, z(t) = i represents that the non-stationary feature s of the t-th training sample n (t) is generated by the i-th hidden state, p(s n (t)|z(t) = i) represents the conditional prior distribution of the non-stationary feature s of the t-th training sample at time t when the hidden state is i n (t), And Respectively represent the mean and covariance of the non-stationary feature s of the t-th training sample at time t when the hidden state is i n (t); Note that the full feature of the t-th training sample Then the conditional prior distribution of the full feature s(t) of the t-th training sample is:

[0090]

[0091] In the formula, p(s(t)|z(t) = i) represents the conditional prior distribution of the full feature s(t);

[0092] Among them, Respectively represent the mean and covariance of the full feature s(t) of the t-th training sample at time t when the hidden state is i.

[0093] Among them, the initial distribution and state transition matrix of the hidden state z are specifically as follows:

[0094] p(z(1) = i) = π i, i = 1, 2, ..., I

[0095] p(z(t) = i|z(t - 1) = j) = a ij , i, j = 1, 2, ..., I, t = 2, 3, ..., N

[0096] In the formula, π i represents the probability that the hidden state at the initial moment is i, and I represents the number of hidden states. p(z(t) = i|z(t - 1) = j) represents the probability that the hidden state at time t - 1 is j and the hidden state at time t is i, and a ij represents the probability of transitioning from hidden state j to hidden state i. And there is A = (a ij ) ∈ R I×I , where A represents the state transition matrix, and a ij is the element in the j-th column of the i-th row of A.

[0097] Among them, the probability of generating the t-th sample when the hidden state at time t is i is:

[0098]

[0099] In the formula, b i (t) = p(x * (t)|z(t) = i) represents the probability of generating the t-th sample x * (t) when the hidden state at time t is i, and B represents the emission matrix.

[0100] In the formula, z(1) is the hidden state at the initial moment, p(z(1) = i) is the prior distribution of the initial hidden state z(1), z(t - 1) represents the hidden state at time t - 1, z(t) represents the hidden state at time t, and p(z(t) = i|z(t - 1) = j) represents the state transition probability distribution of the hidden state.

[0101] In step 3 of this embodiment above, a specific algorithm is used to train the dynamic non-stationary projection structure. The specific training algorithm includes the following process:

[0102] 3.1 Input the training sample set into the projection structure, and calculate the posterior hidden state and posterior joint hidden state of the hidden Markov chain through the forward-backward algorithm in combination with the current model parameters To calculate the above parameters, first define the forward variable and backward variable as follows:

[0103]

[0104]

[0105] In the formula, α i$\alpha(t)$ represents the forward variable when the hidden state at time $t$ is $i$. denotes the set consisting of the first to the $t$-th training samples. $\beta$ i $\beta(t)$ represents the backward variable when the hidden state at time $t$ is $i$. denotes the set consisting of the $(t + 1)$-th to the $N$-th samples. The calculations of the forward variable and the backward variable can be achieved through the forward algorithm and the backward algorithm, and the specific procedures are as follows in steps a) and b):

[0106] a) Initialization:

[0107]

[0108] $\beta$ i $(N)=1$

[0109] In the formula, $\alpha$ i $(1)$ represents the forward variable when the hidden state at the first time is $i$, $\pi$ i represents the prior probability that the hidden state at the first time is $i$, $b$ i $(1)$ represents the probability that the first sample is generated when the hidden state at the first time is $i$. $\beta$ i $(N)$ represents the backward variable when the hidden state at the $N$-th time is $i$.

[0110] b) Recursion:

[0111]

[0112]

[0113] In the formula, $\alpha$ i $(t + 1)$ represents the forward variable when the hidden state at time $t + 1$ is $i$, $b$ j $(t + 1)$ represents the probability that the $(t + 1)$-th sample is generated when the hidden state at time $t + 1$ is $j$; $\beta$ i $(t)$ and $\beta$ i $(t + 1)$ respectively represent the backward variables when the hidden state at the $t$-th time and the $(t + 1)$-th time is $i$.

[0114] After the above two steps, based on the forward variable, the backward variable, and the training samples, the posterior hidden state and the posterior joint hidden state can be calculated as follows:

[0115]

[0116]

[0117] In the formula, $\gamma$ i $(t)=p(z(t)=i|X$ * ) represents the posterior hidden state distribution at time $t$, $\alpha$ i$\alpha(t)$ represents the forward variable when the hidden state at time $t$ is $i$, and $b$ j $(t)$ represents the probability that the hidden state at time $t$ is $j$ for the $t$-th sample $x$ * $(t)$ occurs; $\xi$ ij $(t)==p(z(t)=i,z(t + 1)=j|X$ * ) represents the posterior joint hidden state distribution at time $t$, $\beta$ i $(t + 1)$ represent the backward variables when the hidden states at the $t$-th and $(t + 1)$-th times are $i$ respectively.

[0118] 3.2 Calculate the first moment of the local full features of the training samples and the second moment of the local full features. Among them, calculate the first moment and the second moment of the local full feature $s$ of the $t$-th training sample according to the current parameter set i $(t)$ as follows:

[0119]

[0120]

[0121] Among them, $s$ i $(t)$ represents the local full feature of the $t$-th training sample when the hidden state at time $t$ is $i$, and $\langle\cdot\rangle$ represents the operation of taking the expectation, represents the covariance matrix of the local full feature $s$ of the $t$-th training sample when the hidden state at time $t$ is $i$. i $(t)$.

[0122] 3.3 Combine the first moment of the local full features and the second moment of the local full features, and use the likelihood function to iterate and calculate repeatedly until the likelihood function converges, and update the current model parameters to obtain the optimal parameter set of the dynamic non-stationary projection structure Complete the training of the dynamic non-stationary projection structure. Among them, each parameter in the optimal parameter set is specifically determined by the following formula:

[0123]

[0124]

[0125]

[0126]

[0127]

[0128]

[0129]

[0130]

[0131] Among them, W s ∈R p×a is the first a columns of the identity matrix I p of dimension p, and W n ∈R p×(p-a) is the last (p - a) columns of the identity matrix I p , and Tr(*) represents the trace of a matrix.

[0132] In step 4 of the present embodiment, by combining the normalized training sample set through a dynamic non - stationary projection structure, based on the weighted Mahalanobis distance discrimination method, multiple characteristic statistics of the analysis model are solved to determine the control limits of the multiple characteristic statistics of the dynamic non - stationary projection structure, which specifically includes the following process:

[0133] First, determine the new local full feature of the t - th training sample and the posterior distribution of the new hidden state specifically as follows:

[0134]

[0135]

[0136]

[0137] In the formula, represents the covariance matrix of the new local full feature of the t - th training sample when the hidden state at time t is i , and respectively represent the first - order moment and covariance matrix of the new local full feature of the t - th training sample when the hidden state at time t is i , and respectively represent the prior mean and covariance matrix of the new stationary feature s s (t) of the t - th training sample, and respectively represent the prior mean and covariance matrix of the non - stationary feature s n (t) of the t - th training sample when the hidden state at the new time t is i represents the noise intensity of the new measurement variable, represents the new emission matrix.

[0138] During the parameter training phase, the samples at all times are already determined. Therefore, the posterior distribution of the hidden state can be determined by the forward-backward algorithm. However, for online testing, only the samples collected in the past and present are available, so only forward filtering can be performed. Therefore, during the statistic training phase, the statistics corresponding to each training sample should also be calculated in the online mode. The forward filtering is as follows:

[0139] 1) First, for the first training sample, the posterior distribution of the new hidden state is:

[0140]

[0141] In the formula, is the prior distribution of the hidden state corresponding to the new first sample, represents the probability of generating the first sample when the hidden state at the new initial time is i, and I represents the number of hidden states.

[0142] 2) For the sample at time t, the posterior distribution of the new hidden state is:

[0143]

[0144] In the formula, represents the posterior probability that the new hidden state at time t is i, represents the posterior probability that the new hidden state at time t - 1 is j, is the element in the i-th column of the j-th row of the new state transition matrix , representing the probability of the new hidden state j transitioning to the hidden state i, represents the probability of generating the t-th sample when the new hidden state at time t is i.

[0145] Then, combining the new local full feature of the t-th training sample and the new posterior probability of the discrete state variable, the optimal estimate of the feature corresponding to the t-th training sample is obtained by means of probability fusion, as follows:

[0146]

[0147] Combining the optimal estimate the residual e(t) of the t-th training sample is determined as follows:

[0148]

[0149] Since the residual e(t) follows a normal distribution with a mean of 0 and a covariance of The Gaussian distribution, so the residuals e = [e(1), e(2), …, e(N)] of the training sample set are a stationary time series. Thus, the change in process uncertainty can be monitored by combining the residual e(t) of the t-th training sample and the designed residual statistic SPE, specifically as follows:

[0150] SPE = [SPE(1), SPE(2), SPE(t)...SPE(N)], t ∈ [1, N]

[0151] Among them, SPE(t) represents the residual of the t-th training sample. The residual SPE(t) of the t-th training sample is calculated by the following formula:

[0152]

[0153] Among them, the new local full feature corresponding to the t-th training sample of the mean vector and the covariance matrix of the sum are calculated as follows respectively:

[0154]

[0155]

[0156] Among them, represents the new covariance matrix of the t-th training sample when the hidden state at time t is i; represents the new covariance of the observation noise of the training sample. Denote the covariance matrix of the new local full feature of the t-th training sample as P i .

[0157] The new local full feature of the t-th training sample when the hidden state at time t is i corresponding to the stationary feature The first moment and covariance matrix are determined according to the following formula:

[0158]

[0159]

[0160] Among them, E[·] represents solving the first moment, and respectively represent the first moment and covariance matrix of the stationary feature corresponding to the new local full feature of the t-th training sample when the hidden state at time t is i.

[0161] Based on the new local full feature of the t-th training sample The corresponding stationary feature Define the local Mahalanobis distance index as follows:

[0162]

[0163] where represents the new local full feature of the t-th training sample when the hidden state at time t is i The corresponding stationary feature The Mahalanobis distance between the first moment of itself

[0164] The stationary feature statistic defined by the weighted Mahalanobis distance Specifically as follows:

[0165]

[0166] where represents the stationary feature statistic of the t-th training sample, specifically as follows:

[0167]

[0168] where represents the posterior probability that the hidden state of the t-th training sample is i

[0169] At the same time, the new local full feature of the t-th training sample when the hidden state at time t is i The corresponding non-stationary feature The first moment and covariance matrix are determined according to the following formula:

[0170]

[0171] where and respectively represent the first moment and covariance of the non-stationary feature corresponding to the new local full feature of the t-th training sample when the hidden state at time t is i represents the prior mean of the non-stationary feature s n (t) of the t-th training sample when the hidden state at time t is i

[0172] Based on the above new local full feature of the t-th training sample when the hidden state at time t is i The corresponding non-stationary feature Define the local Mahalanobis distance index specifically as follows:

[0173]

[0174] where represents the new local full feature of the t-th training sample when the hidden state at time t is i The corresponding non-stationary feature The Mahalanobis distance between the local mean of itself

[0175] Define the non-stationary feature statistic through the weighted Mahalanobis distance Specifically as follows:

[0176]

[0177] Among them, Represents the non-stationary feature statistic of the t-th training sample, specifically as follows:

[0178]

[0179] Finally, combined with the confidence level, use the kernel density estimation method to determine the control limits of the multiple feature statistics of the analysis model, including the control limits of the stationary feature statistics The control limit of the non-stationary feature statistic The control limit of the residual statistic

[0180] In step 5 of the above embodiment of the present invention, the normalized set of variables to be measured Specifically determined according to the following formula:

[0181]

[0182] In the formula, Represents the c-th sample to be measured after normalization

[0183] In step 6 of the above embodiment of the present invention, referring to the stationary feature statistic in step four above The non-stationary feature statistic The calculation method of the residual statistic SPE, calculate the various feature statistics in the chemical operation process to be measured, including the stationary feature statistic to be measured The non-stationary feature statistic to be measured The process residual statistic SPE(c) to be measured. The specific calculation process is as follows:

[0184] For the c-th sample The local full feature, local stationary feature and local non-stationary feature are estimated as:

[0185]

[0186]

[0187]

[0188] In the formula, Represents the c-th sample when the hidden state at time c is i The local full feature represents the c-th sample when the hidden state at time c is i The local stationary feature represents the c-th sample when the hidden state at time c is i The local non-stationary feature

[0189] The posterior distribution of the hidden variable at the initial time is

[0190]

[0191] In the formula, represents the posterior probability that the hidden state at the initial time is i represents the probability of generating the first sample to be measured when the hidden state at the initial time is i

[0192] For the c-th sample The posterior distribution of the hidden variable at time c is

[0193]

[0194] In the formula, represents the posterior probability that the hidden state at time c is i represents the posterior probability that the hidden state at time c - 1 is j represents the probability of generating the c-th sample to be measured when the hidden state at the initial time is i

[0195] Then, combining the new local full feature of the c-th sample to be measured and the new posterior probability of the discrete state variable The optimal estimate of the full feature corresponding to the c-th sample to be measured is obtained by using the probability fusion method Specifically as follows:

[0196]

[0197] Combining the optimal estimate Determine the residual e(c) of the c-th sample to be measured, specifically as follows:

[0198]

[0199] Then, the residual statistic corresponding to the c-th sample to be measured is calculated as follows:

[0200]

[0201] Based on the local stationary feature and non-stationary feature of the c-th sample to be measured, calculate the local Mahalanobis distance index as follows:

[0202]

[0203]

[0204] In the formula, represents the local stationary feature of the c-th sample to be measured when the hidden state at time c is i the Mahalanobis distance between itself and its first moment, represents the local non-stationary feature of the c-th sample to be measured when the hidden state at time c is i the Mahalanobis distance between itself and its first moment

[0205] Then, for the c-th sample to be measured, the corresponding and statistic is calculated as follows:

[0206]

[0207]

[0208] In step 6 of this embodiment, determining whether there is a fault in the chemical process to be measured specifically is:

[0209] If each of the characteristic statistics in the chemical process to be measured is less than the control limit of each of the multiple characteristic statistics determined in step 4, that is and and that is, the stationary characteristic statistic in the chemical process to be measured is less than the control limit of the stationary characteristic statistic and the non-stationary characteristic statistic in the chemical process to be measured is less than the control limit of the non-stationary characteristic statistic independent of quality and the residual statistic SPE(c) in the chemical process to be measured is less than the control limit of the residual statistic then there is no fault in the chemical process to be measured, otherwise there is a fault in the chemical process to be measured, thus completing the fault detection of the chemical process to be measured.

[0210] The following is illustrated based on two embodiments. The data of these two embodiments is from a closed-loop controlled penicillin fermentation process, as Figure 2 shown. 1200 samples of the normal fermentation process are collected, and the sampling interval is 0.1 h. Each sample consists of 10 process variables, as shown in Table 1 specifically.

[0211] Table 1 Monitoring variables of penicillin fermentation process

[0212]

[0213] In two embodiments of the present invention, two typical faults in the penicillin fermentation process are considered respectively, including ramp fault and step fault, as shown in Table 2. For these two faults, 1200 samples are collected for verifying the performance of the algorithm. The faults are introduced starting from the 601st sample and continue until the last sample.

[0214] Table 2 Two typical faults in the penicillin fermentation process

[0215]

[0216] In Embodiment 1, the tested fault is the increased air flow rate fault, the fault type is step fault, and the fault amplitude is 5%. Figures 3 - 8 The monitoring results of different algorithms are summarized, including the stationary subspace analysis method, the probabilistic stationary subspace analysis method, and the dynamic non-stationary projection structure method described in the present invention. As Figure 3 shown, since the stationary subspace analysis method does not consider the uncertainty of the process, the fault undetected rate is the highest, about 10%. At the same time, as Figure 4 and Figure 5 shown, the detection rate of the probabilistic stationary subspace analysis method reaches 97.2%. Since the dynamic non-stationary projection structure models the uncertainty and dynamics of the process and strengthens the dynamic information in the non-stationary characteristics, the model can monitor from both the static and dynamic characteristics of the process, realizing comprehensive fault detection. As Figures 6 - 8 shown, in terms of the fault detection rate, both the residual statistic and the stationary feature statistic of the dynamic non-stationary projection structure method are higher than those of the probabilistic stationary subspace analysis method. In addition, the non-stationary feature statistic of the dynamic non-stationary projection structure method has the highest fault detection rate.

[0217] In Embodiment 2, the tested fault is the rising stirring power fault. Specifically, the stirring power rises by 0.5% per hour. The monitoring results of the three methods are as Figures 9 - 14 shown. In this example, the defect of the stationary subspace analysis is further demonstrated. It can be seen that there is an obvious undetected phenomenon in the Mahalanobis distance, and about two-thirds of the faults are ignored. As Figure 10As shown, since the probability stationary subspace analysis models the process uncertainty, the fault detection rate has been improved to a certain extent, and the residual statistic of this method detects 86.7% of the faults. In contrast, the non-stationary feature statistic of the dynamic non-stationary projection structure can detect 96.7% of the faults and can give a timely prediction at the initial stage of fault development. This is because although the amplitude in the early stage of the fault is not high, under the adjustment of the feedback control, the dynamic characteristics of the process have changed. Therefore, the dynamic non-stationary projection structure successfully separates the stationary trend, non-stationary trend and measurement noise, and emphasizes the dynamic characteristics when extracting non-stationary features, providing a complete monitoring framework for industrial process monitoring, thus enhancing the detection ability of the algorithm for non-stationary processes.

[0218] The above embodiments are used to explain the present invention rather than limit the present invention. Any modifications and changes made to the present invention within the spirit and scope of the claims of the present invention fall within the protection scope of the present invention.

Claims

1. A fault detection method based on a dynamic non-stationary projection structure, characterized in that: The method includes the following steps: Step 1: Compose a training sample set from the measurement variables collected during the fault - free chemical operation process that are relevant to the existence of faults; Step 2: Normalize the measurement variables of all training samples in the training sample set in Step 1 to obtain a normalized training sample set; Step 3: Based on the normalized training sample set, construct a dynamic non - stationary projection structure and train the projection structure; Step 4: Combine the normalized training sample set to determine the control limits of each of the multiple characteristic statistics of the projection structure by solving the multiple characteristic statistics of the projection structure; Step 5: Online collect the measurement variables to be measured during the chemical operation process to be measured and perform normalization processing to obtain a normalized set of variables to be measured; Step 6: Combine the normalized set of variables to be measured and use the dynamic non - stationary projection structure to determine the characteristic statistics during the chemical operation process to be measured, and then judge whether there are faults during the chemical operation process to be measured according to the magnitude relationship between the characteristic statistics during the chemical operation process to be measured and the control limits of each of the multiple characteristic statistics of the dynamic non - stationary projection structure; The specific method for constructing the dynamic non - stationary projection structure is as follows: where \(B\in R\) m×p represents a linear superposition matrix, \(B\) s \(\in R\) m×a is composed of the first \(a\) columns of \(B\), \(B\) n \(\in R\) m×(p-a) is composed of the last \((p - a)\) columns of \(B\), \(s\) s (t)\(\in R\) a×1 represents the stationary features of the \(t\)-th training sample, \(s\) n (t)\(\in R\) (p-a)×1 represents the non-stationary features of the \(t\)-th training sample, and represents the full feature matrix of the \(t\)-th training sample, where \(R\) represents the set of real numbers, \(m\) represents the dimension of the measurement variable, \(p\) is the dimension of the full feature \(s(t)\) of the \(t\)-th training sample, and \(a\) is the dimension of the stationary feature \(s\) s (t) of the \(t\)-th training sample; \(e(t)\) represents the noise of the \(t\)-th training sample after normalization; The specific method for training the projection structure is as follows: First, input the training sample set into the projection structure, and calculate the posterior hidden state and posterior joint hidden state of the hidden Markov chain through the forward-backward algorithm combined with the current model parameters where A represents the state transition matrix, B represents the emission matrix, π i represents the probability that the hidden state at the initial moment is i, μ s represents the mean of the stationary features of the training samples, Σ s represents the covariance of the stationary features of the training samples, represents the mean of the non-stationary features of the training samples when the hidden state is i, represents the covariance of the non-stationary features of the training samples when the hidden state is i, σ 2 represents the noise intensity of the measurement variable; then calculate the first moment of the local full features and the second moment of the local full features of the training samples; then, combined with the posterior hidden state, posterior hidden state transition variable, first moment of the local full features, and second moment of the local full features, use the likelihood function to iterate and calculate repeatedly until the likelihood function converges, and update the current model parameters to obtain the optimal parameter set of the dynamic non-stationary projection structure Complete the training of the dynamic non-stationary projection structure, where are A, B, π i , μ s , Σ s , σ 2 the updated optimal parameters.

2. The fault detection method based on a dynamic non-stationary projection structure according to claim 1, wherein: The training sample set X * is expressed as: X * = [x * (1), x * (2), x * (t)... x * (N)], t ∈ [1, N], where x * (t) represents the measurement variable of the t-th training sample after normalization processing, N represents the number of training samples in the training sample set; each sample corresponds to a moment.

3. The fault detection method based on a dynamic non-stationary projection structure according to claim 1, characterized in that: The 8 parameters in the optimal parameter set correspond one - to - one with the 8 parameters in the model parameter Θ, and are specifically determined by the following formula: where, γ i (t) represents the posterior hidden state distribution at time t, ξ ij (t) represents the posterior joint hidden state distribution at time t, x * (t) represents the sample at time t, i.e., the t-th sample, <s i (t)> and represent the first moment and the second moment of the local full feature s i (t) of the t-th training sample respectively, <·> represents the operation of taking the expectation, s i (t) represents the local full feature of the t-th training sample when the hidden state at time t is i; I represents the number of hidden states, W s ∈R p×a is the first a columns of the identity matrix I p of dimension p, W n ∈R p×(p-a) is the last (p - a) columns of the identity matrix I p and Tr(*) represents the trace of the matrix.

4. The fault detection method based on the dynamic non-stationary projection structure according to claim 3, characterized in that: The multiple characteristic statistics of the projection structure include the stationary characteristic statistic T s 2 , the non-stationary characteristic statistic T n 2 , and the residual statistic SPE; the control limits of the multiple characteristic statistics of the projection structure respectively include the control limit of the stationary characteristic statistic the control limit of the non-stationary characteristic statistic the control limit of the residual statistic The calculation formula for the residual statistic SPE is: SPE = [SPE(1), SPE(2), SPE(t)... SPE(N)], t ∈ [1, N] Where, SPE(t) represents the residual of the t - th training sample and is calculated by the following formula: In the formula: e(t) represents the residual of the t - th training sample; The steady feature statistic T s 2 has the following calculation formula: T s 2 = [T s 2 (1), T s 2 (2), T s 2 (t)...T s 2 (N)], t ∈ [1, N] where, T s 2 (t) represents the stationary feature statistic of the t-th training sample, which is calculated by the following formula: In the formula: represents the posterior probability that the hidden state corresponding to the t-th training sample is i, represents the new local full feature of the t-th training sample when the hidden state at time t is i corresponding stationary feature the Mahalanobis distance between itself and the first moment; The non-stationary feature statistic T n 2 has the following calculation formula: T n 2 = [T n 2 (1), T n 2 (2), T n 2 (t)...T n 2 (N)], t ∈ [1, N] Among them, represents the non-stationary feature statistic of the t-th training sample, and the specific calculation formula is as follows: In the formula: represents the new local full feature of the t-th training sample when the hidden state at time t is i corresponding non-stationary feature Mahalanobis distance between itself and its local mean.

5. The fault detection method based on the dynamic non-stationary projection structure according to claim 4, characterized in that: The control limits of each of the multiple characteristic statistics of the projection structure are determined by combining the confidence level and using the kernel density estimation method.

6. The fault detection method based on a dynamic non-stationary projection structure according to claim 1, wherein: Each characteristic statistic in the chemical process to be measured includes the stationary characteristic statistic to be measured The non-stationary characteristic statistic to be measured The residual statistic SPE(c) to be measured, and its calculation method is the same as that of the stationary characteristic statistic T in Step 4 s 2 and the non-stationary characteristic statistic T n 2 and the residual statistic SPE are the same.

7. The fault detection method based on a dynamic non-stationary projection structure according to claim 6, characterized in that: The specific judgment in Step 6 on whether there are faults during the chemical operation process to be measured is as follows: If all the characteristic statistics in the chemical process to be measured are respectively less than the corresponding control limits of the characteristic statistics determined in step 4, that is, the stationary characteristic statistics in the chemical process to be measured are less than the control limit of the stationary characteristic statistics and the non-stationary characteristic statistics in the chemical process to be measured are less than the control limit of the non-stationary characteristic statistics and the residual statistic SPE(c) in the chemical process to be measured is less than the control limit of the residual statistic then there is no fault in the chemical process to be measured; otherwise, there is a fault in the chemical process to be measured, thus completing the fault detection of the chemical process to be measured.

Images