Robust recursive expectation maximization method for dynamic system parameter identification

Abstract

The invention relates to a robust recursive expectation maximization method for dynamic system parameter identification, which comprises the following steps: selecting a proper model structure according to the characteristics of a modeled object, including the input and output dimensions and order of a model; the method comprises the following steps: collecting input data and detection data of a dynamic system in a sampling period to form a group of data samples; the scaling factor obeys Gamma distribution, and according to the property of conjugate distribution, the scaling factor posterior expectation of the current sampling data is calculated; calculating the posteriori expectation of the logarithm of the scaling factor of the current sampling data according to the property of the conjugate distribution; carrying out recursive calculation on sufficient statistics corresponding to regression parameter vectors of the model, and carrying out online updating on the regression parameter vectors; performing recursive calculation on the sufficient statistical magnitude corresponding to the normal data noise variance, and performing online updating on the normal data noise variance; and performing recursive calculation on the intermediate statistics corresponding to the degree-of-freedom parameters distributed by the students t, solving an optimization equation, and performing online updating on the degree-of-freedom parameters.

Description

Technical Field

[0001] The present invention relates to a robust recursive system model parameter estimation method for detection data contaminated by outliers, and in particular to a robust recursive expectation maximization method for dynamic system parameter identification, belonging to the technical field of industrial process dynamic system model parameter identification. Background Art

[0002] The mathematical model of an industrial process is a quantitative description of the dynamics of an industrial process, reflecting the input-output relationship of a dynamic process in a deterministic way. Based on the mathematical model, we can explore the development law of a dynamic process, predict the changing trend of a dynamic process, and reverse the input conditions of the dynamic process. Therefore, the mathematical model of an industrial process plays an indispensable and important role in production optimization, process control, state estimation, and process monitoring. Obtaining a mathematical model that is accurate or within the allowable error range is a common issue faced by the production and control fields.

[0003] There are many methods for industrial process modeling, among which data-driven methods have attracted much attention in the industry because they are not constrained by precise prior knowledge. Data-driven methods are generally a modeling method between gray-box modeling and black-box modeling. Under the preset model structure, the parameters of the model are estimated and identified based on sampled data. However, due to the complexity of industrial processes, data-driven modeling methods often face the problem of missing information, such as unknown data delays and loss of data information. The lack of information often leads to a decrease in modeling accuracy or model mismatch. Therefore, the problem of missing information has become an obstacle in industrial process modeling, hindering the implementation of modeling work, and the problem of outlier data can also be classified as a problem of missing information. Outliers are usually defined as points outside the three standard deviations (i.e., 3-σ) of normal data. They are a type of abnormal data and appear randomly or concentrated in the sampled data. The causes of outliers are generally due to sensor data acquisition failures, abnormal interference in data transmission, data receiving module failures, etc. Outlier data usually have a great impact on the parameter identification method, making it impossible for the parameter identification method to converge to a stable point, or causing a large deviation in the identification model, leading to the failure of modeling.

[0004] Existing outlier processing methods include deletion method, interpolation method and Huber estimator method. These methods judge outliers based on the mathematical and statistical characteristics of existing batch data. However, all three methods have human intervention, which will affect the modeling accuracy to a certain extent. For the estimation problem under missing information, the most commonly used method at this stage is the offline expectation maximization method, also known as the EM method, which can achieve a good balance between the method's computational complexity and the ability to handle missing information. However, due to the iterative characteristics of the traditional expectation maximization method, it is required to generate corresponding expectations and estimates based on a large number of data samples. Therefore, the traditional robust expectation maximization method is not suitable for robust online estimation of industrial process dynamic models. The model of the industrial process has the characteristics of time-varying, that is, as the internal process characteristics of the industrial process change, the external characteristics of its dynamic model will change and migrate. This time-varying characteristic will cause changes in the parameters of the system dynamic model, which requires special processing. If the dynamic changes of the model are not noticed by the staff before modeling, or if no relevant design is made in the modeling method, the accuracy of the model will be affected. However, offline iterative identification methods based on batch data are difficult to respond to unknown dynamic changes in a timely manner and adjust the parameters of the identification model accordingly.

[0005] Therefore, it is necessary to develop a robust recursive expectation maximization parameter identification method to solve the problem of online modeling of industrial process dynamic systems under outlier pollution. Summary of the invention

[0006] The purpose of the present invention is to overcome the deficiencies of the prior art and provide a robust recursive expectation maximization method for dynamic system parameter identification.

[0007] The purpose of the present invention is achieved through the following technical solutions:

[0008] The robust recursive expectation maximization method for dynamic system parameter identification is characterized by comprising the following steps:

