Linear time-varying parameter modeling method for time-delay industrial systems based on exponential optimal smoothing regularization

Abstract

The application discloses a linear variable parameter modeling method for time-delay industrial systems based on exponential optimal smoothing regularization, and belongs to the technical fields of system identification and industrial automation. The application aims at the problems that the modeling method adopted by the existing industrial systems estimates model parameters and time-delay parameters in steps, accumulates errors and causes low model precision. The application comprises the following steps: a local linear finite impulse response time-delay model is established, and a global model of the time-delay industrial system is obtained; a model identity hidden variable is introduced in a probability framework, a probability density function of global output of the global model is obtained based on local output distribution characteristics; a smoothing matrix is constructed based on an exponential optimal smoothing regularization method, and a prior distribution of local model parameters is obtained; an observation data set and a missing data set are established, an iterative updating formula of global model parameters is obtained based on a generalized expectation maximization algorithm, and global model parameter estimation values are obtained after the algorithm converges, so that the modeling of the global model is realized. The application is used for linear variable parameter modeling of time-delay industrial systems.

Description

Technical Field

[0001] The invention relates to a linear variable parameter modeling method for a time-delay industrial system based on exponential optimal smoothing regularization, and belongs to the technical field of system identification and industrial automation. Background Art

[0002] As industrial systems grow larger and more complex, they often exhibit complex nonlinearities and time-lag characteristics, such as in petrochemical processes, steelmaking, and fermentation. In actual operation, industrial systems typically have several typical operating points, with smooth transitions between adjacent operating points to complete different production tasks. In these situations, the performance of system models generated using traditional nonlinear modeling methods degrades and they lack the compactness and accuracy required.

[0003] At present, the model obtained by the linear variable parameter modeling method is widely used in nonlinear system identification. It constructs a simple global model by weighted combination of local linear models at multiple working points, and then designs subsequent control, optimization, and fault diagnosis.

[0004] In practical applications, the complex nonlinear characteristics of industrial systems lead to poor modeling accuracy. Furthermore, given the common time-delay characteristics of industrial systems, inaccurate time-delay estimation can affect the stability of closed-loop systems and pose safety risks. In existing industrial system models, the step-by-step estimation of model parameters and time-delay parameters prevents global optimization and leads to cumulative errors. Summary of the Invention

[0005] In view of the problem that the modeling method adopted by existing industrial systems estimates model parameters and time-lag parameters in steps, resulting in large cumulative errors and low model accuracy, the present invention provides a linear variable parameter modeling method for time-lag industrial systems based on exponential optimal smoothing regularization.

[0006] The present invention provides a linear variable parameter modeling method for a time-delay industrial system based on exponential optimal smoothing regularization, comprising:

[0007] Based on the nonlinear characteristics of the time-delay industrial system, a local linear finite impulse response time-delay model is established at a given operating point; all local linear finite impulse response time-delay models are weighted combined to obtain a global model of the time-delay industrial system;

[0008] In the probabilistic framework, the model identity latent variables are introduced, and the probability density function of the global output of the global model is obtained based on the local output distribution characteristics of the local linear finite impulse response time-lag model.

[0009] Based on the prior knowledge that the finite impulse response coefficients of time-delayed industrial systems exhibit negative exponential decay, a smoothing matrix is ​​constructed using the exponential optimal smoothing regularization method, and the prior distribution of the local model parameters of the local linear finite impulse response time-delay model is obtained. An observation data set and a missing data set are established, and the global model parameters are iteratively estimated within the framework of the generalized expectation-maximization algorithm. The optimal global model parameters are finally obtained such that the log-likelihood function of the global output of the global model reaches an extreme value, thus achieving global modeling. The iterative estimation method for the global model parameters is as follows:

[0010] A complete dataset is constructed from the observed dataset and the missing dataset, and the log-likelihood function of the complete dataset is calculated. Then, based on the generalized expectation-maximization algorithm framework, the mathematical expectation of the complete dataset with respect to the model identity latent variables and time-lagged variables under the conditions of the observed dataset and the current global model parameter estimation is calculated, and the unknown posterior probabilities in the mathematical expectation are determined.

