A Multimodal Monitoring Method Based on Semi-Supervised Discriminant Hybrid Probabilistic Principal Component Analysis

Abstract

The present invention relates to a multi-modal monitoring method based on semi-supervised discriminant mixture probabilistic principal component analysis, which includes two stages: offline modeling and online modeling. The discriminant information in the labeled and unlabeled data samples is used to divide the multi-modal space, and the global and local geometric information of the data set is utilized to enhance the separability of the modal space, which is beneficial to improving the accuracy of the model. Drawing on the idea of mixture probabilistic principal component analysis (MPPCA), the E-F-M algorithm is proposed to automatically determine the dimension of the low-dimensional space of the mixture model and obtain the optimal model parameters. For each modality, a weighted SVDD model and Bayesian inference are used to construct the monitoring statistic, which fully utilizes the data modality information and solves the problem that the detection accuracy rates of the traditional fault monitoring statistics T<supgt;2< / supgt; and SPR are not high due to the presence of process factors and non-process factors in the process data. This method can not only monitor known process modalities but also detect unknown process modalities. The present invention automatically divides a complex data set into multiple modalities without prior modality division, which is beneficial to expanding the applicable scope of the algorithm.

Description

Technical Field

[0001] The present invention belongs to the technical field of industrial process monitoring, and particularly relates to a multi-modal monitoring method based on semi-supervised discriminant hybrid probabilistic principal component analysis. Background Technique

[0002] With the rapid development of modern industry, the scale of the production process shows a rapid expansion trend. Due to the complexity of the working environment and product specifications, product quality and production safety have become crucial. The inherent non-linearity, multi-modality, time-variation, uncertainty and other characteristics of complex processes make it difficult for mechanism-based models to adapt to the complex characteristics of actual industrial processes. With the wide use of networked instruments in industrial processes and the development of big data technology and artificial intelligence, data-driven process fault detection and diagnosis methods can discover process health knowledge from a large amount of redundant production process data and have become the mainstream technology in the current monitoring field, and have been successfully applied to production processes such as pharmaceuticals, plastics, metallurgy, semiconductor manufacturing, etc. [References 1-5].

[0003] In the past decade or so, multivariate statistical methods have been widely used in process monitoring and fault diagnosis. Common methods mainly include principal component analysis PCA, partial least squares PLS, Fisher discriminant analysis FDA, slow feature analysis, kernel methods and related improvements [References 1, 2, 6-12]. These methods model the normal state from historical data and use Hoteloing T 2and mean square prediction error (SPE) statistics to identify process faults. However, in practice, due to various factors such as changes in production strategies, migration of operating states, and fluctuations in raw materials, process data obey non-Gaussian distributions, and process data have strong multimodal characteristics. In addition, process data often have high-dimensional sparsity and time-varying characteristics, and data reflecting the operating state of the process are often located in a low-dimensional feature space. At present, the solutions are mainly divided into nonlinear modeling methods, independent principal component analysis (ICA), hybrid models, multi-model methods based on ensemble learning, and process monitoring methods based on real-time learning [references 1-3, 12-14]. The nonlinear modeling process monitoring method establishes monitoring statistics based on the global nonlinear process model, and its monitoring accuracy depends on the accuracy of nonlinear characterization. However, it is difficult for global nonlinear processes to achieve satisfactory monitoring performance under different conditions, and the interpretability is poor. The multi-model method adopts a "divide and conquer" strategy to model the local areas of the process using multiple linear models separately, and then uses the model output to fuse as the nonlinear process monitoring result. This type of method not only overcomes the shortcomings of the global nonlinear modeling method, but also can effectively improve the process monitoring performance by considering the importance of local models. ICA can extract high-order and low-order statistics from non-Gaussian distribution data to achieve accurate modeling of the process. However, this type of method is difficult to characterize the changing operating state, the collinearity problem of process variables and the process uncertainty, and it is difficult to achieve satisfactory monitoring performance. The instant method establishes a local model for a certain number of adjacent or similar samples, and then proposes a redundant model. This type of method can track the changes in the process in a timely manner, but there are problems such as large computational complexity and collinearity.

[0004] Since labeling data sample categories requires rich process practice experience and expensive human and material costs, the amount of labeled data is far less than the amount of unlabeled data. How to make full use of the semi-supervised learning process monitoring method of labeled data and unlabeled data to improve monitoring performance has always been an important issue for research in industry and academia. PCA has clear physical meaning, simple fault detection statistics, and low computational complexity. Many scholars have done a lot of theoretical research and application in fault monitoring, positioning, robustness, sensor location, multimodal data modeling, etc. [References 1,8,13,15,16]. However, the existing PCA and its improved algorithms are generally unsupervised algorithms, and the extracted features do not consider the identification information between modes. In addition, considering the differences in the distribution of process data in different local areas of the process and the uncertainty of process data, multi-model hybrid methods have become the main method for nonlinear and multi-modal process monitoring. Based on the above discussion, considering the locality of multi-model data and the effectiveness of feature extraction by discriminant analysis [References 13,16-20], the present invention provides a new multi-modal process monitoring method of semi-supervised hybrid discriminant probability principal component analysis. Summary of the invention

[0005] The object of the present invention is to design a multi-modal process monitoring method based on semi-supervised hybrid discriminant probabilistic principal component analysis for the deficiencies existing in the above-mentioned prior art.

[0006] To achieve the above object, the technical solution of the present invention is as follows: A multi-modal monitoring method based on semi-supervised discriminant hybrid probabilistic principal component analysis includes two stages: offline modeling and online modeling. The specific process of the offline modeling stage is as follows:

[0007] S101. Input the labeled data set X L , y and the unlabeled data set X u , and call the semi-supervised discriminant MPPCA method to train the labeled data set and the unlabeled data set to obtain the training results u k , ψ k , k = 1, 2, …, K, m = 1, 2, …, N u ;

[0008] S102. Set the confidence value γ and the regularization constant C;

[0009] S103. Divide the unlabeled samples into the corresponding modalities. The determination method for the sample belonging to the k * -th modality is

[0010] The corresponding weight is and represents the indicator variable belonging to the k-th modality. If then it means completely belongs to the k-th modality, otherwise it does not belong to the k-th modality at all; represents the random variable belonging to the k-th modality, represents the expectation of the random variable ; Let the samples of the labeled samples and the unlabeled sample set belonging to the k-th class be The probability of each sample belonging to this class is denoted as where the weight of the labeled samples is 1, i = 1, 2, …, N k , N k is the number of samples in this class;