[0009] Step 1) Understand the characteristics of the object being modeled, exclude minor factors, and select the appropriate model structure, input and output dimensions and order of the model according to the physical or chemical mechanism properties of the target industrial process;

[0010] Step 2) within a sampling period, the input data and detection data of the dynamic system are collected to form a set of data samples, and the recursive expectation step and the recursive maximization step are executed without additional preprocessing of outliers;

[0011] Step 3) Execute the first step of recursive expectation. The prior of the scale factor obeys the Gamma distribution. According to the properties of the conjugate distribution, the posterior of the scale factor also obeys another Gamma distribution. According to the properties of the Gamma distribution expectation, calculate the posterior expectation of the scale factor of the current sampled data.

[0012] Step 4) Execute the second step of recursive expectation, and calculate the posterior expectation of the logarithm of the current sampling data scale factor according to the properties of the conjugate distribution;

[0013] Step 5) Execute the first step of recursive maximization, recursively calculate the sufficient statistics corresponding to the parameter vector of the model, and then update the parameter vector online;

[0014] Step 6) Execute the second step of recursive maximization, recursively calculate the sufficient statistics corresponding to the normal data noise variance, and then update the normal data noise variance online;

[0015] Step 7) Execute the third step of recursive maximization, recursively calculate the intermediate statistics corresponding to the degree of freedom parameter of the Student's t distribution, solve the optimization equation, and update the degree of freedom parameter online.

[0016] Furthermore, in the above-mentioned robust recursive expectation maximization method for dynamic system parameter identification, in step 1), an autoregressive model with an additional input is selected to simulate the real system dynamics, and a linear or weakly nonlinear system is approximated without considering colored noise; the model form is as follows:

[0017]

[0018] In formula (1), for k Check data at all times. To detect the dataset; is noisy data with outliers; is the regression information vector, which contains both the input and output of the past time. is the input data set, n a ,n b are the output polynomial order and the input polynomial order respectively; is the regression parameter vector; k is the time index, N is the current sample number;

[0019] Step 2) Collect input data u at time k at the same time k and the noisy detection data y k , forming a set of observable data samples {u k ,y k}, at this time yk is normal data that obeys Gaussian distribution, or outlier data; there is no need to judge whether it is an outlier online, but directly substitute it into the next step;

[0020] Step 3) and step 4) are to obtain the posterior distribution of the missing variables in the expectation step. The key lies in the transformation of the Q function, which is defined as the posterior expectation of the missing variables of the likelihood function of the complete data set:

[0021]

[0022] In formula (2), Θ is the estimated value of the current iteration parameter set, Θ' is the estimated value of the previous iteration parameter, C obs = {Y,U} is the set of observable variables, C mis is a set of missing variables, also known as a set of latent variables; the expected step is to solve the Q function; for the online expectation maximization method, it does not rely on batch data sets for parameter update, but is performed online; that is, the parameters are updated once at each sampling moment to respond to unknown changes in the system dynamics at a faster speed; after observing the nature of the Q function, a recursive Q function is proposed to transform a fixed sample data set into a flowing sample pool; when new data arrives, it is added to the sample pool without preprocessing; the weight of old data is gradually reduced according to the step size design, and the information it contains continues to fade with the increase of time; based on the batch Q function, the recursive Q function is transformed into the following form

[0023]

[0024] In formula (3), γ N is the step size factor, is the recursive Q function value of N at the current moment, is the recursive Q function value at the previous moment, that is, the N-1 moment, is the posterior expectation of the likelihood function of the complete data set of a single sample at time N, which contains the joint likelihood function of the observed variables and the missing variables;

[0025] Write the recursive Q function in detail and we get

[0026]

[0027] The model parameter vector, normal noise variance and Student's t distribution degree of freedom function are updated online recursively.

[0028] Furthermore, in the above-mentioned robust recursive expectation maximization method for dynamic system parameter identification, in step 3), in order to process outliers in the detection data, the Student's t distribution is used to describe the statistical characteristics of the detection data; the Student's t distribution is decomposed into a mixture of infinite Gaussian distributions, and the key to the realization of this mixture is the proportional factor Given the scaling factor r k Under the condition of , the detection noise data with outliers at time k is decomposed into the following integral form:

[0029] p(e k |μ,Σ,υ k )=∫p(e k |μ,Σ,r k )p(r k |υ k )dr k , (5)

[0030] In formula (5), the detection noise satisfies the conditional Gaussian distribution μ is the noise mean, Σ is the noise variance when outliers are not considered;

[0031] At the same time, the scaling factor r k The priori follows the Gamma distribution

