Complex industrial process sparse dynamics modeling method fusing physical constraints

Abstract

The application discloses a kind of complex industrial process sparse dynamics modeling methods of fusion physical constraint, belong to industrial operation optimization control field, including: obtaining training dataset;To the pre-training of continuous state reconstruction model, while obtaining reconstruction result;Based on the obtained reconstruction result and key process variable, construct dynamics function library;Further construct sparse dynamics parameter model, the joint training of model and sparse dynamics parameter model, obtain the initial sparse coefficient matrix corresponding to model;Threshold screening and joint fine-tuning training are carried out to initial sparse coefficient matrix, to obtain optimal continuous state reconstruction model and fixed structure re-optimized sparse coefficient matrix;Finally, the sparse coefficient matrix obtained is further fine-tuned training, obtains the final initial sparse coefficient matrix and corresponding sparse dynamics parameter model.The sparse coefficient matrix obtained by the application is highly consistent with the real coefficient matrix as a whole, and the similarity reaches more than 95%.

Description

Technical Field

[0001] This invention belongs to the field of industrial operation optimization and control, specifically involving a sparse dynamics modeling method for complex industrial processes that integrates physical constraints. Background Technology

[0002] In process industries, equipment manufacturing, energy and chemical engineering, and intelligent manufacturing, industrial processes are typically characterized by a large number of variables, strong system coupling, complex operating states, and significant dynamic evolution. As industrial systems continue to develop towards larger scale, continuous operation, automation, and intelligence, the accuracy and interpretability of system dynamic models are increasingly demanding in terms of state monitoring, fault early warning, operation optimization, and control decisions during production processes. This is especially true in continuous reactions, thermal processes, complex heat and mass transfer processes, and multivariable coupled control objects, where the system state often evolves continuously over time and is influenced by multiple factors such as input disturbances, equipment state changes, environmental fluctuations, and measurement noise. Failure to accurately obtain the dynamic evolution law of the system will not only affect state prediction, soft measurement modeling, and advanced control performance but may also reduce the reliability of fault detection and safety early warning. Existing technologies for modeling dynamic systems in industrial processes can generally be divided into two categories: mechanistic modeling methods and data-driven modeling methods. Mechanistic modeling methods mainly rely on mass conservation, energy conservation, momentum conservation, reaction kinetics, and equipment mechanisms to establish a mathematical model of the system. This type of method has certain advantages in terms of physical interpretability. When the mechanism is known and the parameters can be accurately obtained, it can usually describe the system's operating laws well. However, in real-world complex industrial scenarios, the system often exhibits strong nonlinearity, time-varying characteristics, parameter drift, unmodeled disturbances, and complex coupling relationships, making accurate mechanism modeling difficult. For some industrial objects, there are also problems such as unclear mechanisms, difficult-to-measure parameters, incomplete boundary conditions, or inability to directly observe state variables, which greatly limits the practical application of traditional mechanism models. With the development of industrial field sensors, control systems, and data acquisition technologies, methods based on process data modeling have gradually become a research hotspot. Data-driven modeling methods learn the mapping relationship between input variables, state variables, and key state variables by utilizing historical operating data or online sampling data. This can reduce the dependence on accurate mechanism knowledge to a certain extent and has advantages such as flexible modeling and strong adaptability. However, traditional pure data-driven methods often focus more on fitting accuracy and are insufficient in characterizing the internal dynamic structure, state derivative relationships, and physical consistency of the system. This can easily lead to problems where the static fitting effect is good, but the dynamic evolution ability is insufficient. Especially when measurement noise is high, sampling frequency is limited, and system nonlinearity is strong, models obtained by simply relying on data fitting may lack interpretability and make it difficult to extract the dominant dynamic terms that truly reflect the essential characteristics of the system, thus limiting their further application in predicting the state of complex industrial processes, detecting faults, and discovering mechanisms.

