Industrial process robust identification method based on RPCA-EM algorithm

Abstract

The invention discloses an industrial process robust identification method based on an RPCA-EM algorithm, and belongs to the technical field of system identification and parameter estimation. According to the system, the structure of a multi-model system is represented by adopting a linear parameter change model; in order to process abnormal values in collected data, a low-rank matrix recovery observation model based on a robust principal component analysis algorithm is established. A high-dimensional information matrix is constructed by using multiple batches of measurement data. Although the matrix is usually high-dimensional and low-rank, abnormal values can cause the matrix to become high-dimensional and high-rank. By applying a robust principal component analysis algorithm, an information matrix is recovered to a low-rank form, so that pure data is separated. In the process, a batch of separated pure data is used as an information matrix for parameter identification. The estimation problem of model parameters is deduced under the framework of an expectation maximization algorithm.

Description

Technical Field

[0001] The present invention relates to a robust identification method for industrial processes based on the RPCA-EM algorithm, and belongs to the technical field of system identification and parameter estimation. Background Art

[0002] Signal acquisition is crucial for achieving accurate modeling and control of industrial processes. However, due to data transmission errors, sensor failures, and other factors, the collected data may be contaminated by outliers. To minimize or eliminate the impact of these outliers on control and modeling accuracy, various signal processing algorithms have been developed, including direct elimination, filtering, mean compensation, and statistical feature estimation. Among them, the core idea of direct elimination is to directly delete or replace data points that significantly deviate from the normal range based on certain statistical rules and thresholds. Although this method is simple and efficient, it may mistakenly delete normal data and cannot handle local outliers. The core principle of filtering to separate outliers is to smooth the data through a filter to suppress the influence of outliers. Although it can effectively suppress noise and retain the trend of the data, this method is sensitive to the window size and has poor effects on non-stationary data and high-frequency outliers. Mean compensation replaces the identified outliers by using the mean and standard deviation of the data set, but this method is highly dependent on the mean, and the replaced outliers may cause loss of important information. The core idea of statistical feature estimation is to calculate the statistical features of the data, set thresholds based on the statistical features to judge outliers, and for the identified outliers, options such as deletion, replacement, or retention can be selected. Compared with the above other methods, statistical feature estimation has high computational efficiency, strong adaptability, can retain the overall distribution of the data, and has good accuracy and flexibility.

[0003] In fact, the most ideal method to reduce or eliminate the impact of outliers on control and modeling accuracy is to separate the outliers from the collected signals, which can avoid the impact caused by directly removing or estimating them and will not damage the integrity of the measurement data. To accurately find outliers from the collected signals, methods describing the statistical characteristics of outliers are usually adopted. By assuming that normal data has a certain pattern and outliers deviate from this pattern, data points that deviate from the pattern are identified through statistical analysis, and the data at this point is the outlier.

[0004] The current methods for describing the statistical characteristics of outliers usually use the t-distribution and the Laplace distribution. The t-distribution and the Laplace distribution focus on fitting statistical features to reduce the impact of outliers. The t-distribution assumes that the data comes from a t-distribution, calculates the probability of outliers by estimating the degrees of freedom, is suitable for small samples, and has good flexibility. However, the calculation of the t-distribution is relatively more complex, depends on the choice of degrees of freedom, and is ineffective for extreme outliers. The Laplace distribution identifies outliers by transforming the data into the frequency domain and analyzing the characteristics of the frequency domain. Based on its spike characteristics and heavy-tail characteristics, the data has a higher concentration at the central position and can better describe extreme values or outliers. This method is suitable for time series data and can capture transient outliers. However, due to the complexity of the Laplace transform itself and the large amount of calculation, it introduces additional complexity and has high requirements for data stationarity.

[0005] Therefore, there is an urgent need for a method to eliminate the influence of outliers on the control and modeling accuracy without compromising the integrity of the measurement data. Summary of the Invention

[0006] To solve the above problems existing in the current prior art, the present invention provides a robust identification method for industrial processes based on the RPCA-EM algorithm. Robust principal component analysis (RPCA) is a sparse low-rank matrix recovery algorithm for solving convex optimization problems. This algorithm decomposes a given matrix into a low-rank matrix and a sparse noise matrix, effectively separating information data from outliers. However, in the field of system identification, the collected data often only has the dimension of a single information vector and lacks direct high-dimensional features. To solve this problem, multiple batches of input and output data are collected to form information vectors and combined into a high-dimensional N×N matrix (N represents the data length). This constructed high-dimensional matrix contains pure information data and sparse outliers. Based on the RPCA algorithm, this high-dimensional matrix is divided into two parts to obtain pure information data. The linear parameter varying (LPV) model is used to accurately describe nonlinear or time-varying systems by combining linear structures with time-varying parameters, and the expectation maximization (EM) principle is used to establish a parameter identification algorithm to optimize the parameter estimation process and enhance the robustness of the model. An effective identification of a multi-model system containing outliers is achieved through a robust identification method based on the RPCA-EM algorithm.