[0032] υ k represents the degree of freedom parameter of the Student's t distribution, and when υ k As r approaches ∞, k tends to 1, and the corresponding t distribution also tends to Gaussian distribution; while r k is a type of missing data, that is, r k ∈C mis , when other latent variables are not considered,

[0033] Furthermore, in the above-mentioned robust recursive expectation maximization method for dynamic system parameter identification, the proportional factor r in step 3) is k The posterior expectation of the related terms is calculated. It is known that the prior of the proportional factor follows the Gamma distribution.

[0034] According to the properties of the conjugate distribution, the posterior of the scale factor follows the following Gamma distribution

[0035]

[0036] The superscript old indicates the estimated value of the corresponding variable at the previous moment. Represents y k and The square of the Mahalanobis distance between

[0037]

[0038] In step 3), the scaling factor r k The posterior expectation of itself is

[0039]

[0040] And defined as a variable The estimated value of the parameter θ at time k-1 k-1 and the detection value y at time k k Under these conditions, the obtained r k The posterior expectation of

[0041] Step 4), the logarithm of the scaling factor logr k The posterior expectation is

[0042]

[0043] Where ψ(·) is the digamma function, and ψ(·)=Γ'(·) / Γ(·), Γ(·) represents the Gamma function, and the calculation formula is

[0044] Step 5), execute the first step of recursive maximization; the recursive update of parameters in the recursive M step is realized by recursive sufficient statistics; considering the model parameter vector θ, directly calculate the partial derivative of the recursive Q function in equation (4) with respect to the parameter vector θ

[0045]

[0046] Adding subscript N to the parameter vector, we get the update equation for the parameter vector:

[0047]

[0048] In formula (10), the recursive update formula of the denominator sufficient statistics is

[0049]

[0050] In formula (10), the recursive update formula of the numerator sufficient statistics is

[0051]

[0052] Step 6), perform the second step of recursive maximization, regarding the variance of normal noise, and also calculate the partial derivative of the recursive Q function in equation (4) with respect to the variance Σ

[0053]

[0054] The update equation of the variance is obtained as

[0055]

[0056] In formula (14), the denominator sufficient statistics The recursive update formula for the sufficient statistics of the corresponding numerator is

[0057]

[0058] Step 7), about the degree of freedom parameter υ of the Student's t distribution k Update, and find the partial derivative of the recursive Q function with respect to υ in equation (4)

[0059]

[0060] We get the equation for the degrees of freedom parameter

[0061]

[0062] Among them, the recursive intermediate statistics

[0063]

[0064] Solution

[0065] And assign the subscript N to get the estimated value of the degree of freedom parameter at the current moment

[0066] Furthermore, in the above-mentioned robust recursive expectation maximization method for dynamic system parameter identification, the step size factor γ in steps 4), 5), and 6) is N The value of is to achieve, basic convergence, so that 1-γ N ∈(0,1), that is, γ N ∈(0,1); In order to achieve fast convergence of dynamic model parameter vector, noise variance and degree of freedom parameter, the time-varying step size factor is selected and integrated into the online update of parameters; the selection of step size factor is designed, γ N There are two optional update methods, one is γ N =(N+1) -β ,β∈(0.5,1), secondly,γ N =β / (k+β+1),β∈(10,20); for different objectives, try different ways of updating the step length factor to balance the convergence speed in the initial stage and the middle stage.

[0067] Compared with the prior art, the present invention has significant advantages and beneficial effects, which are specifically embodied in the following aspects:

[0068] ① The robust recursive expectation maximization method for dynamic system parameter identification of the present invention is concise, easy to apply, and has low computational complexity; it has strong robustness for outlier data in the detection data; through the derivation of the recursive Q function, an online method with recursive sufficient statistics as the core can be obtained to perform online estimation of parameter vectors, normal noise variance, and Student's t distribution degree of freedom parameters; at the same time, the convergence speed of the method can be adjusted by selecting a time-varying step size parameter, so that it is suitable for different identification objects;

[0069] ② It is used for online identification of parameters of industrial process dynamic system models. It uses Student's t distribution as the basic statistical characteristic description method of sampling data and realizes robust online estimation of parameters based on the recursive Q function. It can identify the parameters of dynamic system models online based on the detection data with outliers, and realize online estimation of the variance of normal noise carried by the detection data and the degree of freedom parameter of Student's t distribution, so as to obtain an accurate mathematical model of the target process.

[0070] ③ Based on the robust characteristics of heavy-tailed distribution, the contaminated detection data with outliers is used to realize the online identification of dynamic system parameters and obtain an accurate system mathematical model; the expectation step and the maximization step are executed sequentially in a recursive operation, and the two steps work together with a clear structure, which can estimate the parameters of the dynamic system, the variance of Gaussian noise, and the degree of freedom parameters of the Student t distribution online. It is very good at dealing with the influence of outlier data on dynamic system modeling and model parameter identification, and is suitable for parameter identification of industrial process dynamic models;