[0011] Calculate the joint probability distribution of the global model parameters and construct the loss function to be optimized for the global model parameters; find the extreme value of the loss function to be optimized with respect to all global model parameters, obtain the update formula of the global model parameters, and calculate the updated value of the global model parameters until the calculated global model parameters meet the convergence conditions, determine the optimal global model parameters, and calculate the optimal time-lag variables based on the optimal global model parameters.

[0012] According to the linear variable parameter modeling method of the time-delay industrial system based on exponential optimal smoothing regularization of the present invention, the local linear finite impulse response time-delay model is:

[0013] y km =φ k-τ T θ m +e k ,

[0014] Where y km is the local output of the system at the kth moment at the mth given operating point, m = 1, 2, 3, ..., M, M is the number of given operating points; θ m is the local model parameter at the mth given working point, e k is the output measurement noise at the kth moment; φ k-τ is the regression vector at the kth moment, φ k-τ =[u k-τ-1 u k-τ-2 ...u k-τ-n ] T ,u k is the input of the local linear finite impulse response time-delay model at the kth moment, n is the model order, and τ is the time-delay variable;

[0015] The global model of the time-delay industrial system is obtained by weighted combination of all local linear finite impulse response time-delay models:

[0016]

[0017] Where y k is the global output of the global model at the kth moment, α km is the weight of the local output of the mth local linear finite impulse response time-delay model at the kth moment:

[0018]

[0019] Where ω km is an exponential function:

[0020]

[0021] Where h m is the scheduling variable value corresponding to the mth given operating point, H k is the scheduling variable value of the system working mode at the kth moment, o m is the effective width of the mth local linear finite impulse response time-delay model.

[0022] According to the time-delay industrial system linear variable parameter modeling method based on exponential optimal smoothing regularization of the present invention, the method for obtaining the probability density function of the global output of the global model is:

[0023] Assuming the output measurement noise e k is Gaussian white noise with mean 0 and variance σ 2 Gaussian distribution, the output measurement noise e k The probability density function p(e k |σ 2 )for:

[0024] p(e k |σ 2 )=N(e k |0,σ 2 ),

[0025] Where N represents Gaussian distribution;

[0026] Then the local output y km is subject to the mean φ k-τ T θ m , the variance is σ 2 Gaussian distribution, local output y km The probability density function p(y km |θ m ,σ 2 )for:

[0027] p(y km |θ m , σ 2 )=N(y km |φ k-τ T θ m , σ 2 );

[0028] In the probabilistic framework, the global output y is introduced k The corresponding model identity latent variable I k , and I k Take an integer value between 0 and M, and the global output y k The probability density function p(y k |θ 1，2，...M , σ 2 , o 1，2，...，M )for:

[0029]

[0030] Where p(I k =m)=α km , and p(y k |I k =m,θ m , σ 2 )=N(y km |φ k-τ T θ m , σ 2 ).

[0031] According to the time-delay industrial system linear variable parameter modeling method based on exponential optimal smoothing regularization of the present invention, the smoothing matrix D is constructed:

[0032]

[0033] in:

[0034] d1=2c(λ-1)[1,-1,0 (n-2)×1 ],

[0035] d i =cλ i [0 (i-2)×1 ,-1,2,-1,0 (n-i-2)×1 ], i∈(2, n-1),

[0036] d n =2c(λ-1)[0 (n-2)×1 , 1, -1],

[0037] Where c is the first parameter to be learned, and λ is the second parameter to be learned;

[0038] For the mth local linear finite impulse response time-delay model, the corresponding smoothing matrix is ​​D m , the parameter to be learned is c m , the second parameter to be learned is λ m ;

[0039] Assumption D m θ m The mean is The variance is I n×n The normal distribution of , then:

[0040] p(D m θ m |c m , a m ,λ m )=N(D m θ m |A,I n×n ),

[0041] Where a m I is the parameter to be learned three; n×n is the n-order unit matrix;

[0042] in

[0043] Further obtain the prior distribution of local model parameters:

[0044]

[0045] According to the time-delay industrial system linear variable parameter modeling method based on exponential optimal smoothing regularization of the present invention, the observation data set C is established. obs :

[0046] C obs ={Y, U, H},

[0047] Where Y is the global output data, Y={y k} k＝1，...L , where L is the number of time points; U is the model input data: U = {u k} k＝1，...L , H is the scheduling variable data at the sampling time, H={H k} k＝1，...L ;

[0048] Create missing data set C mis :

[0049] C mis ={I,τ},

[0050] In the formula I={I k} k＝1，...L ;

[0051] Based on the observation dataset C obs and missing dataset C mis Estimate the global model parameters Θ and the time-lagged variables τ:

[0052] Θ={θ m , o m , c m , a m ,λ m ,σ 2} m＝1，2，...，M .

[0053] Beneficial effects of the present invention: The method of the present invention is used for the identification problem of time-delay industrial systems, and can effectively avoid the problems of step-by-step estimation of model parameters and time-delay values ​​and over-fitting of model parameters in traditional methods. The method of the present application first establishes a linear variable parameter time-delay model for the time-delay industrial system; then, under the probability framework, a statistical characterization of the output of the linear variable parameter time-delay model is established, and according to the prior knowledge that the finite impulse response coefficient of the system exhibits a negative exponential decay, a smoothing matrix is ​​constructed based on the exponential optimal smoothing regularization method to obtain the prior distribution of the parameters of the local linear finite impulse response time-delay model, and a mathematical description of the modeling problem is established; then, based on the generalized expectation maximization algorithm framework, the mathematical expectation of the model identity latent variables and time-delay of the complete data set is calculated under the given observation data set and the current model parameter estimation, and the joint probability distribution of all model parameters is calculated, and the two are combined to construct the loss function to be optimized; finally, the extreme value of the loss function with respect to all unknown parameters is obtained to obtain the iterative update formula of all parameters. The method of the present invention can effectively improve the parameter estimation accuracy of the time-delay industrial system, which is of great significance to ensuring the control accuracy of the actual time-delay industrial system.

[0054] The method of the present invention describes the characteristics of the time-delay industrial system through a linear variable parameter time-delay model structure, and finally obtains a linear variable parameter time-delay model that accurately describes the dynamic characteristics of the time-delay industrial system. It realizes the joint optimization estimation of model parameters and time lag, effectively improves the accuracy of the identification results, and has important value for the identification theory of time-delay industrial systems and their industrial applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 Schematic diagram of the process of the linear variable parameter modeling method for time-delay industrial system based on exponential optimal smoothing regularization according to the present invention;

[0056] Figure 2 This is a schematic diagram of a continuous stirred reactor device in a specific embodiment;

[0057] Figure 3 is a diagram of training input data, scheduling variable data, and output data in a specific embodiment;

[0058] Figure 4 This is a comparison diagram of the simulated global output obtained by the global model obtained by the method of the present invention and the real system output;

[0059] Figure 5 This is a comparison diagram of the simulated global output obtained by the global model obtained by the method of the present invention in cross-validation and the real system output. DETAILED DESCRIPTION

[0060] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative efforts shall fall within the scope of protection of the present invention.

[0061] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0062] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, but they are not intended to limit the present invention.

[0063] Specific implementation method 1. Combination Figure 1 As shown, the present invention provides a linear variable parameter modeling method for a time-delay industrial system based on exponential optimal smoothing regularization, comprising:

[0064] Based on the nonlinear characteristics of the time-delay industrial system, a local linear finite impulse response time-delay model is established at a given operating point; all local linear finite impulse response time-delay models are weighted combined to obtain a global model of the time-delay industrial system;