[0003] In recent years, Physical-Informed Neural Networks (PINNs), as a modeling method combining deep learning and physical constraints, can introduce system dynamic constraints during neural network training, enabling the model to maintain a certain degree of physical consistency while fitting observed data. This type of method has shown strong potential in continuous system state reconstruction, partial differential equation solving, and parameter identification. For industrial process data, PINN can achieve smooth fitting of the system state trajectory by constructing continuous mappings and obtain state derivative information using automatic differentiation techniques, thus providing support for subsequent dynamic relationship modeling. However, when relying solely on PINN for modeling, the model usually tends to obtain a black-box continuous approximation result. Although it has strong function fitting capabilities, it is difficult to directly form explicit, compact, and interpretable dynamic equations. Furthermore, under conditions of weak prior information or insufficient mechanistic information, further extracting dynamic structures with practical physical meaning from PINN results remains challenging. To improve the interpretability of dynamic models, the Sparse Identification of Nonlinear Dynamics (SINDy) method is used to screen dominant dynamic terms from a candidate function library and establish a sparse representation relationship between the system's state derivatives and the candidate functions. This type of method eliminates redundant terms through sparse constraints, and can discover the explicit dynamic structure of the system to a certain extent, exhibiting good model simplicity and interpretability. However, the SINDy method is typically sensitive to the estimation quality of state derivatives. When the original observation data contains noise, sampling errors, or the state trajectory is not smooth enough, the derivatives obtained by direct numerical differencing often have large errors, thus affecting the candidate function library screening results and parameter identification accuracy, and even leading to misidentification, missed identification, or structural instability. Furthermore, if the candidate function library construction and sparse coefficient identification process lacks continuous state reconstruction and derivative stabilization processing, it will further weaken the applicability of this method in complex industrial scenarios. On the other hand, for complex dynamic systems, model training solely through local derivative matching may still suffer from the problem of "consistent local derivatives but gradual accumulation of overall time-progression errors." In other words, while some models may closely approximate the real system at the single-point derivative level, state errors amplify continuously during multi-step recursion or continuous-time integration, ultimately leading to a decline in the model's long-term predictive performance and impacting the accuracy of system simulation and dynamic prediction. Therefore, how to further transform the identified dynamic equations into time-progressive models and utilize state integration errors to re-optimize model parameters is a crucial issue for improving the model's numerical consistency and dynamic prediction capabilities.While some existing methods can simulate models using numerical integration, they are usually used as apocalyptic verification tools rather than feeding the integration error further into the parameter optimization process. As a result, their ability to correct the global consistency of dynamic parameters is limited. Summary of the Invention

[0004] In view of this, the present invention provides a complex industrial process dynamics modeling method that integrates continuous state reconstruction, sparse dynamic structure identification and time integral consistency optimization, in order to solve the problems of strong mechanism dependence, unstable derivative estimation, insufficient model interpretability and low dynamic prediction accuracy in the prior art.

[0005] The purpose of this invention is to provide a sparse dynamics modeling method for complex industrial processes that integrates physical constraints, aiming to overcome many shortcomings of existing technologies in the modeling and dynamics identification of complex industrial processes. Specifically, addressing the common problems of existing methods such as insufficient expression of process mechanisms, difficulty in accurately obtaining key system parameters, poor stability of derivative estimation processes due to noise interference, lack of clear physical meaning and interpretability of dynamic structures, low accuracy of multi-step dynamic prediction, and limited generalization ability, this invention achieves high-precision identification of the potential dynamic laws of complex industrial processes by introducing physical constraint mechanisms, sparse structure expression ideas, and dynamic consistency optimization strategies.

[0006] A sparse dynamics modeling method for complex industrial processes that incorporates physical constraints includes the following steps:

[0007] (1) Collect key process variables and key state variables during the normal operation of the target industrial process, and record the corresponding sampling time to obtain the training dataset;

[0008] (2) Use the obtained training dataset to pre-train the continuous state reconstruction model to obtain the state reconstruction results corresponding to the key state variables;

[0009] (3) Based on the obtained state reconstruction results and key process variables, construct a dynamic function library and obtain the time derivative of the state reconstruction results;

[0010] (4) Based on the dynamic function library and the time derivative of the state reconstruction result, construct the sparse dynamic expression between the two, obtain the corresponding sparse coefficient matrix, and jointly train the pre-trained continuous state reconstruction model and the sparse coefficient matrix to obtain the initial sparse coefficient matrix corresponding to the model.

[0011] (5) Threshold screening is performed on the initial sparse coefficient matrix to remove non-dominant terms. Then, joint fine-tuning training is performed on the continuous state reconstruction model after joint training and the non-zero terms in the sparse coefficient matrix after threshold screening to obtain the optimal continuous state reconstruction model and the sparse coefficient matrix of fixed structure re-optimization.

[0012] (6) Further fine-tune and train the obtained sparse coefficient matrix to obtain the final sparse coefficient matrix and its corresponding sparse dynamic parameter model.

[0013] Furthermore, the key process variables are key process variables used to characterize the dynamic evolution law within the target industrial process; for example, the key process variables include one or more process variables characterizing material changes, energy changes, temperature changes, concentration changes, pressure changes, liquid level changes, or flow rate changes.

[0014] Furthermore, the key state variables are system response variables or output variables obtained through a measuring device. For example, the key state variables include one or more of the following: outlet flow rate, outlet concentration, outlet temperature, pressure response, liquid level response, component content, conversion rate, yield, or other measured variables that are related to the system operation results and can characterize the dynamic evolution process of the system.

[0015] The determination of key process variables and key state variables is generally based on the target industrial process. For example, for some continuous reaction systems, key process variables generally include feed concentration, feed temperature, cooling water inlet temperature, cooling water flow rate, etc.; key state variables of the system generally include material concentration in the reactor, material temperature in the reactor, cooling jacket temperature, etc.

[0016] Furthermore, after obtaining the key process variables and key state variables, these variables need to be processed for prediction. For example, in order to eliminate the differences in the scale of different variables, the key process variables and key state variables of the system need to be normalized and transformed into a standardized form for use in the training of the subsequent continuous state reconstruction model.