[0071] ④ It can realize the online parameter identification and estimation of the industrial process dynamic model when the detection data contains outliers. There is no need to remove or interpolate the outliers in the data before modeling. It simplifies the data preprocessing procedure and improves the accuracy of parameter identification and process modeling. For most industrial process dynamic systems, it can tolerate less than 5% outlier contamination and is insensitive to the specific value of the outliers.

[0072] ⑤ The degree of freedom parameter of the Student's t distribution can be estimated online autonomously by recursively deducing intermediate statistics, which solves the problem that the degree of freedom of the Student's t distribution needs to be given manually in online parameter identification, and realizes the online estimation of the normal noise variance outside the outliers; the convergence speed and accuracy of the method can be adjusted by selecting the time-varying step size factor, making it suitable for different industrial process dynamic systems.

[0073] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or be understood by practicing the specific embodiments of the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures particularly pointed out in the written description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0075] Figure 1: Flow chart of the robust recursive expectation maximization method in the present invention;

[0076] Figure 2 : The probability density function diagram of the heavy-tailed Student's t distribution is adopted in the present invention;

[0077] Figure 3 : The weight scaling of outlier data under Student's t distribution in the present invention;

[0078] Figure 4 : Convergence curve diagram of identification model parameters involved in the present invention;

[0079] Figure 5 : The method involved in the present invention is effective for a certain system;

[0080] Figure 6 : The self-verification output tracking curve of the method involved in the present invention;

[0081] Figure 7 : Cross-validation output tracking curve of the method involved in the present invention. DETAILED DESCRIPTION

[0082] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. The components of the embodiments of the present invention generally described and shown in the drawings here can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed invention, but merely represents the selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making creative work belong to the scope of protection of the present invention.

[0083] It should be noted that similar reference numerals and letters represent similar items in the following drawings, so once an item is defined in one drawing, it does not need to be further defined and explained in the subsequent drawings. At the same time, in the description of the present invention, directional terms and order terms are only used to distinguish the description and cannot be understood as indicating or implying relative importance.

[0084] like Figure 1 As shown, the robust recursive expectation maximization method for dynamic system parameter identification includes the following steps:

[0085] Step 1) Understand the characteristics of the object being modeled, exclude minor factors, select the appropriate model structure, model input and output dimensions and order according to the physical or chemical mechanism properties of the target industrial process;

[0086] Step 2) within a sampling period, the input data and detection data of the dynamic system are collected to form a set of data samples, and the recursive expectation step and the recursive maximization step are executed without additional preprocessing of outliers;

[0087] Step 3) Execute the first step of recursive expectation. The prior of the scale factor obeys the Gamma distribution. According to the properties of the conjugate distribution, the posterior of the scale factor also obeys another Gamma distribution. According to the properties of the Gamma distribution expectation, calculate the posterior expectation of the scale factor of the current sampled data.

[0088] Step 4) Execute the second step of recursive expectation, and calculate the posterior expectation of the logarithm of the current sampling data scale factor according to the properties of the conjugate distribution;

[0089] Step 5) Execute the first step of recursive maximization, recursively calculate the sufficient statistics corresponding to the parameter vector of the model, and then update the parameter vector online;

[0090] Step 6) Execute the second step of recursive maximization, recursively calculate the sufficient statistics corresponding to the normal data noise variance, and then update the normal data noise variance online;

[0091] Step 7) Execute the third step of the recursive maximization step, recursively calculate the intermediate statistics corresponding to the degree of freedom parameter of the Student's t distribution, solve the optimization equation, and update the degree of freedom parameter online.

[0092] Specifically, in step 1), an autoregressive model with external input is selected to simulate the real system dynamics, and a linear or slightly nonlinear system is approximated without considering colored noise; the model form is as follows:

[0093]

[0094] In formula (1), is the detection data at time k, To detect the dataset; is noisy data with outliers; is the regression information vector, which contains both the input and output of the past time. is the input data set, n a ,n b are the output polynomial order and the input polynomial order respectively; is the regression parameter vector; k is the time index, N is the current sample number;

[0095] Step 2) Collect input data u at time k simultaneously k and the noisy detection data y k , forming a set of observable data samples {uk ,y k}, at this time yk is normal data that obeys Gaussian distribution, or outlier data; there is no need to judge whether it is an outlier online, but directly substitute it into the next step;

[0096] Step 3) and step 4) are to obtain the posterior distribution of the missing variables in the expectation step. The key lies in the transformation of the Q function, which is defined as the posterior expectation of the missing variables of the likelihood function of the complete data set:

[0097]

[0098] In formula (2), Θ is the estimated value of the current iteration parameter set, Θ' is the estimated value of the previous iteration parameter, C obs = {Y,U} is the set of observable variables, C mis is a set of missing variables, also known as a set of latent variables; the traditional expectation maximization method uses the Q function as the objective function, and on the basis of batch data, after iterative execution of the expectation and maximization steps, the objective function converges; the expectation step is to solve the Q function; for the online expectation maximization method, it does not rely on batch data sets for parameter updates, but is performed online; that is, the parameters are updated once at each sampling moment to respond to unknown changes in the system dynamics at a faster speed; after observing the nature of the Q function, a recursive Q function is proposed to transform a fixed sample data set into a flowing sample pool; when new data arrives, it is added to the sample pool without preprocessing; the weight of old data is gradually reduced according to the step size design, and the information it contains fades with the increase of time; based on the batch Q function, the recursive Q function is transformed into the following form

[0099]

[0100] In formula (3), γ N is the step size factor, is the recursive Q function value of N at the current moment, is the recursive Q function value at the previous moment, that is, the N-1 moment, is the posterior expectation of the likelihood function of the complete data set of a single sample at time N, which contains the joint likelihood function of the observed variables and the missing variables;

[0101] Write the recursive Q function in detail and we get

[0102]

[0103] The model parameter vector, normal noise variance and Student's t distribution degree of freedom function are updated online recursively.

[0104] Among them, step 3), in order to deal with outliers in the detection data, the Student's t distribution is used to describe the statistical characteristics of the detection data; the Student's t distribution is decomposed into a mixture of infinite Gaussian distributions, and the key to the realization of the mixture is the proportional factor Given the scaling factor r k Under the condition of , the detection noise data with outliers at time k is decomposed into the following integral form:

[0105] p(e k |μ,Σ,υ k )=∫p(e k |μ,Σ,r k )p(r k |υ k )dr k , (5)

[0106] In formula (5), the detection noise satisfies the conditional Gaussian distribution

[0107] μ is the noise mean, Σ is the noise variance when outliers are not considered;

[0108] At the same time, the scaling factor r k The priori follows the Gamma distribution, expressed as

[0109] υ k represents the degree of freedom parameter of the Student's t distribution, and when υ k As r approaches ∞, k tends to 1, and the corresponding t distribution also tends to Gaussian distribution; while r k is a type of missing data, that is, r k ∈C mis , when other latent variables are not considered,

[0110] Furthermore, the proportional factor r in step 3) k The posterior expectation of the related terms is calculated. It is known that the prior of the proportional factor follows the Gamma distribution.

[0111] According to the properties of the conjugate distribution, the posterior of the scale factor follows the following Gamma distribution

[0112]

[0113] The superscript old indicates the estimated value of the corresponding variable at the previous moment. Represents y k and The square of the Mahalanobis distance between

[0114] In step 3), the scaling factor r k The posterior expectation of itself is

[0115]

[0116] And defined as a variable k-1 time parameter estimate θ k-1 and the detection value y at time k k Under these conditions, the obtained r k The posterior expectation of

[0117] Step 4), the logarithm of the scaling factor logr k The posterior expectation is

[0118]

[0119] Where ψ(·) is the digamma function, and ψ(·)=Γ'(·) / Γ(·), Γ(·) represents the Gamma function, and the calculation formula is

[0120] Step 5), execute the first step of recursive maximization; the recursive update of parameters in the recursive M step is realized by recursive sufficient statistics; considering the model parameter vector θ, directly calculate the partial derivative of the recursive Q function in equation (4) with respect to the parameter vector θ

[0121]

[0122] Adding subscript N to the parameter vector, we get the update equation for the parameter vector:

[0123]

[0124] In formula (10), the recursive update formula of the denominator sufficient statistics is

[0125]

[0126] In formula (10), the recursive update formula of the numerator sufficient statistics is

[0127]

[0128] Step 6), perform the second step of recursive maximization, regarding the variance of normal noise, and also calculate the partial derivative of the recursive Q function in equation (4) with respect to the variance Σ

[0129]

[0130] The update equation of the variance is obtained as

[0131]

[0132] In formula (14), the denominator sufficient statistics The recursive update formula for the sufficient statistics of the corresponding numerator is

[0133]

[0134] Step 7), about the degree of freedom parameter υ of the Student's t distribution k Update, and find the partial derivative of the recursive Q function with respect to υ in equation (4)

[0135]

[0136] We get the equation for the degrees of freedom parameter

[0137]

[0138] Among them, the recursive intermediate statistics

[0139]

[0140] Solution