[0007] The technical solution is as follows:

[0008] Step 1: Analyze the characteristics of multi-stage dynamic changes with outliers in the industrial process through the linear parameter varying (LPV) model:

[0009] The model order in the linear parameter varying (LPV) model is n a and n b The local model expression is:

[0010]

[0011] Among them, k represents the sampling time, m represents the m-th local model, T represents the transpose, and v k represents Gaussian white noise with a mean of zero and a variance of σ 2 , the information vector φ k and the parameter vector θ m are expressed as:

[0012]

[0013] Among them, represents the real vector space with (n a +n b ) rows and 1 column. u k and y k are the input and output data collected at time k respectively. a m,i and b m,j are the system parameters of the i-th and j-th in the m-th local model respectively, where i = {1, 2,..., n a}, j = {1, 2,..., n b}.

[0014] Based on the local model, considering the global characteristics, a weighted fusion strategy is introduced to synthesize the outputs of all local models; consider the following linear parameter varying model:

[0015]

[0016] Among them, w k is the scheduling variable, is the weight function, M represents the total number of local models, all local models are weighted, and y k is the global output.

[0017] The migration of the operating point usually changes smoothly. Therefore, a normalized exponential function is selected as the weight factor to achieve a smooth transition between local models. The weight function is expressed as:

[0018]