[0065] In the probabilistic framework, the model identity latent variables are introduced, and the probability density function of the global output of the global model is obtained based on the local output distribution characteristics of the local linear finite impulse response time-lag model.

[0066] Based on the prior knowledge that the finite impulse response coefficients of time-delayed industrial systems exhibit negative exponential decay, a smoothing matrix is ​​constructed using the exponential optimal smoothing regularization method, and the prior distribution of the local model parameters of the local linear finite impulse response time-delay model is obtained. An observation data set and a missing data set are established, and the global model parameters are iteratively estimated within the framework of the generalized expectation-maximization algorithm. The optimal global model parameters are finally obtained such that the log-likelihood function of the global output of the global model reaches an extreme value, thus achieving global modeling. The iterative estimation method for the global model parameters is as follows:

[0067] A complete dataset is constructed from the observed dataset and the missing dataset, and the log-likelihood function of the complete dataset is calculated. Then, based on the generalized expectation-maximization algorithm framework, the mathematical expectation of the complete dataset with respect to the model identity latent variables and time-lagged variables under the conditions of the observed dataset and the current global model parameter estimation is calculated, and the unknown posterior probabilities in the mathematical expectation are determined.

[0068] Calculate the joint probability distribution of the global model parameters and construct the loss function to be optimized for the global model parameters; find the extreme value of the loss function to be optimized with respect to all global model parameters, obtain the update formula of the global model parameters, and calculate the updated value of the global model parameters until the calculated global model parameters meet the convergence conditions, determine the optimal global model parameters, and calculate the optimal time-lag variables based on the optimal global model parameters.

[0069] Furthermore, for the time-delay industrial system, a linear variable parameter time-delay model is established, specifically:

[0070] According to the nonlinear characteristics of time-delay industrial systems, the system models are different at different operating points. However, at a given operating point, the local dynamic characteristics of the system can be described by a linear finite impulse response time-delay model. The local linear finite impulse response time-delay model is:

[0071] y km =φ k-τ T θ m +e k ,

[0072] Where y km is the local output of the system at the mth given operating point at the kth moment, m = 1, 2, 3, ..., M, M is the number of given operating points, and also the number of local models; θ m is the local model parameter at the mth given working point, e k is the output measurement noise at the kth moment; φ k-τ is the regression vector at the kth moment, φ k-τ =[u k-τ-1 u k-τ-2 ...u k-τ-n ] T ,u k is the input of the local linear finite impulse response time-delay model at the kth moment, n is the model order, and τ is the time-delay variable;

[0073] The global model of the time-delay industrial system is obtained by weighted combination of all local linear finite impulse response time-delay models:

[0074]

[0075] Where y kis the global output of the global model at the kth moment, α km is the weight of the local output of the mth local linear finite impulse response time-delay model at the kth moment:

[0076]

[0077] Where ω km is an exponential function:

[0078]

[0079] Where h m is the scheduling variable value corresponding to the mth given operating point, H k is the scheduling variable value of the system working mode at the kth moment, o m is the effective width of the mth local linear finite impulse response time-delay model. km is an exponential function, with weight α km is a normalized exponential function. From the weight expression, we can see that α km The value range of is [0 1], that is, the system transitions smoothly between different working points.

[0080] Next, in the probabilistic framework, we establish a statistical characterization of the output of the linear variable parameter time-lag model:

[0081] In this embodiment, the method for obtaining the probability density function of the global output of the global model is:

[0082] Assuming the output measurement noise e k is Gaussian white noise with mean 0 and variance σ 2 Gaussian distribution, the output measurement noise e k The probability density function p(e k |σ 2 )for:

[0083] p(e k |σ 2 )=N(e k |0,σ 2 ),

[0084] Where N represents Gaussian distribution;

[0085] From the mathematical expression of the local linear finite impulse response time-delay model, we can know that the local output y km is subject to the mean φ k-τ T θ m , the variance is σ 2 Gaussian distribution, local output y km The probability density function p(y km |θ m,σ 2 )for:

[0086] p(y km |θ m ,σ 2 )=N(y km |φ k-τ T θ m ,σ 2 );

[0087] In the probabilistic framework, the global output y is introduced k The corresponding model identity latent variable I k , I k is the model identity to which the k-th data point belongs, that is, which local model generates the k-th data point, I k Take an integer value between 0 and M, and the global output is y k The probability density function p(y k |θ 1，2，...M ,σ 2 , O 1，2，...，M )for:

[0088]

[0089] Where p(I k =m)=α km , and p(y k |I k =m,θ m ,σ 2 )=N(y km |φ k-τ T θ m ,σ 2 ).

[0090] According to the prior knowledge that the finite impulse response coefficients of the system exhibit negative exponential decay, a smoothing matrix is ​​constructed based on the exponential optimal smoothing regularization method to obtain the prior distribution of the parameters of the local linear finite impulse response time-delay model. The specific process is as follows:

[0091] The square smoothing method is used to construct the smoothing matrix D:

[0092]

[0093] in:

[0094] d1=2c(λ-1)[1,-1,0 (n-2)×1 ],

[0095] d i =cλ i [0 (i-2)×1 ,-1,2,-1,0(n-i-2)×1 ], i∈(2, n-1),

[0096] d n =2c(λ-1)[0 (n-2)×1 , 1, -1],

[0097] Where c is the first parameter to be learned, and λ is the second parameter to be learned;

[0098] Since the decay rate of the finite impulse response coefficients corresponding to each local model is different, different parameters to be learned should be set when constructing the smoothing matrix D. For the mth local linear finite impulse response time-delay model, the corresponding smoothing matrix is ​​D m , the parameter to be learned is c m , the second parameter to be learned is λ m ;

[0099] Assumption D m θ m The mean is The variance is I n×n The normal distribution of , then:

[0100] p(D m θ m |c m , a m ,λ m )=N(D m θ m |A,I n×n ),

[0101] Where a m I is the parameter to be learned three; n×n is the n-order unit matrix;

[0102] in

[0103] Further obtain the prior distribution of local model parameters:

[0104]

[0105] Build a mathematical description of the model problem, specifically:

[0106] Create an observation dataset C based on the collected data obs :

[0107] C obs ={Y, U, H},

[0108] Where Y is the global output data, Y={y k} k＝1，...L , where L is the number of time points; U is the model input data: U={u k}k＝1，...L , H is the scheduling variable data at the sampling time, H={H k} k＝1，...L ;

[0109] Create missing data set C mis :

[0110] C mis ={I,τ},

[0111] In the formula I={I k} k＝1，...L ;

[0112] Based on the observation dataset C obs and missing dataset C mis Estimate the global model parameters Θ and the time-lagged variables τ:

[0113] Θ={θ m , o m , c m , a m ,λ m ,σ 2} m＝1，2，...，M .

[0114] Furthermore, based on the generalized expectation maximization algorithm framework, the complete dataset is calculated with respect to the model identity latent variable I given the observed dataset and the current model parameter estimates. 1，2，...，L And the mathematical expectation of the time-lag variable τ:

[0115] From the observation dataset C obs and missing dataset C mis Construct the complete dataset:

[0116] {C obs , C mis}={Y,U,H,I,τ},

[0117] Calculate the log-likelihood function logp(Y, U, H, I, τ|Θ) for the complete data set:

[0118]

[0119] Where C1 is a constant term independent of the global model parameter Θ, C1 = log(U, H, τ|Θ).

[0120] Furthermore, the mathematical expectation Q(Θ|Θ) of the complete dataset with respect to the model identity latent variables and time-lag variables under the conditions of the observed dataset and the current global model parameter estimation is calculated. s ):

[0121]

[0122] Where E represents the expectation, d max is the maximum value of the time-lag variable τ, and the value range of τ is [0 d max ] Integer value; Θ s Represents the global model parameters calculated at the sth iteration, where s is the number of current iteration steps;