[0017] Furthermore, the continuous state reconstruction model employs a multilayer perceptron neural network, including an input layer, several hidden layers, and an output layer. Each layer uses a linear mapping, and except for the output layer, nonlinear activation is performed after the linear mapping. The activation function for each layer except the output layer is the hyperbolic tangent function (Tanh). The input layer takes preprocessed key process variables and sampling times (i.e., the corresponding time variables) as input. State fitting is achieved through multilayer linear mapping and nonlinear activation functions. The continuous state reconstruction function is used to characterize the continuous evolution trajectory of the target industrial process state variables over time.

[0018] Furthermore, when pre-training the continuous state reconstruction model, the error between the key state variables predicted by the network output and the actual observed key state variables is used as the pre-training loss function.

[0019] After obtaining the continuous state reconstruction results of the target industrial process, the automatic differentiation method is used to calculate the time derivative (state derivative) of the state reconstruction results (i.e., the key state variables predicted by the continuous state reconstruction model). Based on the reconstruction results of the key process variables and key state variables, a candidate dynamic function library is constructed, and a sparse linear expression relationship between the state derivative and the candidate function library is established, i.e., the sparse dynamic expression, and the sparse coefficient matrix corresponding to the expression is obtained. On this basis, a joint loss function composed of data fitting loss, physical residual loss and sparse regularization term is further constructed for subsequent synchronous optimization of neural network parameters and sparse coefficient matrix.

[0020] Furthermore, the elements in the dynamic function library are generally determined based on the physical constraint equations (or dynamic equations) corresponding to the target industrial process itself, including constant terms, and any combination of first-order, second-order, and multi-order terms composed of the inputs (corresponding to key process variables) and outputs (corresponding to the state reconstruction results of key state variables) of the continuous state reconstruction model.

[0021] Furthermore, during the joint training, the loss function consists of a data fitting loss term, a physical residual loss term, and a sparse regularization term, along with their corresponding weight coefficients. The data fitting loss term represents the error between the neural network output state and the actual observed state; the physical residual loss term represents the error between the derivative of the neural network output state with respect to time and the predicted value of the sparse dynamic parameter model; and the sparse regularization term is the sparse regularization loss function, i.e., the L1 norm corresponding to the sparse coefficient matrix. By minimizing the joint loss function, the sparse coefficient matrix is ​​updated synchronously, enabling the model to gradually identify the initial sparse dynamic structure of the system while fitting the observed data. When optimizing the model in this step, optimization is performed directly based on the continuous state reconstruction model obtained after pre-training, further improving optimization efficiency and quality.

[0022] Furthermore, the sparse coefficient matrix is ​​used to characterize the contribution of each candidate function term to the dynamic equations of different key state variables. Each column in the sparse coefficient matrix corresponds to the dynamic equation of a key state variable, and each element represents the contribution coefficient of the corresponding candidate function term to the dynamic equation of the corresponding key state variable.

[0023] Furthermore, when performing threshold filtering on the sparse coefficient matrix: original sparse coefficients less than a preset threshold are assigned a value of 0; original sparse coefficients greater than or equal to the preset threshold remain unchanged. After threshold filtering is completed, positions that have been set to zero are no longer allowed to participate in the structure search. Instead, the positions of non-zero terms in the coefficient matrix are fixed, and the sparse coefficient matrix of the continuous state reconstruction model and the fixed structure is optimized only for the non-zero terms.

[0024] Furthermore, during the joint fine-tuning training, the loss function consists of a data fitting loss term, a physical residual loss term, and corresponding weight coefficients. In the joint fine-tuning optimization process, an optimization objective composed of the data fitting term and the physical residual term is adopted, and no further sparse regularization term is added. By minimizing and re-optimizing the aforementioned loss function, the network parameters and the retained non-zero dynamic parameters are jointly fine-tuned, thereby obtaining a simpler, more stable, and physically interpretable sparse dynamic equation.

[0025] Furthermore, during the fine-tuning, a fourth-order Runge-Kutta numerical integration method (RK4) is employed for one-step state progression. The one-step prediction error function is used as the basis for parameter optimization to further refine the parameters in the sparse dynamics model. Specifically, the one-step prediction error is the error between the actual next-time decentralized key system state variables and the predicted key system state variables obtained through the RK4 one-step progression. By minimizing the one-step prediction error function, the parameters in the sparse dynamics model are further refined, ensuring that the dynamics model not only satisfies local derivative constraints but also more accurately reflects the system's dynamic behavior at the state evolution level.

[0026] Furthermore, the complex industrial process system is an industrial process employing a continuous stirred tank reactor.

[0027] This invention also provides a method for predicting key state variables and sparse dynamic parameters of a complex industrial process system, comprising: collecting key process variables during the operation of the target industrial process; using the key process variables and the corresponding sampling time as inputs; and using the model obtained by the modeling method described in any of the above technical solutions to predict the key state variables and sparse dynamic parameters.