[0141] And assign the subscript N to get the estimated value of the degree of freedom parameter at the current moment

[0142] Step 4), 5), and 6) in step length factor γ N The value of is to achieve, basic convergence, so that 1-γ N ∈(0,1), that is, γ N ∈(0,1); In order to achieve fast convergence of dynamic model parameter vector, noise variance and degree of freedom parameter, the time-varying step size factor is selected and integrated into the online update of parameters; the selection of step size factor is designed, γ N There are two optional update methods, one is γ N =(N+1) -β ,β∈(0.5,1), secondly,γ N =β / (k+β+1),β∈(10,20); for different objectives, try different ways of updating the step length factor to balance the convergence speed in the initial stage and the middle stage.

[0143] Figure 2 The probability density function of the Student's t distribution is shown in the figure. From the figure, we can observe that the Student's t distribution has a thick tail relative to the Gaussian distribution. The probability density function of the Student's t distribution is

[0144]

[0145] in, represents the mean of process noise, represents the process noise variance, represents the degree of freedom parameter of Student's t distribution, It represents the dimension of the detection data. Depending on the problem, the detection data can be 1-dimensional or multi-dimensional. Γ(·) represents the Gamma function.

[0146] Figure 3 The scaling factor r introduced by the Student's t distribution shown k The scale factor acts on the variance when processing outliers. The outliers are pulled back to within the 3-σ range by the scale factor and no longer have a huge destructive effect on the parameter estimation.

[0147] In the Matlab environment, a simulation case of a continuous stirred reactor is established. The mechanism model of the continuous stirred reactor is as follows:

[0148]

[0149] In formula (20) and (21), process flow rate q(t) = 100 L / min, component A feed concentration C A0 (t) = 1 mol / L, feed temperature T0(t) = 350.0K, inlet coolant temperature T c0 (t) = 350.0K, specific heat C p ,C pc =1cal / (gK); heat transfer term hA = 7×10 5 cal / (minK); reactor volume V = 100L, liquid density ρ,ρ c =1×10 3 g / L, reaction heat -ΔH=-2×10 5 cal / mol, activation energy term (E / R) = 1×10 4 K, reaction rate constant k0 = 7.2 × 10 10 min -1 The input-output dynamic characteristics between the coolant flow rate and the concentration of component A product in the simulated fermentation reactor were modeled and approximated using a first-order single-input single-output model with a constant term, and the sampling period was 1 min.

[0150] Figures 4 to 7 ,The performance of the robust recursive expectation maximization method in the present invention, when the coolant flow rate is switched from 98L / min to 105L / min, the dynamics of the system changes, simulating the unknown changes of the system, and using the robust recursive expectation maximization method for identification. The noise-to-signal ratio of normal noise to signal is about 25dB, while the amplitude of the outlier data is about 10 times the amplitude of the normal detection data fluctuation, and the number of outliers changes from 0% of the total data to 5%. In γ N =(N+1) -β,β∈(0.5,1) step size selection mode, Figure 4 The items in the parameter vector converge with the increase of sampling time. Selecting β = 0.7, it can be seen that in the initial stage of calculation, the parameters converge quickly within about hundreds of sampling cycles. When the system dynamics change, the system parameters also change. The robust online expectation maximization method re-estimates the parameters, and the parameters re-converge after about 3000 sampling points. Reducing β can increase the convergence speed, but the stability of parameter estimation will be partially affected. Increasing β can improve parameter stability, but the convergence speed of the parameter estimation process will decrease.

[0151] The effective use range under different data contamination ratios is taken as an example of a specific continuous stirred reactor system. The proportion of outliers in the data is continuously increased from 0, and each increase is 0.5%, looking for the point where the system self-verification mean square error fails. This method can tolerate about 4% outliers without significantly affecting the parameter estimates of the system.

[0152] Figure 6 and Figure 7 The comparison between the method of the present invention and the recursive least squares method is shown respectively. Under the condition that 2% of the detection data are contaminated by outliers, the two methods estimate the parameters of the system dynamic model respectively, and make N-step predictions of the detection values ​​of the system based on the parameter estimation values. As can be seen from the figure, the method of the present invention well tracks the detection values ​​of the system after excluding outliers for both self-validation and cross-validation N-step predictions. In contrast, the recursive least squares method will inevitably have some deviations.

[0153] The robust recursive expectation maximization method for dynamic system parameter identification of the present invention is concise, easy to apply, and has low computational complexity; it has strong robustness for outlier data in the detection data; through the derivation of the recursive Q function, an online method with the recursive sufficient statistic as the core can be obtained to perform online estimation of the parameter vector, the normal noise variance, and the degree of freedom parameter of the Student's t distribution; at the same time, the convergence speed of the method can be adjusted by selecting a time-varying step size parameter, so that it is suitable for different identification objects.