[0019] Among them, {o m |m = 1,..., M} is the effective width of the local model, and the value range is [o min , o max . The effective width o m reflects the effective working range of the m-th local model. Let H = {H m |m = 1,..., M} be M different operating points of the industrial process.

[0020] Global output y k It includes the pure output and the interference from outliers, and the expression is:

[0021]

[0022] where y pk represents the pure output data, represents the outliers.

[0023] Step 2: Use the RPCA algorithm to transform the global information matrix of the multi-model system with outliers into a low-rank matrix and a sparse noise matrix through a convex optimization problem;

[0024] Step 2.1: Implement the decomposition of the matrix through a convex optimization problem;

[0025] Use the robust principal component analysis algorithm to transform the pure output and the output containing outliers in the output data into a low-rank matrix Z and a sparse noise matrix E through a convex optimization problem; for the information matrix X ∈ R (e×f) which contains f sampling observations, and each observation contains e variables, and a target convex optimization function problem needs to be solved:

[0026] min rank(Z)+γ‖E‖0(1)

[0027] s.t.X = Z + E

[0028] where Z ∈ R (e×f) is the low-rank matrix, and E ∈ R (e×f) is the sparse noise matrix containing outliers. ‖·‖0 is the l0 norm of the matrix, and γ is the balance factor parameter, which is usually initialized as:

[0029]

[0030] However, the problem in formula (1) is NP-hard (non-deterministic polynomial hard problem), and this problem is solved by relaxing the objective function of the optimization problem. The nuclear norm is the convex hull of the matrix rank function and is used to approximately replace the rank; the l1 norm is the convex hull of the l0 norm and is used to approximate the l0 norm.

[0031] Therefore, the objective function in formula (1) is converted into a convex optimization function, which is expressed as follows:

[0032] min‖Z‖ * +γ‖E‖1(2)

[0033] s.t.X = Z + E

[0034] where ‖·‖ *denotes the nuclear norm, and ‖·‖1 denotes the l1 norm; thus, the information matrix X is separated into a low-rank matrix Z containing pure data and a sparse noise matrix E containing outliers.

[0035] Step 2.2: Introduce the soft shrinkage operator and the singular value shrinkage operator to calculate the low-rank matrix Z and the sparse noise matrix E;

[0036] The present invention forms an information vector {φ k ,|k = 1, 2, …, N} by collecting input and output data of length N. Since matrix factorization based on RPCA requires a high-dimensional matrix, multiple batches of information vectors are collected to form an N×N high-dimensional information matrix X containing sparse outliers, where {d = 1, 2, …, D} represents the batches of sampled data.

[0037] Let X be:

[0038]

[0039] where, n a and n b are the model orders, N is the length of the input and output data to form the information vector, and D is the maximum value of the batches of sampled data.

[0040] Introducing the soft shrinkage operator and the singular value shrinkage operator is required to calculate the low-rank matrix Z and the sparse noise matrix E.

[0041] Define: Q = (q1, q2, …, q f ), Q ∈ R e×f , Q is a real matrix,

[0042] Soft shrinkage operator:

[0043]

[0044] If there is an optimization problem:

[0045]

[0046] Then there is a soft shrinkage operator that can solve this problem;

[0047] Singular value shrinkage operator:

[0048]

[0049] If there is an optimization problem:

[0050]

[0051] Then there is a singular value shrinkage operator that can solve this problem;

[0052] Among them, ‖·‖1 represents the l1 norm, ∈ represents the penalty factor, and ‖·‖ F represents the Frobenius norm, and sign(q i ) represents the sign function of q i , q i ∈Q, γ represents the balance factor; ‖·‖ * represents the nuclear norm, is obtained through the singular value decomposition of matrix X. U and V are the orthogonal matrices composed of left singular vectors and right singular vectors respectively, and T represents the transpose.

[0053] Step 2.3: Use inexact augmented Lagrangian multiplier (IALM) to solve the convex optimization problem in Step 2.1; define the Lagrangian equation expression as:

[0054]

[0055] where β is the penalty parameter in IALM, G is the Lagrangian multiplier matrix used to integrate the constraint objective into the objective function, and <·> represents the matrix inner product operation.

[0056] Optimize the low-rank matrix Z and the sparse noise matrix E respectively,

[0057] For:

[0058] Z min = argmin z L(Z, E, g, β)

[0059] The optimization result is:

[0060]

[0061] For:

[0062] E min = argmin E L(Z, E, G, β)

[0063] The optimization result is:

[0064]

[0065] The optimization iteration formulas for the sparse noise matrix E and the low-rank matrix Z are obtained by the soft shrinkage operator and the singular value operator respectively as:

[0066]

[0067] where h represents the number of iterations, and γ represents the balance factor.

[0068] Step 3: Construct an Expectation-Maximization (EM) framework to update the parameters of the Linear Parameter-Varying (LPV) model;

[0069] The above algorithm decomposes the process data X into a low-rank matrix Z and a sparse noise matrix E. The low-rank matrix Z contains pure multi-batch data, and the sparse noise matrix E contains outlier data. Select a batch of data in Z as the information matrix, and use the Expectation-Maximization (EM) principle to establish an Expectation-Maximization algorithm to update the parameters of the linear parameter-varying model. The iterative calculation formula for model parameter estimation is:

[0070] Step 3.1: Expectation step (E-step): Calculate the Q-function, that is, the expectation of the latent variable under the current parameter estimation;

[0071] Construct the Q-function:

[0072]

[0073] where S is the loop variable for iterative calculation, obeys the Gaussian distribution, is the normalized weighting function, is the posterior output of the weight;

[0074]

[0075] Step 3.2: Maximization step (M-step); Maximize the Q-function to update the model parameters by taking the partial derivatives of the Q-function with respect to each unknown variable and setting them to zero

[0076] The specific steps are as follows:

[0077]

[0078]

[0079] The iterative formula for variance is

[0080]

[0081] Step 3.3: Iterative update; Repeat the expectation step and the maximization step until the parameter estimation converges.

[0082] Step 4: Use the model obtained by identifying the parameters obtained in Step 3 to predict the output data of the industrial process.

[0083] Input the data to be predicted into the model obtained by identifying the parameters obtained in Step 3 to predict the output data of the industrial process.

[0084] The beneficial effects of the present invention are:

[0085] The present invention proposes a robust identification method for industrial processes based on the RPCA-EM algorithm. The data information matrix is decomposed into a low-rank matrix and a sparse noise matrix through the Robust Principal Component Analysis (RPCA) algorithm, and then the parameter estimation of the model is iteratively calculated through the Expectation-Maximization (EM) principle. By processing the data through low-rank matrix decomposition, outliers can be accurately separated from the measurement data, avoiding the interference of outliers on model identification and parameter estimation, and enabling the algorithm to maintain high robustness when facing outliers in the data. At the same time, the algorithm can also achieve a separation efficiency of 99.8% under different outlier ratios. Even when the outlier ratio is relatively high, it can accurately recover the low-rank characteristics of the multi-batch information matrix, ensuring the stability and reliability of the algorithm in various abnormal situations. Compared with other algorithms for processing outlier signals, the algorithm proposed in the present invention significantly reduces the relative error of parameter estimation; and it can accurately estimate the time-varying delay, and its estimation accuracy only drops by about 1% to 2% under different outlier ratios, which is far better than other methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0087] Figure 1 It is a schematic diagram of scheduling variable data and model parameters provided in the second embodiment of the present invention;