[0011] S104. Solve the optimization problem shown based on the latent variable set corresponding to the data set and calculate the center of the hypersphere in the latent space and the radius k = 1, 2, …, K.

[0012] S105. Calculate according to the following formula the reconstruction residual e′ corresponding to the sample x′ in the set k，i k，i ,

[0013]

[0014] In the formula, I is the identity matrix of the corresponding dimension, uk are respectively the latent variable loading matrix and the mean of the k-th sub-model, ψ k is the noise covariance matrix of the k-th sub-model, and k = 1, 2,..., K;

[0015] According to the residual data set Solve the optimization problem, and then calculate the center of the hypersphere in the residual space according to the optimal parameters and the radius k = 1, 2,..., K;

[0016] The specific process of the online modeling stage is as follows:

[0017] S201. For the newly input data x new , calculate its expectation projected onto the k-th (k = 1, 2,..., K) latent variable and the sample reconstruction error e new , k, k = 1, 2,..., K;

[0018] S202. Solve the posterior probability belonging to the abnormal state, and solve the posterior probability that the sample reconstruction error e new，k belongs to the abnormal state, k = 1, 2,..., K;

[0019] S203. Calculate the global Bayeisan fault monitoring index BFDI st (x train ) and BFDI e (x train );

[0020] S204. Calculate the confidence limits BFDI e，lim and BFDI st，lim of the global fault monitoring index from all training samples;

[0021] S205. Calculate the global Bayeisan fault monitoring index BFDI new of x st (x new ) and BFDI e (x new ); ​

[0022] S206. Determine whether a fault has occurred. If BFDI st (x new ) > BFDI st，lim or BFDI e (x new ) > BFDI e，lim , then the process has a fault; otherwise, the process is in a normal state.

[0023] The further refined technical solution of the present invention is as follows:

[0024] Preferably, in step 101, the specific process of training the semi-supervised discriminant MPPCA method on the labeled dataset and the unlabeled dataset is as follows:

[0025] (1) Calculate the mean and variance of the samples, and normalize the data; set the dimension d' of the low-dimensional latent space, and the initial latent space dimension is d = D - 1 (D is the dimension of the original input space); set the initial value a l = 1 and l = 1, 2,..., d', the maximum number of iterations maxiter, the initial value t' of the number of iterations is 1, the algorithm stop threshold ε, and the number of modes K;

[0026] (2) Use the MPPCA algorithm to cluster all samples, and take the mean and variance of each low-dimensional subspace sample as the initial values of the prior distribution parameters of the latent variable t of this class. The MPPCA model parameters are used as the load matrix P k , the initial value of the class prior probability . The mean of the samples in each subspace and the variance of the reconstruction error are used as the mean u k of the k-th class latent vector and the initial value of the noise variance Ψ k (set as a diagonal matrix);

[0027] (3) E-step: Calculate the corresponding expected value of the latent variable and the second-order statistic using the current model parameters;

[0028] (4) F-step: Update the projection discriminant information matrix P';

[0029] (5) M-step: Update the model parameters P k , α l , u k , ψ k , k = 1, 2,..., K;

[0030] (6) If the value of the log-likelihood function satisfies |L(X, Y, X u |Θ t'-1 ) - L(X, Y, X u |Θ t')| ≤ ε or the number of iterations satisfies t' ≥ maxiter, then go to step (7); otherwise, set t' = t' + 1 and return to step (3);

[0031] (7) For a l = 1, sort them in descending order, and retain the corresponding first d' latent variables and model parameter values u k , ψ k as the optimal model parameters.

[0032] Preferably, in the step S104, the dataset corresponding latent variable t k,i , then the weighted SVDD model parameters can be obtained from the following optimization problem,

[0033]

[0034]

[0035] where L(α k ) is a function of the function optimization variable α k , N' k represents the number of data samples in the k-th class mode, α k,i、 α k,j are all elements in α k , t' k,j , t' k,i are all latent variables corresponding to the dataset , represents the weighting of the slack variable in the original SVDD optimization problem, and i, j = 1, 2,..., N' k , C is the regularization constant; let represent the support vector set, let be the number of support vectors of the k-th SVDD model, and the support vector set located on the spherical boundary is denoted as is the boundary support vector of the k-th SVDD model, and the center of the hypersphere is calculated in the following form and the radius

[0036]

[0037]

[0038] where, is any boundary support vector of the k-th class and h k，i″ , h k，i″ are respectively elements of the set , is the set The number of elements in

[0039] Preferably, in the step S105, according to the residual data set Solve the optimization problem similar to the formula (34),

[0040]

[0041]