[0154] It is used for online identification of parameters of industrial process dynamic system models. It uses Student's t distribution as the basic statistical characteristic description method of sampling data and realizes robust online estimation of parameters based on the recursive Q function. It can identify the parameters of dynamic system models online based on the detection data with outliers, and realize online estimation of the variance of normal noise carried by the detection data and the degree of freedom parameter of Student's t distribution, so as to obtain an accurate mathematical model of the target process.

[0155] Based on the robust characteristics of heavy-tailed distribution, the contaminated detection data with outliers is used to realize the online identification of dynamic system parameters and obtain an accurate system mathematical model; the expectation step and the maximization step are executed sequentially in a recursive operation, and the two steps work together with a clear structure, which can estimate the parameters of the dynamic system, the variance of Gaussian noise, and the degree of freedom parameters of the Student t distribution online. It can well deal with the influence of outlier data on dynamic system modeling and model parameter identification, and is suitable for parameter identification of industrial process dynamic models.

[0156] It can realize the online parameter identification and estimation of the industrial process dynamic model when the detection data contains outliers. There is no need to remove or interpolate the outliers in the data before modeling. It simplifies the data preprocessing procedure and improves the accuracy of parameter identification and process modeling. For most industrial process dynamic systems, it can tolerate less than 5% outlier contamination and is not sensitive to the specific value of the outliers.

[0157] The degree of freedom parameter of Student's t distribution can be estimated online autonomously by recursively deducing intermediate statistics, which solves the problem that the degree of freedom of Student's t distribution needs to be manually given in online parameter identification, and realizes the online estimation of the variance of normal noise outside outliers.

[0158] The convergence speed and accuracy of the method can be adjusted by selecting the time-varying step size factor, making it suitable for different industrial process dynamic systems.

[0159] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention. It should be noted that similar numbers and letters represent similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in the subsequent drawings.

[0160] The above is only a specific implementation of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of changes or substitutions within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

[0161] It should be noted that, in this article, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "comprise a ..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

Claims

1. A robust recursive expectation maximization method for dynamic system parameter identification, characterized by: The following steps are involved: Step 1) Understand the characteristics of the object being modeled, and select the appropriate model structure, input and output dimensions and order of the model based on the physical or chemical mechanism properties of the target industrial process and excluding minor factors; Step 2) within a sampling period, the input data and detection data of the dynamic system are collected to form a set of data samples, and the recursive expectation step and the recursive maximization step are executed without additional preprocessing of outliers; Step 3) Execute the first step of recursive expectation. The prior of the scale factor obeys the Gamma distribution. According to the properties of the conjugate distribution, the posterior of the scale factor also obeys another Gamma distribution. According to the properties of the Gamma distribution expectation, calculate the posterior expectation of the scale factor of the current sampling data. Step 4) Execute the second step of recursive expectation, and calculate the posterior expectation of the logarithm of the current sampling data scale factor according to the properties of the conjugate distribution; Step 5) Execute the first step of recursive maximization, recursively calculate the sufficient statistics corresponding to the parameter vector of the model, and then update the parameter vector online; Step 6) Execute the second step of recursive maximization, recursively calculate the sufficient statistics corresponding to the normal data noise variance, and then update the normal data noise variance online; Step 7) Execute the third step of recursive maximization, recursively calculate the intermediate statistics corresponding to the degree of freedom parameter of the Student's t distribution, solve the optimization equation, and update the degree of freedom parameter online.

2. The robust recursive expectation maximization method for dynamic system parameter identification according to claim 1, characterized in that: Step 1), select an autoregressive model with external input to simulate the real system dynamics, and approximate the linear or weakly nonlinear system without considering colored noise; the model form is as follows: In formula (1), is the detection data at time k, To detect the dataset; is noisy data with outliers; is the regression information vector, which contains both the input and output of the past time. is the input data set, n a ,n b are the output polynomial order and the input polynomial order respectively; is the regression parameter vector; k is the time index, N is the current sample number; Step 2) Collect input data u at time k at the same time k and the noisy detection data y k , forming a set of observable data samples {u k ,y k }, at this time yk is normal data that obeys Gaussian distribution, or outlier data; there is no need to judge whether it is an outlier online, but directly substitute it into the next step; Step 3) and step 4) are to obtain the posterior distribution of the missing variables in the expectation step, and the key lies in the transformation of the Q function. The Q function is defined as the posterior expectation of the likelihood function of the complete data set with respect to the missing variable: In formula (2), Θ is the estimated value of the current iteration parameter set, Θ' is the estimated value of the previous iteration parameter, C obs = {Y,U} is the set of observable variables, C mis is a set of missing variables, also known as a set of latent variables; the expected step is to solve the Q function; for the online expectation maximization method, it does not rely on batch data sets for parameter update, but is performed online; that is, the parameters are updated once at each sampling moment to respond to unknown changes in the system dynamics at a faster speed; after observing the nature of the Q function, a recursive Q function is proposed to transform a fixed sample data set into a flowing sample pool; when new data arrives, it is added to the sample pool without preprocessing; the weight of old data is gradually reduced according to the step size design, and the information it contains fades with the increase of time; based on the batch Q function, the recursive Q function is transformed into the following form In formula (3), γ N is the step size factor, is the recursive Q function value of N at the current moment, is the recursive Q function value at the previous moment, that is, the N-1 moment, is the posterior expectation of the likelihood function of the complete data set of a single sample at time N, which contains the joint likelihood function of the observed variables and the missing variables; Write the recursive Q function in detail and we get The model parameter vector, normal noise variance and Student's t distribution degree of freedom function are updated online recursively.