[0088] Figure 2 It is a schematic diagram of the distribution of an information matrix containing 5% outliers provided in the second embodiment of the present invention;

[0089] Figure 3 It is a schematic diagram of the normalized weights of a local model under the condition of 5% outliers provided in the second embodiment of the present invention;

[0090] Figure 4 It is a schematic diagram of the input and output data of a continuous stirred tank reactor in the presence of outliers provided in the second embodiment of the present invention;

[0091] Figure 5 It is a schematic diagram of the input and output data of a continuous stirred tank reactor after separating outliers provided in the second embodiment of the present invention;

[0092] Figure 6 It is a schematic diagram of the comparison of predicted value estimates under different processing algorithms with 15% outliers provided in the second embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0093] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0094] Embodiment 1

[0095] This embodiment provides a robust identification method for industrial processes based on the RPCA-EM algorithm, which separates the outliers in industrial process data and uses the separated data to identify a linear parameter varying model. Finally, the output data of the industrial process is predicted using the model with updated parameters. The technical solution adopted is as follows:

[0096] Step 1: Analyze the dynamic process of a multi-model system with outliers through a linear parameter varying (LPV) model.

[0097] In the linear parameter varying (LPV) model, the model order is n a and n b The local model expression is:

[0098]

[0099] where k represents the sampling time, m represents the m-th local model, T represents the transpose, v k represents Gaussian white noise with a mean of zero and a variance of σ 2 , the information vector φ k and the parameter vector θ m The expressions are:

[0100]

[0101] where T represents the transpose, u k and y k are the collected input and output data respectively, a m,i , b m,j are the system parameters of the i-th and j-th in the m-th local model respectively, i = {1, 2,..., n a}, j = {1, 2,..., n b}.

[0102] Based on the local model, considering the global characteristics, a weighted fusion strategy is introduced to synthesize the outputs of all local models. Consider the following linear parameter varying model:

[0103]

[0104] where w k is the scheduling variable, is the weight function, which weights all local models, and y k is the global output.

[0105] The migration of the working point usually changes smoothly. Therefore, a normalized exponential function is selected as the weight factor to achieve a smooth transition between local models. The weight function is expressed as:

[0106]