[0042] In the formula, is the vector to be optimized in the formula (34'), is the coordinate in, e' k，i , e' k，j represent the data in the residual data set, i, j = 1, 2,..., N' k ;

[0043] Then, according to the optimal parameter and the formulas (35') to (36'), calculate the center and radius of the hypersphere in the residual space, k = 1, 2,..., K.

[0044]

[0045]

[0046] Wherein, is an arbitrary boundary support vector of the k-th class and are the support vector set and the boundary support vector set respectively, has the same meaning as the formula (34'), are both elements of, is the set the number of elements in, The superscript T represents the transpose of a matrix or vector.

[0047] Preferably, in the step S201, for the newly input data x new , calculate its projection to the expectation of the latent variable in the k-th (k = 1, 2,..., K) subspace and the sample reconstruction error e new，k , k = 1, 2,..., K,

[0048]

[0049]

[0050] In the formula, I is the identity matrix of the corresponding dimension, u k are the latent variable load matrix and mean of the k-th sub-model, respectively. ψ k is the noise covariance matrix of the k-th sub-model, and k = 1, 2, …, K.

[0051] Preferably, in the step S202, solve according to the formulas (40) to (41) the posterior probability of the abnormal state,

[0052]

[0053]

[0054] wherein, respectively represent the prior knowledge about the abnormal operating state and the normal state in the k-th latent space, is the posterior probability of the abnormal state, is the conditional probability distribution of under the abnormal state, is the prior distribution of is the conditional probability distribution of under the normal state;

[0055] Similarly, for the residual e new,k , its posterior probability p(Fau|e new,k ) of belonging to the abnormal state is calculated by the following formula,

[0056]

[0057]

[0058] wherein, p(e new,k |Fau) is the conditional distribution probability of the residual e new,k under the abnormal state, p(e new,k ) is the prior distribution of the residual e new,k , p(e new,k |Nor) is the conditional probability distribution of the residual e new,k under the normal state, are the prior probabilities of the k-th mode on the residual space regarding the fault and the normal state, respectively, k = 1, 2, …, K.

[0059] Preferably, in the step S203, calculate the global Bayeisan fault monitoring index BFDIs of all training samples according to the formula in step 6 above, that is, BFDI st (x train ) and BFDI e(x train ),

[0060]

[0061]

[0062] wherein, k = 1, 2, …, K, p(k) = N' k / N, is the posterior probability distribution belonging to mode k, p(k|e new，,k ) is the posterior probability distribution of e train,k belonging to mode k, k' represents the k'-th mode, k' = 1, 2, …, K, represents the marginal probability, p(k') represents the prior probability of mode k', p(e train，,k' ) represents the marginal probability of e train，,k' .

[0063] Preferably, in the step S204, the global fault monitoring index confidence limits BFDI e,lim and BFDI st,lim are calculated from all the training samples of formulas (46) to (47),

[0064]

[0065]

[0066] Preferably, in the step S205, the global Bayesian fault monitoring index BFDI new of x st (x new ) and BFDI e (x new ) are calculated according to formulas (44) to (45),

[0067]

[0068]

[0069] wherein,

[0070]

[0071] p(k) = N' k / N, is the posterior probability distribution belonging to mode k, p(k|e new，,k ) is the posterior probability distribution of e new,k belonging to mode k, k' represents the k'-th mode, denote the marginal probability, p(k') represents the prior probability of mode k', p(e new，,k' ) represents the marginal probability of e new，,k' , where k' = 1, 2, …, K.

[0072] Compared with the prior art, the present invention has the following advantages:

[0073] (1) The method of the present invention uses the discriminant information in the labeled and unlabeled data sample sets to divide the multi-modal space, and utilizes the global and local geometric information of the data set to enhance the separability of the modal space, which is beneficial to improving the accuracy of the model;

[0074] (2) The present invention draws on the idea of the mixture probability principal component analysis MPPCA, and proposes to use the E-F-M algorithm to automatically determine the dimension of the low-dimensional space of the mixture model and obtain the optimal model parameters;

[0075] (3) For each mode, a weighted SVDD model and Bayesian inference are used to construct the monitoring statistic, which fully utilizes the data modal information to solve the problem that the traditional fault monitoring statistics T 2 and SPE have low detection accuracy rates due to the inclusion of process factors and non-process factors in the process data. This method can not only monitor known process modes but also detect unknown process modes.

[0076] In summary, the present invention automatically divides a complex data set into multiple modes without prior mode division, which is beneficial to expanding the application scope of the algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] The present invention will be further described below with reference to the accompanying drawings.

[0078] Figure 1 is the probability structure diagram of semi-supervised mixture probability principal component analysis in the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0079] The present invention will be described in detail below in conjunction with the embodiments and the accompanying drawings of the specification.

[0080] Embodiment 1 Semi-supervised Mixture Discriminant Probability Principal Component Analysis

[0081] Assume that there are K classes of data, the labeled observed data is and the unlabeled data y i ∈ {1, …, K} is the class to which the data belongs, N l and N u are the data volumes of the labeled samples and unlabeled samples respectively. Let Assume that the latent variables and observed variables have the following linear model:

[0082] x = Pt + u + ε (1)

[0083] Where, is the weight matrix similar to the PCA loading matrix and satisfies u is the bias, is the latent vector and d < D. Its physical meaning is similar to the score vector of PCA. ε is the observation noise. To better represent the class of data, the class indicator vector is defined here Where z k ∈ {0, 1} is binary data. z k = 1 indicates that the observed data x comes from class k th satisfies For the labeled data (x, y), denote the class indicator vector as When y = k Otherwise

[0084] For the unlabeled sample data x u , introduce its class indicator vector satisfies Since z u is an unknown latent variable, the prior probability that the sample x u comes from the k-th model is:

[0085]

[0086] Here k = 1, 2, …, K. Assume that the observation noises are independent of each other and follow the probability distribution:

[0087] ε|z k = 1 ~ N(0, ψ k ) (3)

[0088] Where, ψ k is a diagonal matrix and does not need to follow the assumption of consistency of the noise level of each variable. z k can represent or Assume that the prior distribution of the latent variable t in each class latent space is a Gaussian distribution, that is:

[0089] p(t|z k = 1) = N(μ k, Σ k ) (4)

[0090] Where, N(μ k , Σ k ) represents the Gaussian distribution function. μ k and Σ kThey are the mean and covariance matrix of the k-th type of latent vector respectively. Then the conditional probability of the observed data x with respect to t and z k is as follows:

[0091] p(x|t,z k =1)=N(P k t k ,ψ k ) (5)

[0092] The marginal probability of the observed data x with respect to the k-th class is:

[0093] p(x|z k =1)=∫p(x|t,z k =1)p(t|z k =1)t = dt = N(m' k ,S' k ) (6)

[0094] where m' k =P k μ k , Generally speaking, the observed data is generally located in a high-dimensional space, and its corresponding hidden variables are often located in a low-dimensional feature space, satisfying d << D.

[0095] Suppose the prior probabilities of the labeled samples and unlabeled samples from the k-th sub-model are respectively and satisfying and Then the probability of the labeled sample x indicating vector is:

[0096]

[0097] For the unlabeled sample x u , the following form can also be obtained:

[0098]

[0099] where z u represents the class indicator variable of x u , and satisfying

[0100] Then for the labeled sample (x, y), the marginal distribution function is

[0101]

[0102] Note that here is known, and its value is determined by the label y value. And for the unlabeled data x u, whose marginal distribution is modeled using a Gaussian mixture model, i.e.,

[0103]

[0104] introduce the class indicator latent variable z u , then the marginal distribution can be transformed into the following form:

[0105]

[0106] p(x u |Θ) = ∫p(x u |Θ, t u , z u )p(t u |Θ, z u )p(z u )dt u dz u (12)