[0123] The unknown posterior probability term p(I k =m,τ=d|C obs ,Θ s ) is calculated as follows:

[0124]

[0125] in:

[0126]

[0127]

[0128] Assume that the time-lag variable τ follows a uniform distribution, so p(τ=d) is a constant;

[0129] Then the posterior probability of the lagged variable p(τ=d|C obs ,Θ s )for:

[0130]

[0131] In this embodiment, the method for calculating the joint probability distribution of global model parameters is:

[0132]

[0133] Where C2=log(o m , c m , a m ,λ m , σ 2 );

[0134] Construct the loss function to be optimized for the global model parameters

[0135]

[0136] Going further, let’s consider optimizing the loss function About {θ m , o m , c m , a m ,λ m , σ 2} m＝1，2，...，M Find the extreme value and get the update formula of the global model parameters:

[0137]

[0138] After each model parameter update, s=s+1 is set and the next iterative calculation is performed until the algorithm meets the convergence conditions.

[0139] In this embodiment, the convergence condition of the global model parameters is that the relative change between the updated value of the current global model parameter and the global model parameter value obtained in the previous iteration is not greater than an arbitrarily small constant ε:

[0140]

[0141] Taking the updated values ​​of the current global model parameters as the optimal global model parameters;

[0142] After the algorithm converges, the optimal time lag variable is τ * It can be taken as the time lag value corresponding to the maximum posterior probability, that is:

[0143]

[0144] Specific embodiment: A continuous stirred reactor is used to verify the effectiveness of the method of the present invention.

[0145] The experimental setup of the system is shown in the figure Figure 2 According to the law of conservation of mass, the law of conservation of energy and the Arrhenius reaction kinetic equation, the following nonlinear differential equations can be used to describe the dynamic process of the irreversible exothermic reaction A→B in a continuous stirred tank reactor:

[0146]

[0147] Where C A (t) is the concentration of reactant A, in mol / L; q(t) is the feed flow rate, in mol / L; V is the reactor volume, C A0 (t) is the initial concentration of reactant A, k0 is the reaction rate constant, E / R is the activation energy system, T(t) is the reaction temperature in the reactor, unit is K; T0(t) is the initial reaction temperature in the reactor, ΔH is the reaction heat, ρ is the reactant liquid density, C p is the specific heat of the reactants, ρ c is the coolant density, C pc is the specific heat of the coolant, q c (t) is the cooling water flow rate, in L / min. A is the heat exchange coefficient, T C0 (t) is the coolant inflow temperature. c The coolant outlet temperature.

[0148] In this embodiment, the feed flow rate q(t) is selected as the input variable, and the concentration of reactant A C A (t) is the output variable, cooling water flow q c (t) is the scheduling variable that determines the system operating conditions and dynamic characteristics.

[0149] The scheduling variable changes periodically, and the function is q c (t) = 2sin(0.1πt) + 100. The input feed flow rate q(t) is designed to be a random binary sequence in the interval [90, 110], with a data length of 1500 sampling points and a sampling period of 1 minute. Since material transportation is a large time-delay link, it is assumed that the system has an input time delay of 10 minutes. The input data, scheduling variable data, and output data of the continuous stirred tank collected for model identification are as follows: Figure 3 shown.

[0150] The algorithm of the present invention is used to identify the linear variable parameter time-delay model of the continuous stirred reactor. The comparison between the model simulation output obtained by identification and the real system output is as follows: Figure 4 As shown. Figure 4 It can be seen that the model simulation output fits the real system output very well, which shows that the identified model has high accuracy. In order to further verify the accuracy of the identification model, cross-validation is required to re-collect the input data, scheduling variable data and output data of 1200 sampling points and compare the degree of fit between the model simulation output and the real system output. The results are as follows: Figure 5 As shown. Figure 5 It can be seen that the model simulation output in the cross-validation is well matched to the real system output. After two verifications, it is shown that the identification model can well describe the real dynamic characteristics of the continuous stirred reactor.