[0028] Compared with existing technologies, the beneficial effects of this invention are as follows: This invention combines continuous state reconstruction, physical constraints, sparse dynamic structure identification, and numerical integral consistency optimization. It can identify the potential dynamic laws of a system even when the mechanism information of complex industrial processes is incomplete, system parameters are difficult to measure accurately, and variable coupling relationships are complex, thus reducing the dependence of the modeling process on complete mechanism equations. By using a continuous state reconstruction model to smoothly fit the key state variables of the system and using an automatic differentiation method to obtain the state derivatives, this invention avoids the problem of traditional numerical difference methods being susceptible to noise and sampling errors, improving the stability of derivative estimation. By constructing a physically meaningful dynamic function library and introducing sparse constraints, this invention can screen out dominant terms from candidate dynamic terms, eliminate redundant terms, and obtain sparse dynamic expressions with concise structure and clear physical meaning, improving the interpretability of the model. Furthermore, this invention stabilizes the non-zero structure of the sparse coefficient matrix through threshold screening and fixed structure fine-tuning, reducing the impact of residual non-dominant terms on the dynamic identification results. Simultaneously, it employs a fourth-order Runge-Kutta numerical integration method for one-step state advancement and fine-tunes the sparse dynamic parameters using one-step prediction errors. This ensures that the resulting model not only satisfies local derivative constraints but also exhibits better numerical consistency and dynamic prediction capabilities at the state evolution level. The results of the embodiments show that the sparse coefficient matrix identified by this invention maintains a high degree of consistency with the true coefficient matrix, with a similarity exceeding 95%. This effectively improves the structural accuracy, parameter reliability, and physical interpretability of complex industrial process dynamics modeling, providing a reliable model foundation for industrial process state prediction, fault detection, operation optimization, and control decision-making. Attached Figure Description

[0029] Figure 1 This is a flowchart of the method of the present invention;

[0030] Figure 2 Industrial diagrams for CSTR;

[0031] Figure 3 ~ Figure 5 Compare the output of the pre-trained continuous state reconstruction model (MLP multilayer perceptron) with real data;

[0032] Figure 6 To automatically calculate the correlation coefficient between the derivative and the finite difference derivative;

[0033] Figure 7 and Figure 8 The evolution curves of total loss, data loss, physical constraint loss, and sparsity loss in step 3 are shown.

[0034] Figure 9 The result of the sparse coefficient matrix obtained in step 3;

[0035] Figures 10-12 A comparison chart of the full simulation results and real data of the CSTR system;

[0036] Figures 13-15 A scatter plot comparing the single-step predicted values ​​with the actual values ​​of key process variables in CSTR;

[0037] Figure 16 This is the final sparse coefficient matrix result. Detailed Implementation

[0038] Step 1: During the operation of the target industrial process, collect key process variables, key state variables, and corresponding times at each sampling moment. The key process variables are those used to characterize the dynamic evolution of the target industrial process; the key state variables are system response variables obtained through measuring devices. The key process variables include at least one or more process variables characterizing material changes, energy changes, temperature changes, concentration changes, pressure changes, liquid level changes, or flow rate changes; the system key state variables include one or more measured variables such as outlet flow rate, outlet concentration, outlet temperature, pressure response, liquid level response, component content, conversion rate, yield, or other measured variables related to the system operation results and capable of characterizing the dynamic evolution of the system. Let the k-th sampling moment be t. k Let the target industrial process be at discrete sampling time t k The state variable vector at point is:

[0039] ;

[0040] Where, x k Let x represent the vector of key process variables at the k-th sampling time, where k = 1, 2, ..., N, and N is the total number of sampling points or the total number of sampling times. i,k For x k The i-th element in the expression represents the i-th key process variable at time t. k The values ​​of are: i = 1, 2, ..., n; n represents the dimension of the state variables, i.e., the total number of state variables; R represents the set of real numbers.

[0041] Let the corresponding key state variable vector be:

[0042] ;

[0043] Among them, u k This represents the system's key state variable vector at the k-th sampling time; u j , k For u k The j-th element (j=1,2,3,…,m) in the table represents the j-th critical system state variable at time t. kThe value of ; m represents the dimension of the key state variables, that is, the total number of key state variables.

[0044] If a total of N sampling points are collected, the key process variable sample matrix X and the key state variable sample matrix U can be constructed as follows:

[0045] ;

[0046] Key process variable vector x k Let u be the k-th vector in X, corresponding to the key process variable vector at the k-th sampling time. k Let be the k-th vector in U, which corresponds to the key state variable vector at the k-th sampling time.

[0047] The corresponding sampling time series (i.e., the time variable vector) T is represented as:

[0048] ;

[0049] t k Let be the k-th element in T, corresponding to the k-th time variable, which corresponds to the k-th sampling time mentioned above;

[0050] When the sampling period is fixed, we have:

[0051] ;

[0052] in, This represents the sampling time interval.

[0053] To eliminate differences in the units of measurement of different variables, key process variables and key state variables can be normalized. The standardized form of key state variables can be written as:

[0054] ;

[0055] in:

[0056] and Let represent the mean and standard deviation of the i-th key process variable, respectively;

[0057] and Let represent the mean and standard deviation of the j-th key state variable, respectively;

[0058] This represents the standardized form of the i-th key process variable at the k-th sampling time.

[0059] Let j be the standardized form of the j-th key state variable at the k-th sampling time.

[0060] This leads to the normalized key state variable vector at time k. and key process variable vector :

[0061] ;

[0062] ;