[0107] where {o m |m = 1, …, M} is the effective width of the local model, and the value range is [o min , o max . The effective width o m reflects the effective working range of the m-th local model. Let H = {H m} m＝1,…,M be M different working points of the industrial process.

[0108] The global output y k of the model contains the pure output and the interference from outliers, and the expression is:

[0109]

[0110] where y pk represents the information data, represents the outliers.

[0111] Step 2: Decomposition of the information matrix based on the robust principal component analysis algorithm;

[0112] Step 2.1: Use the RPCA algorithm to transform the global output of the multi-model system with outliers into a low-rank matrix and a sparse noise matrix through a convex optimization problem;

[0113] Use the robust principal component analysis algorithm to transform the pure output and the output containing outliers in the output data into a low-rank matrix Z and a sparse noise matrix E through a convex optimization problem. For the information matrix X ∈ R (e×f) with f sampling observations, and each observation contains e variables. To represent the information matrix X as the sum of a low-rank matrix and a sparse noise matrix, a target convex optimization function problem needs to be solved:

[0114] min rank(Z) + γ‖E‖0(1)

[0115] s.t. X = Z + E

[0116] where Z ∈ R (e×f) is the low-rank matrix, and E ∈ R (e×f) is the sparse noise matrix containing outliers. ‖·‖0 is the l0 norm of the matrix, and γ is the balance factor parameter, usually initialized as:

[0117]

[0118] However, the problem in formula (1) is NP-hard (Non-deterministic Polynomial hard problem), and this problem is solved by relaxing the objective function of the optimization problem. The nuclear norm is the convex envelope of the matrix rank function and is used to approximately replace the rank; the l1 norm is the convex envelope of the l0 norm and is used to approximate the l0 norm.

[0119] Therefore, the objective function in formula (1) is converted into a convex optimization function, which is expressed as follows:

[0120] min‖Z‖ * +γ‖E‖1(2)

[0121] s.t. X = Z + E

[0122] where, ‖·‖ * represents the nuclear norm, ‖·‖1 represents the l1 norm, so as to separate the information matrix X into a low-rank matrix Z containing pure data and a sparse noise matrix E containing outliers.

[0123] The global output y is separated into information data and outliers through the Robust Principal Component Analysis (RPCA) algorithm k Not only is the direct deletion of data avoided, but valuable information is also retained.

[0124] Step 2.2: Introduce a soft shrinkage operator and a singular value shrinkage operator to calculate the low-rank matrix Z and the sparse noise matrix E;

[0125] The present invention forms an information vector {φ k , |k = 1, 2, …, N} by collecting input and output data of length N. Since the matrix decomposition based on RPCA requires a high-dimensional matrix, multiple batches of information vectors are collected to form an N×N high-dimensional information matrix X containing sparse outliers, where, {d = 1, 2, …, D}, representing the batches of sampled data.

[0126] Let the information matrix X be:

[0127]

[0128] where, n a and n b are the model orders, E is the length of the input and output data to form the information vector, and D is the maximum value of the batches of sampled data.

[0129] Introducing a soft shrinkage operator and a singular value shrinkage operator is required to calculate the low-rank matrix Z and the sparse noise matrix E.

[0130] Define: Q = (q1, q2, …, q f ), Q ∈ Re×f , Q is a real matrix,

[0131] Soft shrinkage operator:

[0132]

[0133] If there is an optimization problem:

[0134]

[0135] Then there is a soft shrinkage operator that can solve this problem;

[0136] Singular value shrinkage operator:

[0137]

[0138] If there is an optimization problem:

[0139]

[0140] Then there is a singular value shrinkage operator that can solve this problem;

[0141] Among them, ‖·‖1 represents the l1 norm, ∈ represents the penalty factor, ‖·‖ F represents the Frobenius norm, sign(q i ) represents the sign function of q i where q i ∈Q, γ represents the balance factor; ‖·‖ * represents the nuclear norm, Obtained by the singular value decomposition of matrix X, U and V are the orthogonal matrices composed of left singular vectors and right singular vectors respectively, and T represents the transpose.

[0142] Step 2.3: Use inexact augmented Lagrangian multiplier (IALM) to solve the convex optimization problem in Step 2.1; define the Lagrangian equation expression as:

[0143]

[0144] Among them, β is the penalty parameter in IALM, G is the Lagrangian multiplier matrix used to integrate the constraint objective into the objective function, and <·> represents the matrix inner product operation.

[0145] Optimize the low-rank matrix Z and the sparse noise matrix E respectively,

[0146] For:

[0147] Z min = argmin z L(Z, E, G, β)

[0148] Optimized to obtain:

[0149]

[0150] For:

[0151] E min = argmin E L(Z, E, G, β)

[0152] Optimized to obtain:

[0153]

[0154] The optimized iteration formulas for the sparse noise matrix E and the low-rank matrix Z are obtained by the soft shrinkage operator and the singular value operator respectively as follows:

[0155]

[0156] where h represents the number of iterations and γ represents the balance factor.

[0157] Step 3: Construct an Expectation-Maximization (EM) framework to update the model parameters;

[0158] The above algorithm decomposes the process data X into a low-rank matrix Z containing pure multi-batch data and a sparse noise matrix E containing outlier data. Select a batch of data in Z as the information matrix, and use the Expectation-Maximization (EM) principle to construct an Expectation-Maximization algorithm to update the parameters of the parameter change model. The iterative calculation formula for model parameter estimation is:

[0159] Step 3.1: Expectation step (E-step): Calculate the Q-function, that is, the expectation of the latent variable under the current parameter estimation;

[0160] Construct the Q-function:

[0161]

[0162] where S is the loop variable for iterative calculation, obeys the Gaussian distribution, is the normalized weighting function, is the posterior output of the weight,

[0163]

[0164] Step 3.2: Maximization step (M-step); Maximize the Q-function to update the model parameters by taking the partial derivatives of the Q-function with respect to each unknown variable and setting them to zero

[0165] The specific steps are as follows:

[0166]

[0167] The iterative formula for variance is

[0168]

[0169]

[0170] Step 3.3: Iterative update; repeat the expectation step and the maximization step until the parameter estimation converges.

[0171] Step 4: Verification and evaluation;

[0172] Separate the outliers that appear in the industrial process using the optimized parameters obtained in Step 3.

[0173] Example 2

[0174] First, explain the basic concepts involved in this example. In the chemical industry, a Continuous Stirred Tank Reactor (CSTR), as a key production device, ensures uniform mixing of the feed materials in the reactor by virtue of its continuous stirring function, thereby efficiently carrying out chemical reactions to produce products that meet quality standards. Based on the conservation of mass and heat, the following continuous-time nonlinear differential equations are derived to describe its process dynamics:

[0175]

[0176] where, C At is the concentration of compound A in the output at time t, C A0t is the concentration of compound A in the feed at time t, T t is the temperature of the reactor at time t, q t is the feed flow rate at time t, T 0t is the feed temperature at time t, q ct is the cooling water flow rate at time t, T e0t represents the cooling water temperature at time t, and t represents time; other process variables and their steady-state values are shown in Table 1. In the present invention, the output temperature and the cold water flow rate are selected as the output and input data respectively for modeling. Since the cooling water flow rate has an important influence on the process dynamics, it is selected as the scheduling variable h. The operating point is set as H = [97, 100, 103].

[0177] Table 1 CSTR system parameters and their steady-state values

[0178]

[0179]

[0180] This embodiment provides a robust identification method for the process data of a stirred tank reactor based on the RPCA-EM algorithm. This method separates the outliers generated during the operation of the stirred tank reactor based on the robust identification method for industrial process data based on the RPCA-EM algorithm provided in Embodiment 1, then updates the linear change model parameters using the clean data after separating the outliers, and finally predicts the output of the stirred tank reactor through the updated parameters. Taking the LPV system shown below as an example:

[0181]

[0182] Let the sampling data length N be 1300, and the input data be u k Take a generalized binary signal, v k be Gaussian noise with a mean of zero and a variance σ 2 = 0.1 2 and set the changes of three operating points to be w1 = 0.1, w2 = 0.5, and w3 = 0.9 respectively. The variation of the model parameters with the scheduling variable is as shown in Figure 1 As shown, as the operating point changes, the parameters of the model change with the change of the scheduling variable, showing the characteristics of linear parameter variation.

[0183] In this simulation, 325 batches of data are collected to form an information matrix. Each batch of data randomly contains 5% outliers. An information vector matrix containing 5% outliers is as shown in Figure 2 As shown, it can be seen that the outliers deviate significantly from the normal sampling points. The normalized weights of each local model in the case of containing 5% outliers are as shown in Figure 3 As shown, it can be seen that at different operating points, the weights of the model are different, and each stage corresponds to a local model as the main contributing model.

[0184] After identification, the model weights under each working trajectory can be clearly identified, indicating that the selection of the main model is different at different operating points, and thus the models that play a role in each stage are different.

[0185] The outlier separation method for the stirred tank reactor based on RPCA-EM separates the outliers from the information vector. The input and output data of the stirred tank reactor before and after separation are as shown in Figure 4 and Figure 5 As shown; Figure 4 represents the output data and input data containing outliers, Figure 5 represents the input and output data after the outliers are separated. It can be seen that the method proposed in the present invention has good outlier separation efficiency.

[0186] To evaluate the effectiveness of the algorithm in handling outliers in a multi-model system, a comparative analysis of the parameter estimation ability under different degrees of outlier contamination was conducted. Specifically, the performance of the system was evaluated when the outliers accounted for 5%, 10%, and 15% of the data, as shown in Table 2.

[0187] Table 2 Various indicators of the RPCA algorithm for separating outliers

[0188]

[0189] Among them, S represents the number of iterations required to reach the convergence condition, Δ represents the proportion of the separated outliers, and ‖E‖0 represents the l0 norm of matrix E.

[0190] As shown in Table 1, when the outlier ratios are 5%, 10%, and 15%, the number of outliers in the information matrix of dimension 1300×1300 is 42,250, 85,400, and 126,750 respectively. The RPCA algorithm based on IALM iteratively calculates the separation under these different outlier ratios, and the separation efficiency in each case is 99.8%. Despite the continuous increase in the proportion of outliers, the rank of the multi-batch information matrix always returns to 4, indicating that the method proposed in the present invention can accurately identify the number of columns in the original system information matrix. This proves the robustness of the separation efficiency of the algorithm under different outlier ratios.

[0191] To intuitively evaluate the performance of the algorithm in practical applications, the root mean square error (RMSE) was used to measure the difference between the observed value and the estimated value. The calculation method of RMSE is as follows:

[0192]

[0193] Among them, N represents the number of collected data, y k represents the output value, represents the estimated output value; the results of parameter estimation after 300 iterations under different outlier ratios are shown in Table 3;

[0194] Table 3 Parameter estimation after 300 iterations under different outlier ratios

[0195]

[0196]

[0197] From the results in Table 2, it is observed that as the outlier ratio increases, the deviation between the estimated results of different parameters at different ratios and the true values is small, and the relative error δ of parameter estimation only increases slightly, about 0.01%. The results show that the RPCA-EM algorithm proposed in the present invention is not affected by the outlier ratio and exhibits high robustness.

[0198] To further verify the effectiveness of the algorithm proposed in the present invention, different signal processing algorithms were used to process the information vector containing 15% outliers, and the processing results are as follows Figure 6 shown. The errors of the unprocessed output (UP) and the result of directly deleting outliers (RO) relative to the true value are relatively large, while the method proposed in the present invention can well fit the output result of the true value, has a strong prediction ability for the true value output, and has high accuracy.

[0199] Table 4 Comparison of parameter estimation with 15% outliers under different processing algorithms

[0200]

[0201] As can be seen from the results in Table 4, both the unprocessed output (UP) or directly deleting outliers (RO) will cause a significant increase in the relative error δ of parameter estimation; while using the method proposed in the present invention, the relative error of parameter estimation is significantly reduced, from 71.79% of the unprocessed output to 2.79%; the RMSE of the predicted output is significantly reduced.

[0202] In this embodiment, the information matrix is decomposed into a low-rank matrix containing the main component data of the system and a sparse noise matrix containing outliers through the RPCA algorithm, so as to separate the outliers from the normal data, avoid the interference of outliers on model recognition and parameter estimation, and enable the algorithm to still maintain high robustness when facing outliers in the data; and this algorithm can also maintain high separation efficiency and accuracy under different outlier ratios, and has strong stability; compared with other methods, it significantly reduces the relative error of parameter estimation.

[0203] Some steps in the embodiments of the present invention can be implemented by software, and the corresponding software program can be stored in a readable storage medium, such as an optical disc or a hard disk, etc.

[0204] The above are only the preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A robust identification method for industrial processes based on the RPCA-EM algorithm, characterized in that The method includes: Step 1: Analyze the multi-stage change characteristics with outliers in the industrial process through a linear parameter varying model; Step 2: Use the RPCA algorithm to transform the global information matrix of the multi-model system with outliers into a low-rank matrix and a sparse noise matrix through a convex optimization problem; Step 3: Update the parameters of the linear parameter varying model using the data of the low-rank matrix obtained in Step 2 through the expectation maximization framework; Step 4: Input the data to be predicted into the model with updated parameters in Step 3 to predict the output of the industrial process data.

2. The method according to claim 1, wherein The said Step 1 includes: The model order in the linear parameter varying model is n a and n b The local model expression is as follows: where k represents the sampling time, m represents the m-th local model, T represents the transpose, and v k represents Gaussian white noise with zero mean and variance σ 2 , the information vector φ k and the parameter vector θ m are expressed as: Among them, represents a real vector space with (n a + n b ) rows and 1 column. u k and y k are the input data and output data in the industrial process collected at time k respectively. a m,i and b m,j are the i-th and j-th system parameters in the m-th local model respectively, where i = {1, 2, …, n a}, j = {1, 2, …, n b}; Based on the local model, considering the global features, introducing a weighted fusion strategy to synthesize the outputs of each local model; considering the following linear parameter varying model: Among them, w k is a scheduling variable, is a weight function, y m,k represents the output data of the m-th model, M represents the total number of local models, and all local models are weighted, y k is the global output; Weight function The expression is as follows: Among them, {o m | m = 1, …, M} is the effective width of the local model, and the value range is [o min , o max , and the effective width o m reflects the effective working range of the m-th local model. Let H = {H m} m＝1,…,M be M different operating points of the industrial process.

3. The method according to claim 2, characterized in that The said Step 2 includes: Step 2.1: Achieve the decomposition of the matrix through a convex optimization problem; For an information matrix \(X\in\mathbb{R}\) containing outliers or sparse noise (e×f) with \(f\) sampled observations, where each observation contains \(e\) variables, to represent the information matrix \(X\) as the sum of a low-rank matrix and a sparse noise matrix, a target convex optimization problem needs to be solved: min rank(Z)+γ‖E‖0 (1) s.t.X=Z+E where Z ∈ R (e×f) is a low-rank matrix, E ∈ R (e×f) is a sparse noise matrix containing outliers; ‖·‖0 is the l0 norm of matrix E, and γ is a balancing factor, initialized as: The nuclear norm is the convex envelope of the matrix rank function and is used to replace the rank; the l1 norm is the convex envelope of the l0 norm and is used to replace the l0 norm. By relaxing the objective function of the optimization problem, the objective function in formula (1) is converted into a convex optimization function, and its expression is: min‖Z‖ * +γ‖E‖1 (2) s.t.X=Z+E where, ‖·‖ * represents the nuclear norm, ‖·‖1 represents the l1 norm, so as to separate the information matrix X into a low-rank matrix Z containing pure data and a sparse noise matrix E containing outliers.

4. The method according to claim 3, characterized in that, The said Step 2 includes: Step 2.2: Introduce a soft thresholding operator and a singular value thresholding operator to calculate the low-rank matrix Z and the sparse noise matrix E; The present invention collects information vectors {φ of length N k | k = 1, 2, …, N}, and forms an N×N high-dimensional information matrix X containing sparse outliers by collecting multiple batches of information vectors where {d = 1, 2, …, D} represents the batches of sampled data; Let X be: where n a and n b are the model orders, E is the length of the input and output data to form the information vector, and D is the maximum value of the sampling data batches; Introduce a soft thresholding operator and a singular value thresholding operator to calculate the matrix Z and the sparse noise matrix E; Definition: Q = (q1, q2, …, q f ), Q ∈ R e×f , where Q is a real matrix, Soft thresholding operator: If there is an optimization problem: Then solve this optimization problem through the soft thresholding operator; Singular value thresholding operator: If there is an optimization problem: Then solve this optimization problem through the singular value thresholding operator; Among them, ‖·‖1 represents the l1 norm, ∈ represents the penalty factor, γ represents the balance factor, and ‖·‖ F represents the Frobenius norm, and sign(q i ) represents the sign function of q i , where q i ∈Q; ‖·‖ * is the nuclear norm, obtained by the singular value decomposition of matrix X, where U and V are the orthogonal matrices composed of left singular vectors and right singular vectors respectively, and T represents the transpose.

5. The method according to claim 4, wherein The said Step 2 includes: Step 2.3: Use the inexact augmented Lagrangian multiplier method to iteratively solve the convex optimization problem in 2.1; Define the Lagrangian equation expression as: where β is the penalty parameter in IALM, G is the Lagrange multiplier matrix used to integrate the constraint objective into the objective function, and <·> represents the matrix inner product operation; Optimize the low-rank matrix Z and the sparse noise matrix E respectively, For: Z min = argmin z L(Z, E, G, β) The optimized result is: For: E min = argmin E L(Z, E, G, β) The optimized result is: The optimized iteration formulas for the sparse noise matrix E and the low-rank matrix Z obtained by the soft thresholding operator and the singular value operator respectively are: where h represents the number of iterations, and γ represents the balance factor.

6. The method according to claim 5, wherein In the said Step 3: Select a batch of data in the low-rank matrix Z as the information matrix, and use the expectation maximization principle to construct an expectation maximization algorithm to update the parameters of the linear parameter varying model; The said Step 3 includes: Step 3.1: Expectation step, construct the Q-function, that is, the expectation of the latent variable under the current parameter estimation; The iterative calculation formula for the model parameter estimation is: where S is the loop variable for iterative calculation, obeys a Gaussian distribution, is the normalization weighting function, is the posterior output of the weights; Step 3.2: Maximization step; update the model parameters by obtaining the partial derivatives of the Q function with respect to each unknown variable and setting them to zero The specific steps are: The iterative formula for the variance is: Step 3.3: Iterative update; repeat the expectation step and the maximization step until the parameter estimation converges.

7. The method according to claim 6, wherein The said Step 4 includes: Predict the output data in the industrial process through the model with updated parameters in Step 3.

8. The method according to claim 7, characterized in that The said industrial process includes chemical processes, automobile manufacturing, iron and steel smelting, and food processing.

9. An outlier separation method in a stirred tank reactor based on the RPCA-EM algorithm, characterized in that, The method separates outliers generated during the operation of the stirred tank reactor based on the method according to any one of claims 1-7.

10. The method according to claim 9, wherein The input data during the operation of the stirred tank reactor is the cold water flow rate, and the output data is the output temperature.

Images