[0107] Note that the complete likelihood function form of the unlabeled samples is where represents the latent variable corresponding to x u in the k-th model. Let T u and T represent the latent variables of the labeled samples and unlabeled samples output by the mixture model respectively, and T k represent the latent variables of the labeled samples and unlabeled samples by the k-th model respectively, and Z represents the class label indicator vector of the labeled samples and Z = [z1, z2, …, z N . In Figure 1 , the quality data acquisition period is N'T s , where T s is the data sampling period. As Figure 1 shown, in the actual industrial process, generally quality data is scarcer than unlabeled data, and obtaining its quality data often requires spending manpower and material resources. Not only is its period much larger than that of unlabeled data, but the number of quality data samples is also very scarce. According to Figure 1 the probability structure diagram of the supervised mixture probabilistic principal component analysis algorithm shown, the complete likelihood function form of the model is:

[0108]

[0109] Thus, the optimal model parameters Θ * can be obtained by maximizing equation (13). The above optimization problem is a non-linear optimization problem, and its solution method is similar to the EM algorithm used by MPPCA. Then the above complete likelihood function with respect to the model latent variables T u, T, Zu The expected function form is

[0110]

[0111] It can be iteratively optimized through the E-step and M-step.

[0112] E-step: Calculate the posterior distribution of the latent variable

[0113] The E-step of the EM algorithm is to fix the model parameter Θ and calculate the posterior probabilities of the latent variables corresponding to the labeled data and unlabeled data and p(t m,k |x m , Θ k , z m,k = 1) and They can calculate the posterior probability of the class indicator latent variable from the Bayesian posterior probability formula, that is

[0114]

[0115] Note that for simplicity here is abbreviated as is the prior of the k-th class. The conditional distribution probability can be calculated by Equation (6).

[0116] The model latent variable t m,k and The calculation form of the posterior distribution probability is as follows:

[0117]

[0118]

[0119] Note that the class indicator variable z in the formula m,k and means that z m,k = 1, For simplicity, abbreviations are made in the text. The unlabeled data is denoted as The labeled data is denoted as If z m,k = 1 holds, then the function Otherwise it is equal to 0. Since the prior of the model latent variable and the conditional distributions of the labeled data and unlabeled data are all Gaussian distributions, then the posterior probability of the latent variable corresponding to the labeled data must also be a Gaussian distribution, and its first-order statistic and second-order statistic forms are as follows:

[0120]

[0121]

[0122] Among them, Let According to equations (18) - (19), it can be known that the label sample x m The mean and variance of the corresponding hidden vector t m,k in the k-th model are μ m,k and The posterior probability of the latent variable corresponding to the unlabeled data must also be a Gaussian distribution, and its first-order and second-order statistics are in the form of:

[0123]

[0124]

[0125] Among them, Let According to equations (18) - (19), it can be known that the label sample The mean and variance of the corresponding hidden vector in the k-th model are

[0126] In practice, it is necessary to control the dimension of the latent space and expect the loading matrix P k to contain data discrimination information and improve the separability between different categories of data in the latent space. Here, the hyperparameter variable α = [α1, α2,..., α d is introduced to control the dimension of the hidden vector. For this purpose, the conditional probability of P and the discrimination information containing hyperparameters is defined, and each hyperparameter controls a column of P, and its form is:

[0127]

[0128] Here, the matrix is a matrix containing data discrimination information, which is given in the F-step of the algorithm. When p ki - p' i is very large, its corresponding control parameter α i is very small, indicating that p ki deviates too much from p' i , and it should be deleted from the hidden space loading variable. In practice, the threshold of α i is set to select the effective latent variable dimension value. The objective function of the above optimization problem is transformed into the following form:

[0129]

[0130] The purpose of adding the second term is to expect P k not to deviate too much from the matrix P' containing discrimination information. Therefore, an F-step is added to the traditional EM algorithm to calculate the projection matrix containing the Fisher discrimination information of the dataset.

[0131] F-step: Update the discriminant projection matrix P'

[0132] In the F-step, the within-class scatter matrix and between-class scatter matrix are calculated using the results of the E-step.

[0133] Based on the calculated class posterior probabilities above, the within-class scatter matrix and between-class scatter matrix of the data are

[0134]

[0135]

[0136] Here, u' k is the mean of the k-th class, N = N l + N u , is the overall mean of the sample set, is the number of samples belonging to the k-th class, and the mean of the k-th class is calculated in the following form, i.e.,

[0137]

[0138] Based on the Fisher discriminant analysis principle, the projection matrix P' can be obtained from the following optimization problem, i.e.,

[0139]

[0140] Obviously, the above optimization problem is a generalized eigenvalue problem. The projection matrix P' consists of the eigenvectors corresponding to the first d largest eigenvalues of the matrix . Note that if S w is singular, generally S w + γI is used instead of S w , where the parameter γ is a very small positive number.

[0141] M-step: Update of model parameters

[0142] According to the calculation results of the above E-step and F-step, the M-step of the EM algorithm updates the model parameters by solving the optimization problem function shown in (21). For the model parameter π k , the form of its optimization problem is given according to the form of the likelihood function, i.e.,

[0143]

[0144] Solving the above optimization problem easily gives:

[0145]

[0146] The update form of Θ can be solved through the following optimization problem, i.e.,

[0147]

[0148]

[0149] Among them, A = diag{α1, α2, …, α d}。

[0150]

[0151]

[0152] Here represents the number of samples of the k-th class in the labeled sample set, represents the mathematical expectation of, and other similar symbols have the same meaning, represents the mathematical expectation of the complete likelihood function with respect to the latent variable Z u , T u , T.

[0153] Based on the above, the implementation process of the semi-supervised discriminative mixture probabilistic principal component analysis algorithm is as follows:

[0154] Algorithm 1: EFM learning algorithm of the model

[0155] (1) Calculate the mean and variance of the samples, and normalize the data; set the dimension d' of the low-dimensional latent space, and the initial latent space dimension is d = D - 1, where D is the dimension of the original input space. Let the initial value of the latent variable dimension control parameter a l = 1 (l = 1, 2, …, d'), the maximum number of iterations maxiter, the initial value of the iteration number is set to t' = 1, the algorithm stop threshold ε, and the number of modes K;

[0156] (2) Use the MPPCA algorithm to cluster all samples, and take the mean and variance of each low-dimensional subspace sample as the initial values of the prior distribution parameters of the latent variable t of this class. The MPPCA model parameters are used as the initial values of the loading matrix P k and the class prior probability . The mean and reconstruction error variance of the samples in each subspace are used as the mean u k and Ψ k (set as a diagonal matrix) initial values;

[0157] (3) E-step: Calculate the corresponding expected values and second-order statistics of the latent variables through formulas (18) - (19) using the current model parameters;

[0158] (4) F-step: Update the projection discriminant information matrix P' through formula (27);

[0159] (5) M step: Update the model parameters through formulas (29) to (33). P k , α / , u k , ψ k , k = 1, 2, …, K;

[0160] (6) If the log-likelihood function value satisfies |L(X, Y, X u |Θ t'-1 ) - L(X, Y, X u |Θ t' )| ≤ ε or the number of iterations satisfies t' ≥ max iter, then go to step (7); otherwise, t' = t' + 1 and return to step (3);

[0161] (7) Sort a l = 1 in descending order, and retain the corresponding hidden variables and model parameter values of the first d' ones u k , ψ k as the optimal model parameters.

[0162] Example 2 Multimodal Monitoring Based on Semi-Supervised Discriminative Mixture Probabilistic PCA

[0163] Since it is often difficult to separate the actual process modal data, using the hard decision method will inevitably lead to the problem of increased wrong decisions. In addition, the process data is affected by various factors and largely follows the Gaussian distribution and linear characteristics. Therefore, simply using T 2 and SPE and the hard decision method are difficult to achieve satisfactory results. To make full use of the sample information obtained from the training model, the present invention proposes a fault monitoring method based on weighted support vector data description (WSVDD) and Bayesian inference.

[0164] First, use weighted SVDD to model the hidden variables and reconstruction errors of each modal training sample respectively. Let the labeled samples and unlabeled sample sets belonging to the k-th class be The probability that each sample belongs to this class is denoted as i = 1, 2, …, N' k , N' k The number of samples in this class. For the unlabeled samples x′ in k,i its probability of belonging to the k-th class is <z i,k > = p(z i,k = 1|x' k,i , Θ), The probability of the labeled samples in k,i, then the weighted SVDD model parameters can be obtained from the following optimization problem:

[0165]

[0166]

[0167] where C is the regularization constant. The support vectors located on the spherical boundary are used to calculate the center and radius of the hypersphere in the following form:

[0168]

[0169]

[0170] where is any boundary support vector of the k-th class and h k,i” and h k,j” are elements of the set respectively, and is the number of elements in the set

[0171] Similarly, for the sample x' k,i of the k-th class, the corresponding reconstruction residual is:

[0172]

[0173] where k = 1, 2,..., K. For an optimization problem for establishing a weighted SVDD model is

[0174]

[0175]

[0176] In the formula, is the vector to be optimized in the above formula, is the coordinate in, e k,i and e k,j represent the data in the residual dataset, i, j = 1, 2,..., N' k .

[0177] According to the optimal solution obtained by solving the optimization problem, the center and radius of the hypersphere are calculated using the following formulas:

[0178]

[0179]

[0180] wherein, is any boundary support vector of the k-th class and are the support vector set and the boundary support vector set respectively, has the same meaning as formula (34'), both are elements of is the set the number of elements in The superscript T represents the transpose of a matrix or vector.

[0181] For the new input data x new , first calculate its projection onto the mean and sample reconstruction error of the latent variables in the k-th (k = 1, 2,..., K) subspace:

[0182]

[0183]

[0184] In the formula, ψ k is the noise covariance matrix of the k-th sub-model, k = 1, 2,..., K, is the loading matrix of the k-th sub-model.

[0185] To make full use of the classification information of each modality, use the Bayesian inference method to fuse the monitoring results of each local model. First calculate The posterior probability of belonging to the abnormal state is calculated by the following formula:

[0186]

[0187]

[0188] Here, 'Fau' and 'Nor' represent the abnormal and normal operating states respectively. is the prior knowledge about the abnormal operating state and the normal state in the k-th latent space. is the posterior probability of belonging to the abnormal state, is the conditional probability distribution of under the abnormal state, is the prior distribution of is the conditional probability distribution of under the normal state. The residual e new，,k The posterior probability p(Fau|e new，,k ) and the marginal distribution p(e new，,k ) are calculated by the following formulas,

[0189]

[0190]

[0191] where p(e new,k |Fau) is the conditional distribution probability of the residual e new,k under the abnormal state, p(e new,k ) is the prior distribution of the residual e new,k , and p(e new,k |Nor) is the conditional probability distribution of the residual e new,k under the normal state. are the prior probabilities of the k-th mode on the residual space for the faulty and normal states, respectively, k = 1, 2, …, K, k = 1, 2, …, K.

[0192] The conditional probabilities regarding the abnormal state and the normal state are defined as:

[0193]

[0194]

[0195] Here, v1, v2 > 0, From and give the specific forms of p(e new，,k |Fau) and p(e new，,k |Nor). In this way, it is easy to calculate the posterior probabilities of the latent variable and the residual regarding the normal state and p(Nor|e new，,k ). Given the significance level γ, the prior probabilities regarding the abnormal state and the normal state are simply set as γ and 1 - γ, respectively.

[0196] To utilize the local monitoring results corresponding to the modal latent space and the residual space, the following formula gives the fault monitoring index (BFDIs) of global Bayesian inference:

[0197]

[0198]

[0199] Here, P(k) = N' k / N, is the posterior probability distribution belonging to mode k, p(k|e new，,k ) is the posterior probability distribution of e new,k belonging to mode k, k' represents the k'-th mode, k' = 1, 2, …, K, denotes The marginal probability, p(k') represents the prior probability of mode k', p(e train，,k' ) represents the marginal probability of e train，,k' .

[0200] Global fault monitoring index confidence limit BFDI e,lim and BFDI st,lim are calculated from all training samples and are in the form of:

[0201]

[0202]

[0203] For a given significance level γ, if both global fault detection indicators are not greater than their control limits, then the new data x new is considered normal, otherwise the new data is considered abnormal data.

[0204] The implementation process of the comprehensive semi-supervised discriminant MPPCA algorithm and the fault monitoring method based on weighted support vector data description WSVDD and Bayesian inference is as follows:

[0205] Algorithm 2 Fault monitoring algorithm based on weighted support vector data description WSVDD and Bayesian inference

[0206] Offline modeling stage

[0207] S101. Input the labeled data set X L , y and the unlabeled data set X u , and call the semi-supervised discriminant MPPCA method shown in Algorithm 1 to train the labeled data set and the unlabeled data set to obtain the training results u k , ψ k , k = 1, 2,..., K, m = 1, 2,..., N u .

[0208] S102. Set the confidence value γ and the regularization constant C.

[0209] S103. Divide the unlabeled samples into the corresponding modes. The determination method for the sample belonging to the k * th mode is

[0210]

[0211] The corresponding weight is and represents The indicator variable belonging to the k-th modality, if then it means completely belongs to the k-th modality, otherwise it completely does not belong to the k-th modality; means the random variable belonging to the k-th modality, means the random variable the expectation of; Let the samples of the labeled samples and unlabeled sample set belonging to the k-th class be The probability that each sample belongs to this class is denoted as where the weight of the labeled samples is 1, i = 1, 2,..., N' k , N' k is the number of samples in this class.

[0212] S104. Solve the optimization problem shown in formula (34) based on the set of latent variables corresponding to the data set , and calculate the center and radius of the hypersphere in the latent space according to the optimal parameters and formulas (35) - (36), k = 1, 2,..., K;

[0213] The data set corresponding latent variable t' k,i , then the weighted SVDD model parameters can be obtained from the following optimization problem,

[0214]

[0215]

[0216] where, L(α k ) is a function of the function optimization variable α k , N' k represents the number of data samples in the k-th modality, α k,i , α k,j are both elements in t' k,j , t' k,i are both elements in α k , the latent variables corresponding to the data set , represents the weighting of the slack variables in the original SVDD optimization problem, i, j = 1, 2,..., N' k , C is the regularization constant; Let represent the support vector set, let be the number of support vectors of the k-th SVDD model, and the support vector set located on the spherical boundary is denoted as be the boundary support vectors of the k-th SVDD model, and calculate the center of the hypersphere using the following form and the radius

[0217]

[0218]

[0219] wherein, is an arbitrary boundary support vector of the k-th class.

[0220] S105. Calculate according to the following formula the reconstruction residual e′ corresponding to the sample x′ in the set k,i k,i ,

[0221]

[0222] In the formula, I is the identity matrix of the corresponding dimension, u k are respectively the latent variable loading matrix and the mean of the k-th sub-model, ψ k is the noise covariance matrix of the k-th sub-model, and k = 1, 2,..., K;

[0223] According to the residual data set solve the following optimization problem:

[0224]

[0225]

[0226] Then, calculate the center of the hypersphere in the residual space according to the optimal parameters and formulas (35′) - (36′) and the radius k = 1, 2,..., K,

[0227]

[0228]

[0229] wherein, is an arbitrary boundary support vector of the k-th class, are respectively the support vector set and the boundary support vector set.

[0230] Online modeling stage

[0231] S201. For the newly input data x new , calculate its expectation projected onto the k-th (k = 1, 2,..., K) latent variable and the sample reconstruction error e according to formulas (38) and (39) new，,k ​, k = 1, 2, …, K;

[0232] For the new input data x new , calculate its expectation projected onto the k-th (k = 1, 2, …, K) latent variable and the sample reconstruction error e according to equations (38) and (39); new，,k , k = 1, 2, …, K,

[0233]

[0234]

[0235] where ψ k is the noise covariance matrix of the k-th sub-model, and

[0236] is the load matrix of the k-th sub-model. S202. Solve the posterior probability belonging to the abnormal state according to equations (40) - (41), and solve the sample reconstruction error e new，,k belonging to the abnormal state according to equations similar to (40) - (41); k = 1, 2, …, K;

[0237] Solve the posterior probability belonging to the abnormal state according to equations (40) - (41), where

[0238]

[0239]

[0240] where respectively represent the prior knowledge about the abnormal operating state and the normal state in the k-th latent space, is the posterior probability belonging to the abnormal state, the conditional probability distribution of under the abnormal state, is the prior distribution of new,k . Similarly, for the residual e and is the conditional probability distribution of under the normal state; k = 1, 2, …, K.

[0241] S203. Calculate the global Bayeisan fault detection index (BFDIs) BFDI st (x train ) and BFDIe (x train );

[0242] Calculate the global Bayeisan fault monitoring index BFDI of all training samples according to the formula in step 6 above, BFDI st (x train ) and BFDI e (x train ).

[0243]

[0244]

[0245] In the formula, k = 1, 2,..., K, p(k) = N' k / N, is the posterior probability distribution belonging to mode k, p(k|e new，,k ) is the posterior probability distribution that e train,k belongs to mode k, k' represents the k'th mode, k' = 1, 2,..., K, represents 's marginal probability, p(k') represents the prior probability of mode k', p(e train，,k' ) represents the marginal probability of e train，,k' .

[0246] S204. Calculate the global fault monitoring index confidence limit BFDI from all training samples in equations (46) - (47), BFDI e,lim and BFDI st,lim ;

[0247] Calculate the global fault monitoring index confidence limit BFDI from all training samples in equations (46) - (47), BFDI e,lim and BFDI st,lim ,

[0248]

[0249]

[0250] S205. Calculate the global Bayeisan fault monitoring index BFDI of x new according to a similar formula (40) - (41), BFDI st (x new ) and BFDI e (x new );

[0251] Calculate x according to equations (44) - (45), newGlobal Bayesian Fault Detection Index BFDI st (x new ) and BFDI e (x new ),

[0252]

[0253]

[0254] wherein, P(k) = N' k / N, is the posterior probability distribution belonging to mode k, k = 1, 2,..., K.

[0255] S206. Determine whether a fault has occurred. If BFDI st (x new ) > BFDI st,lim or BFDI e (x new ) > BFDI e,lim , then a fault has occurred in the process; otherwise, the process is in a normal state.

[0256] References:

[0257] 1. Liu Yue. Industrial Process Monitoring Based on Ensemble Learning. Ph.D. Thesis, Zhejiang University, 2019.

[0258] 2. Le Yao, Zhiqiang Ge, Weiming Shao, Zhihuan Song, A Novel Scalable Semi-supervised GMM and Its Application for Multimode Process Quality Prediction with Big Data. 2018 IEEE 7th Data Driven Control and Learning Systems Conference May 25 - 27, 2018, Enshi, Hubei Province, China.

[0259] 3. Zhu Jinlin, Ge Zhiqiang, Zhihuan Song. Robust semi-supervised mixture probabilistic principal component regression model development and application to soft sensors. Journal of Process control, 32(2015)25 - 37.

[0260] 4. Bei Wang, Zhichao Li, Zhenwen Dai, Neil Lawrence, Xuefeng Yan. Data-driven mode identification and unsupervised fault detection for nonlinear multimode processes. IEEE Transactions on Industrial informatics, 16(6):3651 - 3660, 2020.

[0261] 5. Xiaogang Deng, Yongxuan Chen, Ping Wang, Yuping Cao. Soft sensor modeling for unobserved multimode nonlinear processes based on modified kernel partial least squares with latent factor clustering. IEEE Access, 8:35864 - 35872, 2020.

[0262] 6. Weiming Shao, Zhiqiang Ge, Zhihuan Song. Quality variable prediction for chemical process based on semi-supervised Dirichlet process mixture Gaussians. Chemical engineering science, 193(2019)394 - 410.

[0263] 7. Zhu Jinlin, Ge Zhiqiang, Zhihuan Song. Robust semi-supervised mixture probabilistic principal component regression model development and application to soft sensors. Journal of Process control, 32(2015)25 - 37.

[0264] 8. Reza Sharifi, Reza Langari. Nonlinear sensor fault diagnosis using mixture of probabilistic PCA models. Mechanical Systems and Signal Processing Volume 85:638 - 650, 2017.

[0265] 9. Kaixiang Peng, Kai Zhang, Bo You, et al. A quality-based nonlinear fault diagnosis framework focusing on industrial multimode batch processes. IEEE Transactions on Industrial electronics, 2016, 63(4):2615 - 2624.

[0266] 10. Kai Zhong, Dewei Ma, Min Han. Distributed dynamic process monitoring based on dynamic slow feature analysis with minimal redundancy maximal relevance. Control Engineering Practice November 2020.

[0267] 11. Jian Huang, Xu Yang, Xuefeng Yan. Slow feature analysis-independent component analysis based integrated monitoring approach for industrial processes incorporating dynamic and static characteristics. Control Engineering Practice September 2020.

[0268] 12. Wei Zhou, Bo Wang, Juan Liu, Jianyang Shi, et al. Fault diagnosis method based dynamic axis nucleation KPLS for pumping unit. IEEE Acess, 2019.

[0269] 13. Yuan Li, Dongsheng Yang. Local component based PCA model for Multimode Process Monitoring. Chinese Journal of Chemical Engineering, 27, 2020.

[0270] 14. Kong Xiangyu, Cao Zehao, Du Baiyang, Luo Jiayu. Quality-related multimode fault detection technology based on partial least squares. Control and Decision, 2019, 34(12): 2547-2557.

[0271] 15. Wei Li, Minjun Peng, Qingzhong Wang. Fault identification in PCA method during sensor condition monitoring in a nuclear power plant. Annals of Nuclear Energy, 2018.

[0272] 16. Xianchao Xiu, Ying Yang, Wanquan Liu. Laplacian regularized robust principal component analysis for process monitoring. Journal of Process Control, 2020.

[0273] 17. Rahul Raveendran, Biao Huang. Two layered mixture Bayesian probabilistic PCA for dynamic process monitoring. Journal of Process Control, 57(2017)148–163.

[0274] 18 Weipeng Lu, Xuefeng Yan. Balanced multiple weighted linear discriminant analysis and its application to visual process monitoring. Chinese Journal of Chemical Engineering, 2020.

[0275] 19. Qiu Tang, Yi Chain Xiaoyu Fang. Industrial process monitoring based on Fisher discriminant global-local preserving projection. Journal of Process Control September, 2019.

[0276] 20. Xiaogang Deng, Xuemin Tian, Chris J. Harris. Fault discriminant enhanced kernel principal component analysis incorporating prior fault information for monitoring nonlinear processes. Chemometrics and Intelligent Laboratory Systems, 2017.

Claims

1. A multi-modal monitoring method based on semi-supervised discriminant mixture probabilistic principal component analysis, characterized in that It includes two stages: offline modeling and online modeling. The specific process of the offline modeling stage is as follows: S101. Input the labeled data set X of K modalities L , y and the unlabeled data set X u , and call the semi-supervised discriminative MPPCA method to train the labeled data set and the unlabeled data set to obtain the training results k = 1, 2, …, K, m = 1, 2, …, N u ; The specific process of training the labeled dataset and the unlabeled dataset by the semi-supervised discriminant MPPCA method is as follows: (1) Calculate the mean and variance of the samples, and normalize the data; set the dimension d' of the low-dimensional latent space, and the initial latent space dimension is d = D - 1; let the initial value of the latent variable dimension control parameter a in the model l = 1 and l = 1, 2, …, d', the maximum number of iterations maxiter, the initial value of the iteration number t' = 1, the algorithm stop threshold ε, and the number of modes K; (2) Use the MPPCA algorithm to cluster all samples, and take the mean and variance of each low-dimensional subspace sample as the initial values of the prior distribution parameters of the latent variable t of this class. The MPPCA model parameters are used as the initial values of the load matrix P k , the prior probability of the class . The mean and the variance of the reconstruction error of the samples in each subspace are used as the mean u k and Ψ k initial values respectively; (3) E-step: Calculate the corresponding expected value of the latent variable and the second-order statistics using the current model parameters; (4) F-step: Update the projection discriminant information matrix P'; (5) M step: Update the model parameters P k , α l , u k , ψ k , k = 1, 2, …, K; (6) If the value of the log-likelihood function satisfies |L(X, Y, X u |Θ t'-1 ) - L(X, Y, X u |Θ t' )| ≤ ε or the number of iterations satisfies t' ≥ maxiter, then go to step (7); otherwise, set t' = t' + 1 and return to step (3); (7) For a l = 1, sort in descending order, and retain the first d' corresponding latent variables and model parameter values u k , ψ k as the optimal model parameters; S102. Set the confidence value γ and the regularization constant C; S103. Divide the unlabeled samples into the corresponding modalities. The method for determining that a sample belongs to the k*-th modality is The corresponding weight is and denotes the indicator variable belonging to the k-th modality. If then it denotes completely belongs to the k-th modality; otherwise, it completely does not belong to the k-th modality; denotes the random variable belonging to the k-th modality, denotes the expectation of the random variable ; Let the samples of the labeled samples and the unlabeled sample set belonging to the k-th class be The probability that each sample belongs to this class is denoted as where the weight of the labeled samples is 1, i = 1, 2, …, N', k , N' k is the number of samples in this class; S104. Solve the optimization problem shown based on the dataset and the corresponding set of latent variables, and calculate the center of the hypersphere in the latent space according to the optimal parameters and radius k = 1, 2, …, K; S105. Calculate according to the following formula the reconstruction residual e' corresponding to the sample x' in the set k,i k,i ,​ where \(I\) is the identity matrix of the corresponding dimension, u k are the latent variable loading matrix and mean of the \(k\)-th sub-model, respectively, ψ k is the noise covariance matrix of the \(k\)-th sub-model, and \(k = 1, 2, \ldots, K\); According to the residual data set Solve the optimization problem, and then calculate the center of the hypersphere in the residual space according to the optimal parameters and the radius k = 1, 2, …, K; The specific process of the online modeling stage is as follows: S201. For the newly input data x new , calculate the expectation of its projection onto the k-th (k = 1, 2, …, K) latent variable and the sample reconstruction error e new，,k , where k = 1, 2, …, K; S202, Solve The posterior probability belonging to the abnormal state, and solve the sample reconstruction error e new，,k The posterior probability belonging to the abnormal state, k = 1, 2, …, K; S203. Calculate the global Bayesian fault detection index BFDI for all training samples st (x train ) and BFDI e (x train ); S204. Calculate the global fault monitoring index confidence limit BFDI from all training samples e,lim and BFDI st,lim ; S205. Calculate x new The global Bayesian fault detection index BFDI of st (x new ) and BFDI e (x new ); S206. Determine whether a failure has occurred. If BFDI st (x new ) > BFDI st,lim or BFDI e (x new ) > BFDI e,lim , the process has failed; otherwise, the process is in a normal state.

2. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 1, wherein In the step S104, the data set corresponding latent variable t' k,i , then the weighted SVDD model parameters can be obtained from the following optimization problem Among them, \(L(\alpha k )\) is a function with respect to the function optimization variable \(\alpha k \), \(N'\) k represents the number of data samples of the \(k\)-th mode, \(\alpha k,i \), \(\alpha k,j are both elements in \(\alpha k \), \(t'\) k,j \), \(t'\) k,i are both the latent variables corresponding to the data set ; \(\omega represents the weighting of the slack variable in the original SVDD optimization problem, and \(i, j = 1, 2, \ldots, N'\) k , \(C\) is the regularization constant; let represent the support vector set, let be the number of support vectors of the \(k\)-th SVDD model, and the support vector set located on the spherical boundary is denoted as The center of the hypersphere and the radius Among them, is any boundary support vector of the k-th class and are respectively elements of the set and is the number of elements in the set .

3. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 2, wherein In the step S105, according to the residual data set solve the optimization problem similar to that shown in formula (34), In the formula, is the vector to be optimized in formula (34'), is coordinates in, e' k,i and e' k,j represent data in the residual dataset, i, j = 1, 2, …, N' k ; Then, according to the optimal parameters and the center of the hypersphere in the residual space is calculated according to equations (35′) to (36′) and the radius k = 1, 2, …, K Among them, is any boundary support vector of the k-th class and are the support vector set and the boundary support vector set respectively, both are elements of is the number of elements in the set ​ 4. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 3, characterized in that, In the step S201, for the newly input data x new , calculate its expectation projected onto the latent variables of the k-th (k = 1, 2, …, K) subspace according to equations (38) and (39) and the sample reconstruction error e new，,k , k = 1, 2, …, K wherein, ψ k is the noise covariance matrix of the k-th sub-model, k = 1, 2, …, K, is the load matrix of the k-th sub-model.

5. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 4, wherein In the step S202, solve according to the formulas (40) to (41). The posterior probability belonging to the abnormal state In the formula, respectively represent the prior knowledge about the abnormal operating state and the normal state in the k-th latent space, is the posterior probability belonging to the abnormal state, under the abnormal state the conditional probability distribution of is the prior distribution of under the normal state the conditional probability distribution of For the residual e new,k , the posterior probability p(Fau|e new,k ) of it belonging to the abnormal state is calculated by the following formula where p(e new,k |Fau) is the conditional distribution probability of the residual e new,k under abnormal conditions, p(e new,k ) is the prior distribution of the residual e new,k , and p(e new,k |Nor) is the conditional probability distribution of the residual e new,k under normal conditions. are the prior probabilities of the k-th mode on the residual space for the faulty and normal states, respectively, where k = 1, 2, …, K.

6. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 5, characterized in that, In the said step S203, calculate the global Bayeisan fault monitoring index BFDI of all training samples st (x train ) and BFDI e (x train ). wherein, k = 1, 2, …, K, p(k) = N' k / N, is the posterior probability distribution belonging to mode k, p(k|e new，,k 0 is e new,k the posterior probability distribution belonging to mode k, k' represents the k'-th mode, k' = 1, 2, …, K, represents the marginal probability of, p(k') represents the prior probability of mode k', p(e train，,k′ ) represents the marginal probability of e train，,k′ .

7. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 6, characterized in that In the step S204, the global fault monitoring index confidence limits BFDI are calculated from all the training samples of formulas (46) to (47). e,lim and BFDI st,lim , 8. The multimodal monitoring method of semi-supervised discriminant hybrid probabilistic principal component analysis according to claim 7, wherein In the step S205, calculate x according to the formulas (44) to (45). new The global Bayesian fault detection index BFDI st (x new ) and BFDI e (x new ). In the formula, For is the posterior probability distribution belonging to modality k, p(k|e new，,k ), where e new,k is the posterior probability distribution belonging to modality k, and k' represents the k'-th modality. denotes the marginal probability of, p(k') represents the prior probability of modality k', and p(e new，,k' ) represents the marginal probability of e new，,k' , where k' = 1, 2, …, K.

Images