[0063] , Corresponding vectors , The i-th and j-th elements in the array.

[0064] To facilitate subsequent neural network fitting, the time variable at each sampling time, and t... k The corresponding key process variable vectors and key state variable vectors are concatenated together:

[0065] ;

[0066] Where D represents the training dataset used for subsequent continuous state reconstruction models.

[0067] Step 2: In Step 1, we completed the data preprocessing. In this second step, we construct a multilayer perceptron neural network (MLP) as a continuous state reconstruction model, and use sampling time t as the basis for the model. k Key process variable vector As network input, the key state variable vector The neural network is pre-trained using the output labels to obtain key state variables of the target industrial process. A continuous mapping function with respect to time.

[0068] Specifically, in step 1 we have already obtained , and Therefore, we can construct the neural network mapping relationship as follows:

[0069]

[0070] ;

[0071] For any sampling time t k ,have:

[0072]

[0073] in, This is the state reconstruction result at time t (i.e., the predicted value of the key state variables) output by the continuous state reconstruction model. For t k The state reconstruction result output by the continuous state reconstruction model at time t. The network parameters are used to reconstruct the continuous state model.

[0074] In this embodiment, the neural network adopts a multilayer perceptron structure with a total of L layers, including an input layer, L-2 hidden layers, and an output layer; state fitting is achieved through multilayer linear mapping and nonlinear activation functions. For the , The layer, internally expressed as:

[0075] ;

[0076] in:

[0077] ; and They represent the first Layer and first Layer output;

[0078] and They represent the first Layer weight matrix and bias vector;

[0079] The activation function is the hyperbolic tangent function Tanh, which is expressed as follows:

[0080] ;

[0081] The output layer is represented as:

[0082] ;

[0083] W (L) and b (L) Let these represent the weight matrix and bias vector of the Lth layer, respectively;

[0084] h (L-1) This represents the output of the (L-1)th layer.

[0085] The state reconstruction result output by the network Compared with actual observed key state variables u k The error between them is used as the pre-training loss function L data (ω), expressed as:

[0086]

[0087] Where ||·|| is the L2 norm;

[0088] The optimal network parameters for the pre-training stage are obtained by minimizing the loss function using an optimization algorithm, and the predicted values ​​are used as the continuous state reconstruction results.

[0089] ;

[0090] The continuous state reconstruction model is used to characterize the continuous evolution trajectory of key state variables of the target industrial process over time.

[0091] Step 3: After obtaining the continuous state reconstruction results of the target industrial process in Step 2, the automatic differentiation method is used to calculate the time derivative of the state reconstruction results (i.e., the state derivative); a candidate dynamic function library is constructed based on key process variables and key state variables, and a sparse linear expression relationship (i.e., sparse dynamic parameter model) is established between the state derivative and the candidate dynamic function library to determine the structure of the sparse coefficient matrix; on this basis, a joint loss function composed of data fitting loss, physical residual loss and sparse regularization term is further constructed for subsequent synchronous optimization of the continuous state reconstruction model parameters and sparse coefficient matrix.

[0092] The continuous state reconstruction result obtained in step 2 is as follows:

[0093]

[0094] in, This represents the system output vector reconstructed by the neural network (i.e., the continuous state reconstruction model in step 2), where m represents the dimension of the key state variables.

[0095] Since the continuous state reconstruction model is represented by a differentiable neural network, its derivative with respect to time can be calculated using an automatic differentiation method, yielding:

[0096]

[0097] in, Let be the time derivative of the state reconstruction result.

[0098] Obtaining neural network parameters After obtaining its time derivative, the next step is to automatically recover the system's underlying dynamic structure from the data. This involves establishing a library of candidate dynamic functions for key processes with physical meaning.

[0099] ;

[0100] in, Represents a candidate dynamics function library. Each row corresponds to a candidate basis function, which is the dynamic equation for a key state variable. It can include constant terms, linear terms, quadratic terms, and state terms. (Note: This candidate basis function can be freely chosen and is not fixed in form. The above formula is only for illustrative purposes. In the actual structure, the constant term may be absent or replaced with other constant terms; the linear term can also be the input.) and predicted value The addition and subtraction operations, etc., need to be determined based on the kinetic equations corresponding to the target industrial process.

[0101] Further establish the sparse dynamics expression:

[0102] ;

[0103] in, This is a sparse coefficient matrix used to characterize the contribution of each candidate function term to the dynamic equations of different key state variables.

[0104] To ensure consistency between the neural network output state and the actual observed state, a data fitting loss function L is defined. data for:

[0105]

[0106] To ensure that the state derivative satisfies the sparse dynamics expression, a physical residual loss function L is defined. phy for:

[0107]

[0108] in This represents the derivative of the state reconstruction result at time k with respect to time. This represents the candidate dynamics function library vector corresponding to time k;

[0109] To ensure that the dynamic equations retain only a small number of dominant terms, a sparse regular loss function L is defined. sparse for:

[0110] ;

[0111] Where ||·|| represents the L1 norm;

[0112] Further construct the joint loss function L joint :

[0113] ;