3. The robust recursive expectation maximization method for dynamic system parameter identification according to claim 2, characterized in that: Step 3), in order to deal with outliers in the test data, the Student's t distribution is used to describe the statistical characteristics of the test data; the Student's t distribution is decomposed into a mixture of infinite Gaussian distributions, and the key to this mixture is the scale factor Given the scaling factor r k Under the condition of , the detection noise data with outliers at time k is decomposed into the following integral form: In formula (5), the detection noise satisfies the conditional Gaussian distribution μ is the noise mean, Σ is the noise variance when outliers are not considered; At the same time, the scaling factor r k The priori follows the Gamma distribution, expressed as υ k represents the degree of freedom parameter of the Student's t distribution, and when υ k As r approaches ∞, k tends to 1, and the corresponding t distribution also tends to Gaussian distribution; while r k is a type of missing data, that is, r k ∈C mis , when other latent variables are not considered, 4. The robust recursive expectation maximization method for dynamic system parameter identification according to claim 3, characterized in that: The scaling factor r in step 3) k The posterior expectation of the related terms is calculated. It is known that the prior of the proportional factor follows the Gamma distribution. According to the properties of the conjugate distribution, the posterior of the scale factor follows the following Gamma distribution The superscript old indicates the estimated value of the corresponding variable at the previous moment. Represents y k and The square of the Mahalanobis distance between In step 3), the scaling factor r k The posterior expectation of itself is And defined as a variable The estimated value of the parameter θ at time k-1 k-1 and k Time detection value y k Under these conditions, the obtained r k The posterior expectation of Step 4), the logarithm of the scaling factor logr k The posterior expectation is Where ψ(·) is the digamma function, and ψ(·)=Γ'(·) / Γ(·), Γ(·) represents the Gamma function, and the calculation formula is Step 5), execute the first step of recursive maximization; the recursive update of parameters in the recursive M step is realized by recursive sufficient statistics; considering the model parameter vector θ, directly calculate the partial derivative of the recursive Q function in equation (4) with respect to the parameter vector θ Adding subscript N to the parameter vector, we get the update equation for the parameter vector: In formula (10), the recursive update formula of the denominator sufficient statistics is In formula (10), the recursive update formula of the numerator sufficient statistics is Step 6), perform the second step of recursive maximization, regarding the variance of normal noise, and also calculate the partial derivative of the recursive Q function in equation (4) with respect to the variance Σ The update equation of the variance is obtained as In formula (14), the denominator sufficient statistics The recursive update formula for the sufficient statistics of the corresponding numerator is Step 7), about the degree of freedom parameter υ of the Student's t distribution k Update, and find the partial derivative of the recursive Q function with respect to υ in equation (4) We get the equation for the degrees of freedom parameter Among them, the recursive intermediate statistics Solution And assign the subscript N to get the estimated value of the degree of freedom parameter at the current moment 5. The robust recursive expectation maximization method for dynamic system parameter identification according to claim 4, characterized in that: Step 4), 5), and 6) in step length factor γ N The value of , in order to achieve basic convergence, makes 1-γ N ∈(0,1), that is, γ N ∈(0,1); In order to achieve fast convergence of dynamic model parameter vector, noise variance and degree of freedom parameter, the time-varying step size factor is selected and integrated into the online update of parameters; the selection of step size factor is designed, γ N There are two optional update methods, one is γ N =(N+1) -β ,β∈(0.5,1), secondly,γ N =β / (k+β+1),β∈(10,20); for different objectives, different updating methods of the step length factor are tried to balance the convergence speed of the initial stage of the method and the online calculation stage after the system dynamics changes.

Images