[0151] Although the present invention is described herein with reference to specific embodiments, it should be understood that these embodiments are merely illustrative of the principles and applications of the invention. It should be understood that many modifications may be made to the illustrative embodiments, and that other arrangements may be devised, without departing from the spirit and scope of the invention as defined by the appended claims. It should be understood that the various dependent claims and features described herein may be combined in ways other than those described in the original claims. It should also be understood that features described in conjunction with individual embodiments may be used in conjunction with other described embodiments.

Claims

1. A linear variable parameter modeling method for time-delay industrial systems based on exponential optimal smoothing regularization, characterized by: Based on the nonlinear characteristics of the time-delay industrial system, a local linear finite impulse response time-delay model is established at a given operating point; all local linear finite impulse response time-delay models are weighted combined to obtain a global model of the time-delay industrial system; In the probabilistic framework, the model identity latent variables are introduced, and the probability density function of the global output of the global model is obtained based on the local output distribution characteristics of the local linear finite impulse response time-lag model. Based on the prior knowledge that the finite impulse response coefficients of time-delayed industrial systems exhibit negative exponential decay, a smoothing matrix is ​​constructed using the exponential optimal smoothing regularization method, and the prior distribution of the local model parameters of the local linear finite impulse response time-delay model is obtained. An observation data set and a missing data set are established, and the global model parameters are iteratively estimated within the framework of the generalized expectation-maximization algorithm. The optimal global model parameters are finally obtained such that the log-likelihood function of the global output of the global model reaches an extreme value, thus achieving global modeling. The iterative estimation method for the global model parameters is as follows: A complete dataset is constructed from the observed dataset and the missing dataset, and the log-likelihood function of the complete dataset is calculated. Then, based on the generalized expectation-maximization algorithm framework, the mathematical expectation of the complete dataset with respect to the model identity latent variables and time-lagged variables under the conditions of the observed dataset and the current global model parameter estimation is calculated, and the unknown posterior probabilities in the mathematical expectation are determined. Calculate the joint probability distribution of the global model parameters and construct the loss function to be optimized for the global model parameters; find the extreme value of the loss function to be optimized with respect to all global model parameters, obtain the update formula of the global model parameters, and calculate the updated value of the global model parameters until the calculated global model parameters meet the convergence conditions, determine the optimal global model parameters, and calculate the optimal time-lag variables based on the optimal global model parameters.

2. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 1, characterized in that: The local linear finite impulse response time-delay model is: y km =φ k-τ T i m +e k , Where y km is the local output of the system at the mth given operating point at the kth moment, m = 1, 2, 3, ..., M, M is the number of given operating points; θ m is the local model parameter at the mth given working point, e k is the output measurement noise at the kth moment; φ k-τ is the regression vector at the kth moment, φ k-τ =[u k-τ-1 u k-τ-2 …u k-τ-n ] T ,u k is the input of the local linear finite impulse response time-delay model at the kth moment, n is the model order, and τ is the time-delay variable; The global model of the time-delay industrial system is obtained by weighted combination of all local linear finite impulse response time-delay models: Where y k is the global output of the global model at the kth moment, α km is the weight of the local output of the mth local linear finite impulse response time-delay model at the kth moment: Where ω km is an exponential function: Where h m is the scheduling variable value corresponding to the mth given operating point, H k is the scheduling variable value of the system working mode at the kth moment, o m is the effective width of the mth local linear finite impulse response time-delay model.

3. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 2, characterized in that: The method to obtain the probability density function of the global output of the global model is: Assuming the output measurement noise e k is Gaussian white noise with mean 0 and variance σ 2 Gaussian distribution, the output measurement noise e k The probability density function p(e k |σ 2 )for: feet k |σ 2 )=N(e k |0,σ 2 ), Where N represents Gaussian distribution; Then the local output y km is subject to the mean φ k-τ T θ m , the variance is σ 2 Gaussian distribution, local output y km The probability density function p(y km |θ m , σ 2 )for: p(y km |θ m ,s 2 )=N(y km |f k-τ T i m ,s 2 ); In the probabilistic framework, the global output y is introduced k The corresponding model identity latent variable I k , and I k Take an integer value between 0 and M, and the global output is y k The probability density function p(y k |θ 1，2，...M ,σ 2 , o 1，2，...，M )for: where p(I k = m) = α km , and p(y k |I k = m, θ m , σ 2 ) = N(y km |φ k-τ T θ m , σ 2 ).

4. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 3, characterized in that: Construct the smoothing matrix D: in: d1=2c(λ-1)[1,-1,0 (n-2)×1 ], d i =cλ i [0 (i-2)×1 ,-1.2,-1,0 (n-i-2)×1 ],i∈(2,n-1), d n =2c(λ-1)[0 (n-2)×1 ,1,-1], Where c is the first parameter to be learned, and λ is the second parameter to be learned; For the mth local linear finite impulse response time-delay model, the corresponding smoothing matrix is ​​D m , the parameter to be learned is c m , the second parameter to be learned is λ m ; Assumption D m θ m The mean is The variance is I n×n The normal distribution of , then: p(D m i m |c m ,a m ,l m )=N(D m i m |A, I n×n ), Where a m is the third parameter to be learned; ×n is the n-order unit matrix; in Further obtain the prior distribution of local model parameters:

5. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 4, characterized in that: Create observation dataset C obs : C obs ={Y,U,H}, Where Y is the global output data, Y={y k } k＝1，...L , where L is the number of time points; U is the model input data: U={u k } k＝1，...L , H is the scheduling variable data at the sampling time, H={H k } k＝1，...L ; Create missing data set C mis : C mis ={I,τ}, Where I = {I k } k＝1，...L ; Based on the observation dataset C obs and missing dataset C mis Estimate the global model parameters Θ and the lagged variables τ: Θ={θ m ,the m ,c m ,a m ,l m ,s 2 } m＝1，2，…，M 。 6. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 5, characterized in that: Construct the complete dataset: {C obs ,C mis }={Y,U,H,I,τ}, Calculate the log-likelihood function logp(Y, U, H, i, τ|Θ) for the complete data set: Where C1 is a constant term independent of the global model parameter Θ, C1 = log(U, H, τ|Θ).

7. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 6, characterized in that: Calculate the mathematical expectation Q(Θ|Θ) of the complete dataset with respect to the model identity latent variables and time-lag variables under the conditions of the observed dataset and the current global model parameter estimation. s ): Where E represents the expectation, d max is the maximum value of the time-lag variable τ, and the value range of τ is [0d max ];Θ s represents the global model parameters calculated at the sth iteration; The unknown posterior probability term p(I k =m,τ=d|C obs ,Θ s ) is calculated as follows: in: Assume that the time-lag variable τ follows a uniform distribution, so p(τ=d) is a constant; Then the posterior probability of the lagged variable p(τ=d|C obs ,Θ s )for:

8. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 7, characterized in that: The method for calculating the joint probability distribution of global model parameters is: Where C2=log(o m , c m , a m ,λ m ,σ 2 ); Construct the loss function to be optimized for the global model parameters 9. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 8, characterized in that: Treat the optimization loss function About {θ m , o m , c m , a m ,λ m , σ 2 } m＝1，2，...，M Find the extreme value and get the update formula of the global model parameters:

10. The method for linear variable parameter modeling of a time-delay industrial system based on exponential optimal smoothing regularization according to claim 9, characterized in that: The convergence condition of the global model parameters is that the relative change between the updated value of the current global model parameters and the global model parameter value obtained in the previous iteration is no greater than the set constant ε: Taking the updated values ​​of the current global model parameters as the optimal global model parameters; The optimal lag variable is τ*:

Images