[0114] in, , and These are the weight coefficients for the data fitting loss, physical residual loss, and sparse regularization term, respectively. By minimizing the joint loss function, the pre-trained continuous state reconstruction model and the sparse coefficient matrix are updated synchronously, enabling the model to gradually identify the initial sparse dynamic structure of the system while fitting the observed data.

[0115] Step 4: After completing the joint training in Step 3, we obtain a further optimized continuous state reconstruction model and an initial sparse coefficient matrix:

[0116] .

[0117] Wherein, the sparse coefficient matrix Each column corresponds to a dynamic equation for a state variable, and each element This represents the contribution coefficient of the j-th candidate function term to the dynamic equation of the i-th state variable. Since the sparse regularization term introduced during joint training is a soft constraint, the coefficient matrix after joint training may still contain several non-dominant terms with small but not entirely zero values. These non-dominant terms may originate from noise interference, correlations between candidate functions, and local optima during training. Directly retaining these terms would hinder the concise expression of the final dynamic equation and weaken the clarity and interpretability of the dynamic laws.

[0118] Therefore, in step 4, the sparse coefficient matrix is... Threshold filtering is performed using the following thresholding rules:

[0119]

[0120] in: This represents the original sparse coefficients obtained after joint training; Indicates the sparsity coefficient after threshold filtering; This indicates a preset threshold.

[0121] After threshold screening, positions that have been set to zero are no longer allowed to participate in the structure search again. Instead, the positions of non-zero terms in the coefficient matrix are fixed, which can be represented by the mask matrix M.

[0122]

[0123] The subsequent coefficient matrix to be optimized is written as:

[0124]

[0125] in: Here are the structure mask matrix, structure and sparse coefficient matrix. correspond; This represents the Hadamard element-wise product; This is the coefficient matrix to be optimized under a fixed sparse structure.

[0126] At this point, the optimization objective is reverted to consist of a data fitting term and a physical residual term, without adding a sparse regularization term. That is, the updated loss function L... refine The corresponding formula is:

[0127]

[0128] in: The data fitting loss term, used to constrain the consistency between the state reconstruction results output by the neural network and the key state variables observed in reality, can be written as:

[0129]

[0130] in: Let be the true state vector at the k-th sampling time. The normalized value corresponding to the true state vector; is the state reconstruction result output by the neural network at the k-th sampling time; N is the total number of samples.

[0131] The consistency between the state derivative obtained by constrained automatic differentiation using the physical residual loss term and the right-hand side of the dynamic equation under the fixed sparse structure can be written as:

[0132]

[0133] The state derivative obtained by automatic differentiation at the k-th sampling time; This represents the candidate function library vector at the k-th sampling time. The coefficient matrix to be optimized after fixing the sparse structure. and These are the corresponding weighting coefficients.

[0134] By minimizing the re-optimization loss function:

[0135]

[0136] For network parameters By jointly fine-tuning the retained non-zero dynamic parameters (i.e., non-zero terms in the sparse coefficient matrix), a more concise, stable, and physically interpretable sparse dynamic equation (or sparse dynamic model) can be obtained.

[0137] Step 5: After completing the sparse coefficient matrix threshold screening and fixed structure re-optimization in Step 4, the sparse dynamic parameter structure of the system has been basically determined. At this point, in order to further improve the numerical consistency and state evolution accuracy of the model in the dynamic prediction process, the identified sparse dynamic equations are constructed into ordinary differential equation models, and a fourth-order Runge–Kutta (RK4) numerical integration method is used for one-step state advancement. The one-step prediction error is used as the basis for parameter optimization to further fine-tune the dynamic parameters.

[0138] Suppose the dynamic equation obtained after sparse identification is:

[0139]

[0140] in, These are the key process variables after the system is decentralized. These are the key state variables after the system is decentralized. Represented by the sparse coefficient matrix The determined dynamic mapping relationship.

[0141] For a given time t k System real key process variables The RK4 method is used to predict the state at the next time step, and its four intermediate slopes are defined as follows:

[0142]

[0143] in, The sampling time interval is given. Based on the four intermediate slopes mentioned above, the key state variables for one-step integral prediction are obtained as follows:

[0144]

[0145] Based on the fourth-order Runge–Kutta propagation results, an error function is constructed between the one-step integral prediction of key state variables and the actual key states:

[0146]

[0147] in, The key state variables after decentralization in the next real time step. These are the key state variables for prediction obtained through a one-step RK4 process. By minimizing the... The parameters (i.e., the sparse coefficient matrix) in the sparse dynamic model are further modified so that the dynamic model not only satisfies the local derivative constraint, but also more accurately reflects the dynamic behavior of the system at the state evolution level.

[0148] Specific application examples:

[0149] This embodiment uses the classic dataset in the field of fault diagnosis—CSTR (Continuously Stirred Reactor) process—for data experiment verification. (See below...) Figure 2 As shown, the feed stream into the CSTR reactor is a mixture of solvent vapor and reactant vapor, and a product stream is output from the reactor. Under PID feedback control, the temperature and concentration of the output product are controlled using the reactant flow rate and cooling water flow rate. The figure shows the measurement location and control strategy: by manipulating the coolant flow rate Q... c To maintain the reaction temperature T. For greater practicality, the controller is configured to control the coolant flow rate Q. c Below and higher The above describes saturation (saturation refers to the physical limitation of the controller output). Saturation is crucial for simulating situations where a fault gradually develops into a severe one that the control system struggles to handle. The physical constraint dynamics equations for the CSTR process are as follows (t is time):

[0150]

[0151] Specific parameters in CSTR can be found in Table 1. In the model, parameters a0 and b0 are both 1.0 under normal operating conditions. By gradually reducing their values ​​to zero, faults such as catalyst failure and heat transfer fouling can be simulated, respectively. For details on the coefficients of these physical equations, please refer to Table 2; for specific parameter values, please refer to Table 3. The CSTR process collects four key process variables x=[C i ,T i ,T ci Q c ] and three key state variables u=[C,T,T c ].

[0152] Table 1 Detailed Explanation of CSTR Parameters

[0153]

[0154] Table 2

[0155]

[0156] Table 3

[0157]

[0158] Figures 3-5 This describes the reconstruction results of key state variables obtained from the continuous state reconstruction model after pre-training using the method in step 2. Figures 3-5 It can be seen that after pre-training in step 2, the model performs well under varying concentrations C, reaction temperatures T, and jacket temperatures T0. cThe predicted results for key state variables show high consistency with the denoised real data, and the overall trend and local fluctuations of each variable can be well reconstructed. This indicates that the neural network can accurately approximate the system state trajectory, thus providing a foundation for subsequent automatic differentiation and physical constraint modeling.

[0159] Although the state fitting results are good, further comparison of the relationship between the derivative obtained by automatic differentiation and the derivative of finite difference is needed. Figure 6 It can be observed that the correlation between the two is low. Figure 6 As shown, , as well as The relatively small correlation coefficients indicate that the derivative information output by the network cannot be directly equivalent to the true derivative of the system. This phenomenon suggests that high accuracy in state-value fitting does not necessarily equate to high accuracy in derivative estimation. The reason for this is that the goal of the first-stage network optimization is primarily to minimize the fitting error at the state-value level, rather than explicitly constraining derivative consistency. Furthermore, derivative calculations are more sensitive to local errors and residual noise; therefore, even if the predicted curve closely approximates the true curve in terms of function values, its derivative may still exhibit significant deviations.

[0160] After obtaining the smooth state trajectory, the process proceeds to step 3, the joint identification stage of physical constraints and sparse dynamics. In this stage, the model no longer aims solely at fitting state values, but also requires the network output to satisfy the system's underlying dynamic laws. Specifically, based on the CSTR mechanism knowledge and the aforementioned physical constraint dynamic equations, a candidate dynamic function library is constructed. Concentration difference terms, state terms, temperature difference terms, and heat transfer-related terms are included in the candidate set to form the aforementioned dynamic function library, thus creating a sparse dynamic expression with clear physical meaning. The candidate dynamic function library and sparse coefficient matrix used in this paper can be represented as follows:

[0161]

[0162]

[0163] Unlike traditional SINDy libraries that directly construct high-dimensional polynomials, this function library is designed based on system mechanism characteristics. Figure 7 and Figure 8 The results of the changes in each loss function during the joint training process; Figure 9 This is the final sparse coefficient matrix. Figure 7 , Figure 8 and Figure 9It can be seen that after joint training, some key parameters are closer to their true values, indicating that the network has been able to extract the main dynamic information from the smooth state trajectory and initially recover the system equation structure. Meanwhile, the total loss and PDE constraint loss (physical constraint loss) continue to decrease during training. However, since this stage still relies on automatic differentiation as a bridge, and the results of automatic differentiation are greatly affected by state fitting errors, residual noise, and correlations between candidate parameters, the identification of some parameters still exhibits instability.

[0164] Although step 3 applied sparsity constraints to the coefficient matrix through regularization terms, regularization is essentially a soft sparsity mechanism, primarily used to compress smaller coefficients, but it cannot guarantee the strict removal of redundant terms. Therefore, in step 4, a threshold sparsity strategy is further employed to perform hard screening on the identified coefficient matrix, that is, coefficients with absolute values ​​below a set threshold are directly set to zero, and fine-tuning joint training continues on this basis. It should be noted that the network parameters and coefficient matrix are not reinitialized in this stage, but the training results from step 3 are directly inherited. This is because after the previous training stage, the model has obtained a reasonable state representation capability and preliminary dynamic structure information. Reinitialization would not only destroy the existing state trajectory expression, but also lead to the loss of previous identification results. After the sparse dynamic structure of the system is basically determined, step 5, the parameter fine-tuning stage based on RK4, is entered. Unlike the previous stages, this stage no longer directly uses the consistency between the automatic differential derivative and the right-hand side of the dynamic equation for parameter updates, but instead calculates the single-step state prediction results through fourth-order Runge-Kutta numerical integration, and uses the single-step prediction error as the basis for parameter optimization.

[0165] Depend on Figures 10-12 You can see the global R 2 The results are not good because for dynamic systems like CSTRs, which have nonlinear and strongly coupled characteristics, small parameter deviations, initial state errors, and model structure mismatches accumulate over time during the full-process rolling simulation and may be gradually amplified. Therefore, the full-process simulation results depend not only on the accuracy of parameter identification but also on the combined effects of various factors such as denoising errors, unmodeled dynamics, and numerical integration errors. In contrast, one-step R... 2 By using the actual current state as the starting point for prediction at each step, the model's ability to characterize local state transitions is more directly reflected. Figures 13-15 As shown. Therefore, when evaluating the parameter identification results, one should not rely solely on the full-process simulation R... 2 Instead of drawing conclusions, a comprehensive analysis should be conducted, taking into account the accuracy of one-step prediction, the relative error of parameters, and physical rationality.

[0166] from Figure 16It can be seen that the sparse coefficient matrix finally identified by the proposed method has a high degree of consistency with the true coefficient matrix in terms of overall structure. The main non-zero terms in the true coefficient matrix are well preserved in the identification results, and their corresponding positions are basically consistent, indicating that the model can accurately identify the dominant candidate terms in the system dynamic equations and effectively eliminate non-dominant and redundant terms. At the same time, the numerical values ​​of the main non-zero coefficients are close to the true parameters, indicating that the proposed method not only achieves the recovery of the dynamic structure but also has good parameter identification capabilities. Especially considering the low correlation between the automatic differential derivative and the finite difference derivative in the early stages and the instability of some parameter identifications, the numerical accuracy of the coefficient matrix is ​​further improved after subsequent threshold screening, fixed structure fine-tuning, and RK4 single-step prediction error fine-tuning. The final result has a similarity of over 95% with the true coefficient matrix, fully demonstrating that the proposed method can effectively improve the structural accuracy, parameter reliability, and physical interpretability of sparse dynamic models. This result further verifies the effectiveness of combining physical constraints, sparse identification, and integral consistency optimization for dynamic modeling of complex industrial processes.

Claims

1. A sparse dynamics modeling method for complex industrial processes that incorporates physical constraints, characterized in that, Includes the following steps: (1) Collect key process variables and key state variables during the normal operation of the target industrial process, and record the corresponding sampling time to obtain the training dataset; (2) Use the obtained training dataset to pre-train the continuous state reconstruction model to obtain the state reconstruction results corresponding to the key state variables; (3) Based on the obtained state reconstruction results and key process variables, construct a dynamic function library and the time derivative of the state reconstruction results; (4) Based on the dynamic function library and the time derivative of the state reconstruction result, construct the sparse dynamic expression between the two to obtain the corresponding sparse coefficient matrix. Jointly train the pre-trained continuous state reconstruction model and the sparse coefficient matrix to obtain the initial sparse coefficient matrix. (5) Threshold screening is performed on the initial sparse coefficient matrix to remove non-dominant terms. Then, joint fine-tuning training is performed on the continuous state reconstruction model after joint training and the non-zero terms in the sparse coefficient matrix after threshold screening to obtain the optimal continuous state reconstruction model and the sparse coefficient matrix of fixed structure re-optimization. (6) Further fine-tune and train the re-optimized sparse coefficient matrix to obtain the final initial sparse coefficient matrix and the corresponding sparse dynamic parameter model; During the fine-tuning, a fourth-order Runge-Kutta numerical integration method is used for one-step state advancement, with the one-step prediction error as the basis for parameter optimization, to further correct the parameters in the sparse coefficient matrix; the one-step prediction error is the error between the real next time-decentralized system key state variables and the predicted system key state variables obtained by the one-step advancement of the fourth-order Runge-Kutta numerical integration method.

2. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints as described in claim 1, characterized in that, The continuous state reconstruction model employs a multilayer perceptron neural network, comprising an input layer, several hidden layers, and an output layer; except for the output layer, the activation function of each layer is the hyperbolic tangent function.

3. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints as described in claim 1, characterized in that, When pre-training the continuous state reconstruction model, the error between the key state variables predicted by the network and the key state variables observed in reality is used as the pre-training loss function.

4. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints as described in claim 1, characterized in that, The elements in the dynamics function library are determined according to the dynamics equations corresponding to the target industrial process itself, including constant terms, and any combination of first-order terms, second-order terms, and third-order terms composed of the inputs and outputs of the continuous state reconstruction model.

5. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints according to claim 1, characterized in that, During the joint training, the loss function consists of a data fitting loss term, a physical residual loss term, and a sparse regularization term, along with their corresponding weight coefficients.

6. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints according to claim 1, characterized in that, When performing threshold filtering on the sparse coefficient matrix: the original sparse coefficients less than the preset threshold are assigned a value of 0; the original sparse coefficients greater than or equal to the preset threshold remain unchanged.

7. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints according to claim 1, characterized in that, During the joint fine-tuning training, the loss function consists of a data fitting loss term, a physical residual loss term, and corresponding weight coefficients.

8. The sparse dynamics modeling method for complex industrial processes incorporating physical constraints according to claim 1, characterized in that, The industrial process described is an industrial process that uses a continuous stirred tank reactor.

Images