A control method for continuous stirred reactor based on Koopman operator

Abstract

The present invention provides a continuous stirred reactor control method based on Koopman operator, which belongs to the technical field of continuous stirred reactor control, and comprises: based on the dynamic model of a continuous stirred reactor (CSTR) system, a nonlinear system is mapped to a high-dimensional linear space, and the finite-dimensional approximation of the Koopman operator is solved by an extended dynamic mode decomposition method, so as to obtain a high-dimensional linear state space model of the CSTR system. On this basis, a model predictive controller objective function including state tracking error, control input and terminal penalty term is designed, and the objective function is converted into a standard linear quadratic programming problem for solution to obtain an optimal control sequence. The system state is dynamically updated through a rolling optimization strategy, and the state and control weighting matrix are adjusted to optimize the controller performance, and finally a high-precision trajectory tracking predictive control of the CSTR system is realized. This scheme can ensure the stability and control accuracy of the reactor system.

Description

Technical Field

[0001] The invention belongs to the technical field of continuous stirred tank reactor control, and in particular relates to a continuous stirred tank reactor control method based on a Koopman operator. Background Art

[0002] As one of the most common and important reactors in the field of chemical engineering, the continuous stirred tank reactor (CSTR) is widely used in many industries such as petrochemicals, fine chemicals, and bioengineering. CSTR has the characteristics of continuous feeding and discharging, which can achieve full mixing and continuous reaction of materials and plays a key role in large-scale industrial production. However, the control of CSTR system has always been a challenging problem, mainly due to its inherent complexity and nonlinear characteristics.

[0003] CSTR systems usually involve complex physical and chemical reaction processes, including heat transfer, material transfer, and chemical reaction kinetics. There is a strong coupling relationship between these processes, which makes the system exhibit significant nonlinear characteristics. For example, the reaction rate changes exponentially with temperature, and the temperature is affected by multiple factors such as reaction heat and cooling effect. In addition, CSTR systems are also multivariable, time-varying, and uncertain, which further increases the difficulty of control. In actual production, external factors such as raw material composition fluctuations, equipment failures, and environmental changes can also cause disturbances to the system, making the control problem more complicated.

[0004] Although the traditional PID control method is widely used in industrial practice, it is often difficult to achieve ideal control effects when dealing with strong nonlinear systems such as CSTR. Especially when the system parameters change or face large disturbances, the performance of the reactor control will drop significantly, making it difficult to ensure the stability and control accuracy of the reactor. Summary of the invention

[0005] In view of this, the present invention provides a continuous stirred reactor control method based on the Koopman operator, which can solve the technical problems in the prior art that when the system parameters change or face large disturbances, the performance of the reactor control will significantly decrease, making it difficult to ensure the stability and control accuracy of the reactor system.

[0006] The present invention is achieved in that:

[0007] A first aspect of the present invention provides a continuous stirred tank reactor control method based on a Koopman operator, comprising the following steps:

[0008] S10, establish a dynamic model of the continuous stirred reactor and determine the state variables and control variables of the system;

[0009] S20, selecting an appropriate set of basis functions, defining a state lifting function, and mapping the original nonlinear CSTR system to a high-dimensional linear space;

[0010] S30, collect data on the evolution of the CSTR system state over time and construct a data set;

[0011] S40, using the extended dynamic mode decomposition method, solve the finite-dimensional approximation of the Koopman operator and obtain the state space matrix of the high-dimensional linear CSTR system;

[0012] S50, based on the obtained high-dimensional linear model, establish a prediction equation for the future state of the CSTR system;

[0013] S60. Design the objective function of the model predictive controller, including state tracking error, control input and terminal penalty terms;

[0014] S70, transforming the objective function into a standard linear quadratic programming problem and solving the optimal control sequence;

[0015] S80, implement the rolling optimization strategy, update the system status at each sampling time, and repeat steps S50 to S70;

[0016] S90, by adjusting the state weighting matrix and the control weighting matrix, the controller performance is optimized to achieve predictive control of high-precision trajectory tracking of the CSTR system.

[0017] Wherein, the step S10 specifically includes:

[0018] Step 101, determining the state variables of the continuous stirred tank reactor system, including reactant concentration and reactor temperature;

[0019] Step 102, determining the control variable of the continuous stirred tank reactor system, that is, the coolant temperature;

[0020] Step 103, listing the differential equation of reactant concentration, taking into account input concentration, reaction consumption and discharge loss;

[0021] Step 104, listing the differential equation of the reactor temperature, taking into account the input heat, reaction heat, cooling heat and discharge heat;

[0022] Step 105, introducing relevant parameters, including reaction rate constant, activation energy, reaction heat, flow rate, reactor volume, feed concentration, feed temperature, gas constant, density, specific heat capacity, heat transfer coefficient and heat transfer area;

[0023] Step 106, normalizing the differential equation to obtain a dimensionless dynamic model;

[0024] Step 107: Verify the accuracy of the established dynamic model to provide a system description for subsequent Koopman operator analysis.

[0025] Furthermore, the step S20 specifically includes:

[0026] Step 201, selecting a thin plate spline radial basis function as a basis function;

[0027] Step 202, determine the number of basis functions through cross-validation, and select 10 to 20 basis functions;

[0028] Step 203: define a state improvement function to map the original state space to a high-dimensional feature space;

[0029] Step 204: Determine the center point of the basis function by uniform sampling or clustering method;

[0030] Step 205: construct a state improvement function, take the original state variables and the control variables as input, and output a high-dimensional feature vector;

[0031] Step 206: Verify the validity of the constructed state improvement function to ensure that it can accurately capture the nonlinear characteristics of the system;

[0032] Step 207: Convert the nonlinear system into a high-dimensional linear system, laying a foundation for the subsequent application of linear control theory.

[0033] Furthermore, the step S30 specifically includes:

[0034] Step 301, set the sampling time interval, and select a sampling period between 0.1 seconds and 1 second according to the dynamic characteristics of the system;

[0035] Step 302, determine the total sampling time, and select 5 to 10 times of the system response time to fully capture the system dynamics;

[0036] Step 303, running a continuous stirred tank reactor system simulation under different initial conditions and control inputs;

[0037] Step 304, recording the state variables of each time step, including reactant concentration and reactor temperature, and the control variable, i.e., coolant temperature;

[0038] Step 305: preprocessing the collected data, including denoising and normalization;

[0039] Step 306: Arrange the processed data into a data set according to the time series;

[0040] Step 307: Verify the quality and coverage of the collected data to ensure that it can fully characterize the dynamic behavior of the system.

[0041] Furthermore, the step S40 specifically includes:

[0042] Step 401, construct a data matrix, including the state improvement vectors at all times;

[0043] Step 402, calculate matrix operations to obtain the intermediate matrix required for Koopman operator approximation;

[0044] Step 403, solving the linear equations to obtain a finite-dimensional approximate matrix of the Koopman operator;

[0045] Step 404, calculating the control matrix by the least square method;

[0046] Step 405: construct an output matrix so that the original state can be reconstructed from the high-dimensional state;

[0047] Step 406, verifying the accuracy of the obtained high-dimensional linear model by comparing the model prediction value with the actual system response;

[0048] Step 407: Adjust model parameters to optimize the performance of the high-dimensional linear system and lay the foundation for subsequent model predictive control.

[0049] Furthermore, the step S50 specifically includes:

[0050] Step 501, determine the prediction time domain, and select 10 to 30 sampling periods;

[0051] Step 502, determine the control time domain, select 5 to 15 sampling periods, and do not exceed the prediction time domain;

[0052] Step 503: using a high-dimensional linear model, recursively calculate the state prediction value in the future prediction time domain;

[0053] Step 504: construct a prediction equation matrix form, and express the state prediction value in the future prediction time domain as a function of the current state and the future control input;

[0054] Step 505: Consider system constraints, including upper and lower limits of state variables and control variables;

[0055] Step 506: introducing an interference model to improve prediction accuracy;

[0056] Step 507: Verify the accuracy of the prediction equation and evaluate the prediction performance by comparing the predicted value with the actual system response.

[0057] Furthermore, the step S60 specifically includes:

[0058] Step 601: define a state tracking error term to measure the deviation between the predicted state and the reference trajectory;

[0059] Step 602: define control input items for limiting the amplitude and change rate of the control amount;

[0060] Step 603: define a terminal penalty term to ensure the stability of the closed-loop system;

[0061] Step 604: introduce a weight matrix to balance the importance of state tracking error, control input and terminal penalty;

[0062] Step 605: Consider system constraints and convert them into penalty terms or hard constraints of the objective function;

[0063] Step 606: construct a complete objective function expression;

[0064] Step 607: Verify the validity of the objective function to ensure that it can comprehensively reflect the control performance requirements of the system.

[0065] Furthermore, the step S70 specifically includes:

[0066] Step 701: Expand the objective function to obtain a quadratic expression for the control sequence;

[0067] Step 702: construct the Hessian matrix and gradient vector of the quadratic programming problem;

[0068] Step 703: convert the system constraints into the form of inequality constraints;

[0069] Step 704: Select an appropriate quadratic programming solution algorithm, such as an interior point method or an effective set method;

[0070] Step 705: Setting the parameters of the solver, including the maximum number of iterations and the convergence threshold;

[0071] Step 706, solving the quadratic programming problem to obtain the optimal control sequence;

[0072] Step 707: Verify the feasibility and optimality of the obtained control sequence to ensure that it satisfies the system constraints and minimizes the objective function.

[0073] Furthermore, the step S80 specifically includes:

[0074] Step 801, obtaining the actual state of the continuous stirred tank reactor system at the current moment;

[0075] Step 802: Based on the current state, use the prediction equation to predict the system behavior in the future prediction time domain;

[0076] Step 803: construct an optimization problem using the designed objective function;

[0077] Step 804, solving the optimization problem to obtain a control sequence;

[0078] Step 805, only the first element of the control sequence is executed as the current control input;

[0079] Step 806: Wait for the next sampling moment;

[0080] Step 807, repeating steps 801 to 806 to achieve rolling optimization control of the continuous stirred tank reactor system.

[0081] Furthermore, the step S90 specifically includes:

[0082] Step 901, setting the initial state weight matrix and control weight matrix;

[0083] Step 902, running a closed-loop control simulation, recording the continuous stirred tank reactor system response and control input;

[0084] Step 903, calculating performance indicators, including the root mean square value of the tracking error and the energy consumption of the control input;

[0085] Step 904: Adjust the elements of the state weight matrix and the control weight matrix according to the performance index;

[0086] Step 905, repeat the simulation and adjustment process until satisfactory control performance is achieved;

[0087] Step 906, verifying the robustness of the obtained controller under different working conditions;

[0088] Step 907: Determine the final state weighted matrix and control weighted matrix to achieve high-precision trajectory tracking control of the continuous stirred tank reactor system.

[0089] Optionally, in step S10, the original kinetic model is converted into a dimensionless form by introducing a dimensionless parameter combination, thereby simplifying the model structure and improving the calculation efficiency.

[0090] Optionally, in step S20, an orthogonal polynomial basis function is used as an alternative basis function, and the basis function type most suitable for the continuous stirred reactor system is selected by comparing the performances of different basis functions.

[0091] Optionally, in step S30, a Latin hypercube sampling method is used to generate initial conditions and control inputs to ensure uniformity and representativeness of the collected data.

[0092] Optionally, in step S40, a sparse recognition method, such as LASSO regression, is used to identify a finite-dimensional approximation of the Koopman operator to improve the generalization ability of the model.

[0093] Optionally, in step S50, an adaptive prediction time domain strategy is introduced to dynamically adjust the length of the prediction time domain according to changes in system status to balance computational complexity and prediction accuracy.

[0094] Optionally, in step S60, a multi-objective optimization method is used to construct an objective function, while considering multiple indicators such as tracking performance, control cost and stability.

[0095] Optionally, in step S70, a real-time solution strategy is introduced, such as a real-time iterative linear quadratic regulator, to improve the real-time performance of the controller.

[0096] Optionally, in step S80, an event triggering mechanism is used to determine whether to perform optimization calculations according to the degree of change in the system state, thereby reducing unnecessary computing overhead.

[0097] Optionally, in step S90, a heuristic algorithm, such as particle swarm optimization or genetic algorithm, is used to optimize the state weighted matrix and the control weighted matrix to improve the parameter adjustment efficiency.

[0098] Optionally, it also includes online adaptive updating of the Koopman operator model, which collects system operation data in real time and regularly updates model parameters to improve the long-term prediction accuracy and control performance of the model.

[0099] Compared with the prior art, the beneficial effects of the continuous stirred reactor control method based on Koopman operator provided by the present invention are:

[0100] 1) Using the Koopman operator, a finite high-dimensional linear dynamic model of CSTR is established to transform the nonlinear description of the system into a linear representation in high-dimensional space. The high-dimensional linearized model can more accurately capture the nonlinear characteristics of the CSTR system and better describe the entire dynamic behavior of the actual engineering system.

[0101] 2) The model predictive controller based on the Koopman operator does not need to rely on the precise mathematical model of the system, and can directly learn the behavior of the system from the collected system data; the controller can avoid non-convex optimization problems by using a high-dimensional linear model, thereby reducing its computational complexity and improving computational efficiency;

[0102] 3) The model predictive control method based on Koopman operator can control the CSTR system more accurately, making the system converge faster and the static deviation smaller. It is beneficial to improve the stability and control accuracy of the system.

[0103] In summary, the solution of the present invention solves the technical problem in the prior art that when the system parameters change or face a large disturbance, the performance of the reactor control will significantly decrease, making it difficult to ensure the stability and control accuracy of the reactor system. BRIEF DESCRIPTION OF THE DRAWINGS

[0104] Figure 1 A flow chart of the method provided by the present invention;

[0105] Figure 2 It is a schematic diagram of the structure of a continuous stirred reactor system;

[0106] Figure 3 Schematic diagram of the system evolution process based on the Koopman operator;

[0107] Figure 4 It is the control structure block diagram;

[0108] Figure 5 This is a graph showing the control results of the discharge concentration of the continuous stirred reactor;

[0109] Figure 6 This is a graph showing the control results of the discharge temperature of the continuous stirred reactor;

[0110] Figure 7 The input signal diagram of the CSTR system under three control methods;

[0111] Figure 8 This is a graph showing the control results of the discharge concentration of the continuous stirred reactor;

[0112] Fig. 9 This is a graph showing the control results of the discharge temperature of the continuous stirred reactor;

[0113] Fig.10 Input signal diagram of CSTR system under three control methods. DETAILED DESCRIPTION

[0114] In order to make the purpose, technical solution and advantages of the embodiments of the present invention more clear, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0115] like Figure 1 As shown, it is a flow chart of a continuous stirred reactor control method based on Koopman operator provided by the present invention, and the method comprises the following steps:

[0116] S10, establish a dynamic model of the continuous stirred reactor and determine the state variables and control variables of the system;

[0117] S20, selecting an appropriate set of basis functions, defining a state lifting function, and mapping the original nonlinear CSTR system to a high-dimensional linear space;

[0118] S30, collect data on the evolution of the CSTR system state over time and construct a data set;

[0119] S40, using the extended dynamic mode decomposition method, solve the finite-dimensional approximation of the Koopman operator and obtain the state space matrix of the high-dimensional linear CSTR system;

[0120] S50, based on the obtained high-dimensional linear model, establish a prediction equation for the future state of the CSTR system;

[0121] S60. Design the objective function of the model predictive controller, including state tracking error, control input and terminal penalty terms;

[0122] S70, transforming the objective function into a standard linear quadratic programming problem and solving the optimal control sequence;

[0123] S80, implement the rolling optimization strategy, update the system status at each sampling time, and repeat steps S50 to S70;

[0124] S90, by adjusting the state weighting matrix and the control weighting matrix, the controller performance is optimized to achieve predictive control of high-precision trajectory tracking of the CSTR system.

[0125] The specific implementation methods of the above steps are described in detail below:

[0126] The specific implementation of step S10 is as follows: First, according to the law of conservation of matter and energy, a kinetic model of a continuous stirred tank reactor (CSTR) is established. This process includes the following sub-steps: determining the state variables of the system, including the reactant concentration C A and reactor temperature T; determine the control variable of the system, namely the coolant temperature T c ; List the reactant concentration C A The differential equation of the reaction kettle temperature T is listed, taking into account the feed concentration, reaction consumption and discharge loss; the differential equation of the reaction kettle temperature T is listed, taking into account the feed heat, reaction heat, cooling heat and discharge heat; relevant parameters are introduced, such as the reaction rate constant k0, activation energy E, reaction heat ΔH, etc.; the differential equation is normalized to obtain a dimensionless kinetic model. The purpose of this step is to provide an accurate system description for the subsequent Koopman operator analysis. The established kinetic model can be expressed as:

[0127]

[0128]

[0129] Where q is the flow rate, V is the reactor volume, C A0 is the feed concentration, T0 is the feed temperature, R is the gas constant, ρ is the density, C p is the specific heat capacity, U is the heat transfer coefficient, and A is the heat transfer area. This model accurately describes the dynamic behavior of the CSTR system and lays the foundation for subsequent control design.

[0130] The specific implementation method of step S20 is: select an appropriate set of basis functions, define a state lifting function, and map the original nonlinear CSTR system to a high-dimensional linear space. This process includes the following sub-steps: select the basis function type, and this method uses thin plate spline radial basis function (TPS-RBF) as the basis function; determine the number of basis functions, usually 10 to 20 basis functions are selected, and the specific number can be determined by cross-validation; define a state lifting function to map the original state space to a high-dimensional feature space; determine the center point of the basis function, which can be selected by uniform sampling or clustering methods; construct a state lifting function, take the original state variables and control variables as input, and output a high-dimensional feature vector. The purpose of this step is to transform the nonlinear system into a high-dimensional linear system, which is convenient for the subsequent application of linear control theory. The state lifting function can be expressed as:

[0131] ψ(x)=[ψ1(x), ψ2(x),…,ψ N (x)] T

[0132] Among them, ψ i (x)=||xc i || 2 log(||xc i ||) is the i-th TPS-RBF basis function, x is the original state vector, c i is the center point of the ith basis function. Through this mapping, the dynamic characteristics of the original nonlinear system can be linearly represented in a high-dimensional space.

[0133] The specific implementation method of step S30 is: collect data on the evolution of the CSTR system state over time and construct a data set. This process includes the following sub-steps: set the sampling time interval, usually selected between 0.1 seconds and 1 second, depending on the dynamic characteristics of the system; determine the total sampling time, generally select 5 to 10 times the system response time to fully capture the system dynamics; run the CSTR system simulation under different initial conditions and control inputs; record the state variables (reactant concentration and temperature) and control variables (coolant temperature) at each time step; pre-process the collected data, including denoising and normalization; organize the processed data into a data set according to the time series. The purpose of this step is to obtain sufficient system dynamic information to provide data support for the subsequent Koopman operator approximation. The data set can be expressed as:

[0134]

[0135] Among them, x k is the state vector at time k, u k is the control input at time k, x k+1is the state vector at time k+1, and M is the total number of sampling points. In this way, the dynamic behavior data of the CSTR system under different working conditions can be obtained, providing a basis for subsequent model identification.

[0136] The specific implementation of step S40 is: using the extended dynamic mode decomposition (EDMD) method to solve the finite-dimensional approximation of the Koopman operator and obtain the state space matrix of the high-dimensional linear CSTR system. This process includes the following sub-steps: constructing data matrices X and Y, where X contains the state lifting vector ψ(x k ), Y contains the state improvement vector ψ(x k+1 ); Calculate the matrix G = X T X and A = X T Y; solve the linear equation system GA = A, and obtain the finite-dimensional approximate matrix A of the Koopman operator; calculate the control matrix B, which can be solved by the least squares method; construct the output matrix C so that the original state can be reconstructed from the high-dimensional state; verify the accuracy of the obtained high-dimensional linear model, which can be evaluated by comparing the model prediction value and the actual system response. The purpose of this step is to obtain a high-dimensional linear representation of the CSTR system, laying the foundation for subsequent model predictive control. The state space representation of the high-dimensional linear system is:

[0137] z k+1 =Az k +Bu k

[0138] x k =Cz k

[0139] Among them, z k is the high-dimensional state vector at time k, A is the state transfer matrix, B is the control matrix, and C is the output matrix. In this way, the nonlinear CSTR system is transformed into a high-dimensional linear system, which facilitates the application of mature linear control theory.

[0140] The specific implementation of step S50 is: based on the obtained high-dimensional linear model, establish a prediction equation for the future state of the CSTR system. This process includes the following sub-steps: determining the prediction time domain N p , usually 10 to 30 sampling periods are selected; determine the control time domain N c , generally less than or equal to N p , you can choose 5 to 15 sampling periods; use high-dimensional linear models to recursively calculate the future N p The state prediction value at each moment; construct the prediction equation matrix form, and pThe state prediction value at each moment is expressed as a function of the current state and future control input; system constraints are considered, such as the upper and lower limits of state variables and control variables; and a disturbance model is introduced to improve the prediction accuracy. The purpose of this step is to provide a prediction basis for subsequent model predictive control. The prediction equation can be expressed as:

[0141] Z=Gx k +ΦU

[0142] Among them, Z is the future N p The high-dimensional state prediction vector at the moment, x k is the current state, U is the future state c In this way, the future behavior of the CSTR system can be predicted based on the current state and future control inputs.

[0143] The specific implementation of step S60 is: design the objective function of the model predictive controller, including state tracking error, control input and terminal penalty term. This process includes the following sub-steps: define the state tracking error term to measure the deviation between the predicted state and the reference trajectory; define the control input term to limit the amplitude and rate of change of the control quantity; define the terminal penalty term to ensure the stability of the closed-loop system; introduce a weight matrix to balance the importance of each term; consider the system constraints and convert them into penalty terms or hard constraints of the objective function; construct a complete objective function expression. The purpose of this step is to form a comprehensive performance indicator to guide the optimization process of the controller. The objective function can be expressed as:

[0144]

[0145] Among them, z k+i|k is the state at time k+i predicted at time k, r k+i is the reference trajectory, u k+i is the control input, Q, R, and P are the state weight matrix, control weight matrix, and terminal penalty matrix, respectively. By adjusting these weight matrices, the control performance and control cost can be balanced.

[0146] The specific implementation method of step S70 is: convert the objective function into a standard linear quadratic programming problem and solve the optimal control sequence. This process includes the following sub-steps: expand the objective function to obtain a quadratic expression for the control sequence U; construct the Hessian matrix H and gradient vector f of the quadratic programming problem; convert the system constraints into the form of inequality constraints AU≤b; select an appropriate quadratic programming solution algorithm, such as the interior point method or the effective set method; set the parameters of the solver, including the maximum number of iterations, convergence threshold, etc.; solve the quadratic programming problem and obtain the optimal control sequence. The purpose of this step is to efficiently solve the control optimization problem and obtain the optimal control strategy. The standard form of the quadratic programming problem can be expressed as:

[0147]

[0148] stAU≤b

[0149] Among them, U is the control sequence, H is the Hessian matrix, f is the gradient vector, and A and b define the constraints. By solving this optimization problem, we can get the control sequence that minimizes the objective function while satisfying the constraints.

[0150] The specific implementation of step S80 is: implement the rolling optimization strategy, update the system state at each sampling time, and repeat steps S50 to S70. This process includes the following sub-steps: at the current time k, obtain the actual system state x k Based on the current state, use the prediction equation of step S50 to predict the future N p The system behavior at each moment; construct an optimization problem using the objective function designed in step S60; solve the optimization problem by the method of step S70 to obtain the control sequence U k ; Only execute control sequence U k The first element u k ; Wait for the next sampling time k+1; Repeat the above process. The purpose of this step is to achieve closed-loop control and continuously adjust the control strategy according to the latest system status. The rolling optimization strategy can be expressed as:

[0151] for k=1,2,...

[0152] Get the actual status x k

[0153] Predict the future state Z = f(x k , U)

[0154] Solving optimization problems

[0155] Execute control input u k

[0156] end

[0157] In this way, real-time optimized control of the CSTR system can be achieved to adapt to changes in system status and external disturbances.

[0158] The specific implementation method of step S90 is: by adjusting the state weighting matrix Q and the control weighting matrix R, the controller performance is optimized to achieve high-precision trajectory tracking control of the CSTR system. This process includes the following sub-steps: setting the initial Q and R matrices; running a closed-loop control simulation, recording the system response and control input; calculating performance indicators, such as the root mean square value of the tracking error and the energy consumption of the control input; adjusting the elements of the Q and R matrices according to the performance indicators; repeating the simulation and adjustment process until satisfactory control performance is achieved; verifying the robustness of the resulting controller under different operating conditions. The purpose of this step is to fine-tune the controller parameters to obtain the best control effect. The optimization process can be expressed as:

[0159]

[0160] stQ>0,R>0

[0161] Where RMSE(e) is the root mean square value of the tracking error, Energy(u) is the energy consumption of the control input, and λ is the balancing factor. In this way, a balance can be achieved between tracking accuracy and control cost to achieve high-performance control of the CSTR system.

[0162] The specific implementation of the present invention is explained in more detail below in conjunction with the specific continuous stirred reactor system structure:

[0163] Depend on Figure 2 It can be seen that the material is at a certain concentration C A0 , temperature T0, flow rate Q enters the reactor with volume V, and a first-order irreversible exothermic reaction occurs in the reactor: A→B, and the output concentration C is A , flow rate Q, temperature T of the material, the concentration here is for the reactant A. The reactor is provided with a cooling layer outside, and a cold (hot) agent can be introduced to adjust the temperature of the reactor to meet the reaction needs. The flow rate of the coolant in the cooling layer is Q C , temperature is T c , the following assumptions need to be made when building the CSTR model:

[0164] 1) The materials in the reactor are completely mixed;

[0165] 2) The total volume of materials in the reactor remains unchanged before and after the reaction;

[0166] 3) The reference value does not change during the reaction;

[0167] Table 1 Meaning of symbols in CSTR system

[0168]

[0169]

[0170] Based on the assumptions and the law of conservation of matter and energy, the kinetic model of the continuous stirred reactor is as follows:

[0171]

[0172]

[0173] The meaning of each parameter in formula (1) is shown in Table 1.

[0174] To simplify the above model (Formula (1)), let x1 = C A ,x2=T,u=T C , Then formula (1) can be rewritten as:

[0175]

[0176] Among them, x1 and x2 are system state variables, and y is the system output.

[0177] Consider a discrete-time nonlinear dynamical system

[0178] x(k+1)=f(x(k),u(k)),x(k)∈M. (3)

[0179] Where x(k)=[x1(k)x2(k)x α (k)] Τ , α is the system dimension, k represents the time step, and f represents the nonlinear mapping of different states in the state space M. Assuming that the observation function of the system is g(x), the Koopman operator κ is defined as a linear operator acting on the observation function g:

[0180] κg(x(k))=g(f(x(k),u(k)))=g(x(k+1)). (4)

[0181] From the definition, we can know that the Koopman operator updates the value of the observation function g by evolving the trajectory of g(x) in the state space. The evolution process is as follows: Figure 3 shown.

[0182] Since the Koopman operator is infinite-dimensional, it is necessary to approximate the infinite-dimensional Koopman operator with finite dimensions in practical applications. The finite-dimensional approximation is to approximate the Koopman operator by selecting a finite number of basis functions in the state space and represent it as a finite-dimensional matrix. Extended Dynamic Mode Decomposition (EDMD) is a method for extracting the dynamic characteristics of a system. It is mainly used to approximate the Koopman operator from finite-dimensional data and analyze the dynamic behavior of the system.

[0183] For the nonlinear dynamic system, that is, equation (3), the state improvement function is defined as follows:

[0184] Ψ(x(k), u(k))=[Ψ(x(k)), Ψ(u(k))] T . (5)

[0185] The promotion status is:

[0186] in, is the basis function, and β is the system dimension in the new space.

[0187] The Koopman operator is used to establish a finite high-dimensional linear dynamic model in the state space. The constructed high-dimensional linear model is as follows:

[0188] z(k+1)=Az(k)+Bu(k),

[0189]

[0190] In the formula, is the state of the high-dimensional system, z(k)=[z1(k) z2(k) … z β (k)] T , is an estimate of the original nonlinear system state x(k). are linear constant matrices, m is the number of columns of matrix B, and n is the number of rows of matrix C.

[0191] According to the true state and control input data of the original nonlinear system, i.e., equation (3), the data set is defined as follows:

[0192] X=[x(1),…,x(P)],

[0193] Y=[x(2),…,x(P+1)],

[0194] U=[u(1),…,u(P)]·(8)

[0195] Among them, x(k+′)=f(x(k), u(k)), k=1, 2,…,.

[0196] Based on the state improvement function (Formula (5)), the datasets X and Y are improved and are X lift and Y lift :

[0197] X lift =[Ψ(x(1)),…,Ψ(x(P))],

[0198] Y lift =[Ψ(x(2)),…,Ψ(x(P+1))]. (9)

[0199] According to the EDMD method, A, B, and C in formula (7) can be obtained by solving the following optimization problem:

[0200]

[0201]

[0202] The solution of formula (10) is:

[0203] [A, B] = Y lift [X lift , U] ↑ ,

[0204]

[0205] Among them, is the Moore-Penrose pseudoinverse.

[0206] Substituting A, B, and C obtained in equation (11) into equation (7), we can obtain the finite-dimensional approximate model of the Koopman operator in the state space.

[0207] For the CSTR system in equation (1), the state variable x(k) = [C A T],C A and T represent the discharge concentration and discharge temperature of the continuous stirred reactor respectively; the control variable u(k)=[T C ], T C The goal of trajectory tracking control of CSTR is to design a controller so that the CSTR system can track the desired trajectory in real time under the influence of the control input; that is, under effective control, the actual output (including concentration and temperature) of CSTR at any time can be infinitely close to the expected value.

[0208] To obtain the Koopman high-dimensional linear model, the thin plate spline radial basis function is selected as the basis function: x1=C A, x2=T, according to the Koopman operator defined in formula (4), the specific application form of the CSTR system in the Koopman operator method is:

[0209]

[0210] The state improvement function is:

[0211]

[0212] The basis function is

[0213]

[0214] Among them, c i is the center point of the i-th basis function. If the center point is [i,i] Τ but:

[0215]

[0216] Collect the data of the evolution of the CSTR system status over time as follows:

[0217]

[0218]

[0219]

[0220] According to formulas (9), (10), and (11), the state space matrix of the high-dimensional linear CSTR system can be obtained.

[0221] The Koopman model of CSTR is shown in formula (7), and the prediction equation for the future state of the high-dimensional system z(k) can be written as:

[0222] z(k+1|k)=Az(k)+Bu(k),

[0223] z(k+2|k)=Az(k+1|k)+Bu(k+1|k),

[0224]

[0225] z(k+N C |k)=Az(k+N C -1|k)+Bu(k+N C -1|k),

[0226]

[0227] z(k+N P |k)=Az(k+NP -1|k)+Bu(k+N P -1|k). (17)

[0228] Where, N C To control the time domain, N P is the prediction time domain, z(k+i|k) and u(k+i|k) represent the predicted state and control input at time k respectively.

[0229] Define y(k) as the output of the CSTR system. From formula (2), we can see that the output of the CSTR system is the state of the system, that is, y(k) = x(k). According to formula (7), the predicted output is expressed as:

[0230]

[0231] Based on formula (18), the prediction equation of the system can be rewritten as:

[0232]

[0233] in:

[0234]

[0235]

[0236]

[0237]

[0238] The reference trajectory is usually expressed in the form of a time series, containing the expected value of the system state or output at each time step. The reference trajectory can be expressed as:

[0239]

[0240] Among them, y ref (k+i) is the expected state or output value at the i-th time step.

[0241] The present invention adopts the following objective function:

[0242]

[0243] stz(k+1)=Az(k)+Bu(k),

[0244]

[0245] k=1,...,N P -1,

[0246] x(k+i)∈Ω.(21)

[0247] Among them: Q, R and F are the state weight matrix, control weight matrix and terminal penalty matrix of the nonlinear system, respectively, and all are positive definite weight matrices.

[0248] In order to improve the solution efficiency, the performance index is transformed into:

[0249]

[0250] Where:

[0251]

[0252]

[0253]

[0254]

[0255] Formula (21) can be solved as a standard linear quadratic programming problem, so that the global optimal solution can be quickly obtained. By solving formula (21), we can get N in the prediction time domain P The optimal control sequence U(k) in . The first element in U(k) is the coolant temperature at time k.

[0256] The specific controller design block diagram is as follows Figure 3 In practical applications, it is necessary to adjust the Q and R matrices so that the system can achieve high-precision trajectory tracking within the specified constraints.

[0257] 3 Stability analysis

[0258] Assumptions:

[0259] (1) The initial state x(0) is feasible

[0260] (2) There is an optimal control sequence at each sampling time:

[0261] U(k)=[u(k),u(k+1),…,u(k+N-1)] Τ . (twenty four)

[0262] (3) At each sampling moment, the corresponding expected trajectory in the future can be obtained:

[0263]

[0264] The control time domain and prediction time domain are defined to be the same, that is, N P =N C , the terminal equality constraint is:

[0265]

[0266] The Lyapunov function of the system is defined as:

[0267]

[0268] It can be obtained that V(k)≥0, then V(k)≥0 is positive definite.

[0269]

[0270] because and Always greater than or equal to 0.

[0271] Therefore, V(k+1)≤V(k), and stability is proved.

[0272] In order to verify the feasibility and effectiveness of the high-dimensional Koopman linear model of the continuous stirred reactor proposed in the present invention, the present invention built a continuous stirred reactor model predictive control system based on the Koopman operator based on MATLAB / simulink. In order to further demonstrate the superiority of the CSTR model predictive controller designed by the present invention, the prediction results of the model of the present invention are compared with the local linearization model predictive control and the nonlinear model predictive control during simulation. The present invention uses the root mean square error (RMSE) to calculate the prediction error, and the formula is as follows:

[0273]

[0274] Table 2 CSTR model parameters

[0275]

[0276] The model parameters of the continuous stirred reactor of the present invention are shown in Table 2. Based on the Runge-Kutta method, the dynamic model (2) of the continuous stirred reactor is discretized, and the data of the dynamic system of the continuous stirred reactor is collected. The sampling interval T s =0.1s. 1000 initial points were selected, and 2500 samples were taken from each initial point, that is, 1000 trajectories were selected for sampling, and each trajectory was generated by running the CSTR kinetic model for 250 seconds. In order to ensure that the collected system data can more comprehensively reflect the system characteristics, the selected initial points should be different. In order to ensure that the data obtained from the simulation experiment is within a reasonable physical range, the state quantities x1 and x2 of the present invention are randomly generated within [0, 0.9] and [0, 170]. The control quantity u is randomly generated within [-50, 50].

[0277] To obtain the Koopman high-dimensional linear model, the thin plate spline radial basis function is selected as the basis function: x0 is a constant vector whose elements are uniformly distributed in a given interval, and the center point of the radial basis function is randomly selected in [0,20]. In data collection, the system state is processed by using the thin plate spline radial basis function to increase the dimension, and the collected state data can accurately reflect the nonlinear dynamics of the system. In the process of increasing the dimension, the thin plate spline radial basis function maps the original state to a high-dimensional space, so that the state after the dimension increase has better linear characteristics. The Koopman operator is constructed by extending the dynamic mode decomposition method, and the fitting accuracy of the model is improved by using the characteristics of the thin plate spline radial basis function, so that the prediction model can accurately describe the dynamic behavior of the system. Finally, the original continuous stirred reactor dynamics system is increased to a 24-dimensional linear space, that is, a 24-dimensional continuous stirred reactor linear model is obtained.

[0278] When the initial value of the state quantity is selected as [20,140], two different trajectories are used to perform trajectory tracking simulation verification. Trajectory I is a constant signal, and trajectory II is a square wave signal. The controller parameters of the entire simulation process are set consistently. The model predictive control (KMPC) method based on the Koopman operator of the present invention is compared with the local linearized model predictive control (LMPC) method in the literature and the nonlinear model predictive control (NMPC) method in the literature to verify the superiority of the method of the present invention. The relevant parameters of the controller are shown in Table 3.

[0279] Table 3 Model predictive controller related parameters

[0280]

[0281] Figures 5 to 7 This is the simulation result of the tracking performance of trajectory I. Figure 5 The three control methods can make the CSTR discharge concentration converge to the expected value smoothly without overshoot, but the KMPC method proposed in this invention has the fastest tracking speed. Figure 6 It is the control effect of CSTR discharge temperature. Obviously, the KMPC method proposed in the present invention has a significant improvement in convergence speed compared with the other two methods (LMPC method and NMPC method), and the convergence time is reduced by about one order of magnitude. Figure 7 Table 4 and Table 5 are the root mean square errors of trajectory I under the three control methods. Although the error of the method of the present invention is high when the discharge concentration tracks trajectory I from 0 to 100 s, the difference is not large. After 100 s, the error is significantly better than that of the other two methods, and the overall control effect is more accurate.

[0282] Table 4 Root mean square error of discharge concentration trajectory I under three control methods

[0283]

[0284] Table 5 Root mean square error of discharge temperature trajectory I under three control methods

[0285]

[0286] Figures 8 to 10 This is the simulation result of the tracking performance of trajectory II. Figure 8 It is the control effect of CSTR discharge concentration. Obviously, the KMPC method proposed in this invention shows the best tracking effect, the fastest tracking speed and the smallest steady-state error; while the LMPC method and the NMPC method show serious lag in tracking performance when the concentration expected value jumps. Fig. 9 It is the control effect of CSTR discharge temperature. Although the three methods can track trajectory II and tend to be stable, the KMPC method proposed in this invention has a significant improvement in steady-state performance. Fig.10 The table 6 and table 7 show the root mean square error of trajectory II under the three control methods. The root mean square error of the discharge concentration and discharge temperature tracking trajectory II of the method of the present invention in 400s-1200s is significantly lower than that of the other two methods, but the control effect of the method of the present invention in the early stage of the discharge temperature needs to be further optimized. In general, the control effect of the method of the present invention is more accurate over a longer time range.

[0287] Table 6 Root mean square error of discharge concentration trajectory II under three control methods

[0288]

[0289] Table 7 Root mean square error of discharge temperature trajectory II under three control methods

[0290]

[0291] Table 8 KMPC controller running time

[0292]

[0293] Comprehensive analysis of the table data shows that the average running time of the KMPC controller is 0.007 seconds under trajectory I and 0.0167 seconds under trajectory II. Although the data collection and processing time is long, the average running time of each cycle is still within an acceptable range, showing the high computational efficiency of KMPC in dealing with dynamic changes. Although the data processing time of each stage of trajectory II is long, the running time of MPC itself is short, indicating that it has the potential for real-time operation in practical applications, and is particularly suitable for industrial control scenarios that require fast response.

[0294] In summary, although the other two control methods and the control method of the present invention can ensure the stability of the closed-loop control system, the method of the present invention exhibits a faster convergence speed, higher control accuracy, and has the best dynamic performance and steady-state performance. It can be seen that the high-dimensional linear model established by Koopman operator identification has better accuracy and can capture the key nonlinearity of the actual CSTR system, thereby avoiding possible deviations and infeasible solutions of the set values ​​under MPC conditions.

[0295] Specifically, the principle of the present invention is:

[0296] 1. High-precision nonlinear system characterization:

[0297] The present invention uses Koopman operator theory to map the CSTR system from a low-dimensional nonlinear space to a high-dimensional linear space. This mapping can capture the nonlinear dynamic characteristics of the system without explicit nonlinear modeling. By selecting an appropriate set of basis functions, the complex dynamic behavior of the CSTR system can be more accurately characterized in a high-dimensional space. Compared with traditional linearization methods, this method can maintain the accuracy of the model in a larger operating range, thereby providing a more reliable basis for subsequent controller design.

[0298] 2. Improved computing efficiency:

[0299] After converting the nonlinear system into a high-dimensional linear system, the control problem can be converted into a standard linear quadratic programming problem. This conversion significantly reduces the computational complexity of the solution process. Compared with directly dealing with nonlinear optimization problems, linear quadratic programming problems can be quickly solved using mature numerical optimization algorithms, greatly improving computational efficiency. This feature makes this method more suitable for real-time control applications, and can complete the calculation of the control quantity within a shorter sampling period, thereby improving the system's response speed and control accuracy.

[0300] 3. Strengthen model prediction capabilities:

[0301] The present invention combines the Koopman operator method with the model predictive control (MPC) strategy, giving full play to the advantages of both. The Koopman operator provides an accurate prediction of the future state of the system, while MPC uses this prediction information to optimize the control sequence. This combination not only improves the accuracy of the prediction, but also can effectively deal with system constraints and multi-objective optimization problems. Compared with the traditional MPC method, this method shows better prediction performance and control effect when dealing with strongly nonlinear systems.

[0302] 4. Enhance system robustness:

[0303] The method based on the Koopman operator has a strong data-driven characteristic and can learn and extract dynamic features from the system operation data. This characteristic makes the controller have good adaptability to system parameter changes and external disturbances. Even when there is uncertainty in the system model or slow changes occur, the controller can still maintain stable performance. In addition, by updating the approximation of the Koopman operator online, the system's anti-interference ability and adaptability can be further improved.

[0304] 5. Improve control performance:

[0305] The method of the present invention can achieve faster transient response and smaller steady-state error of the CSTR system. Since the Koopman operator can accurately capture the nonlinear dynamic characteristics of the system, the controller can more accurately predict and compensate for the nonlinear behavior of the system. This enables the system to reach the target state faster and achieve more accurate regulation near the target state.

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

Claims

1. A continuous stirred reactor control method based on Koopman operator, characterized in that: The following steps are involved: S10, establish a dynamic model of the continuous stirred reactor and determine the state variables and control variables of the system; S20, selecting an appropriate set of basis functions, defining a state lifting function, and mapping the original nonlinear CSTR system to a high-dimensional linear space; S30, collect data on the evolution of the CSTR system state over time and construct a data set; S40, using the extended dynamic mode decomposition method, solve the finite-dimensional approximation of the Koopman operator and obtain the state space matrix of the high-dimensional linear CSTR system; S50, based on the obtained high-dimensional linear model, establish a prediction equation for the future state of the CSTR system; S60. Design the objective function of the model predictive controller, including state tracking error, control input and terminal penalty terms; S70, transforming the objective function into a standard linear quadratic programming problem and solving the optimal control sequence; S80, implement the rolling optimization strategy, update the system status at each sampling time, and repeat steps S50 to S70; S90, by adjusting the state weighting matrix and the control weighting matrix, the controller performance is optimized to achieve predictive control of high-precision trajectory tracking of the CSTR system.

2. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 1, characterized in that: The step S10 specifically includes: Step 101, determining the state variables of the continuous stirred tank reactor system, including reactant concentration and reactor temperature; Step 102, determining the control variable of the continuous stirred tank reactor system, that is, the coolant temperature; Step 103, listing the differential equation of reactant concentration, taking into account input concentration, reaction consumption and discharge loss; Step 104, listing the differential equation of the reactor temperature, taking into account the input heat, reaction heat, cooling heat and discharge heat; Step 105, introducing relevant parameters, including reaction rate constant, activation energy, reaction heat, flow rate, reactor volume, feed concentration, feed temperature, gas constant, density, specific heat capacity, heat transfer coefficient and heat transfer area; Step 106: Normalize the differential equation to obtain a dimensionless dynamic model.

3. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 2, characterized in that: The step S20 specifically includes: Step 201, selecting a thin plate spline radial basis function as a basis function; Step 202, determine the number of basis functions through cross-validation, and select 10 to 20 basis functions; Step 203: define a state improvement function to map the original state space to a high-dimensional feature space; Step 204, determining the center point of the basis function by uniform sampling or clustering method; Step 205: construct a state improvement function, take the original state variables and the control variables as input, and output a high-dimensional feature vector; Step 206: Verify the validity of the constructed state improvement function to ensure that it can accurately capture the nonlinear characteristics of the system; Step 207: Convert the nonlinear system into a high-dimensional linear system.

4. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 3, characterized in that: The step S30 specifically includes: Step 301, set the sampling time interval, and select a sampling period between 0.1 seconds and 1 second according to the dynamic characteristics of the system; Step 302, determine the total sampling time, and select 5 to 10 times of the system response time to fully capture the system dynamics; Step 303, running a continuous stirred tank reactor system simulation under different initial conditions and control inputs; Step 304, recording the state variables of each time step, including reactant concentration and reactor temperature, and the control variable, i.e., coolant temperature; Step 305: preprocessing the collected data, including denoising and normalization; Step 306: Arrange the processed data into a data set according to time series.

5. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 4, characterized in that: The step S40 specifically includes: Step 401, construct a data matrix, including the state improvement vectors at all times; Step 402, calculate matrix operations to obtain the intermediate matrix required for Koopman operator approximation; Step 403, solving the linear equations to obtain a finite-dimensional approximate matrix of the Koopman operator; Step 404, calculating the control matrix by the least square method; Step 405: construct an output matrix so that the original state can be reconstructed from the high-dimensional state; Step 406 : Verify the accuracy of the obtained high-dimensional linear model by comparing the model prediction value with the actual system response.

6. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 5, characterized in that: The step S50 specifically includes: Step 501, determine the prediction time domain, and select 10 to 30 sampling periods; Step 502, determine the control time domain, select 5 to 15 sampling periods, and do not exceed the prediction time domain; Step 503: using a high-dimensional linear model, recursively calculate the state prediction value in the future prediction time domain; Step 504: construct a prediction equation matrix form, and express the state prediction value in the future prediction time domain as a function of the current state and the future control input; Step 505: Consider system constraints, including upper and lower limits of state variables and control variables; Step 506: Introduce an interference model to improve prediction accuracy.

7. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 6, characterized in that: The step S60 specifically includes: Step 601: define a state tracking error term to measure the deviation between the predicted state and the reference trajectory; Step 602: define control input items for limiting the amplitude and change rate of the control amount; Step 603: define a terminal penalty term to ensure the stability of the closed-loop system; Step 604: introduce a weight matrix to balance the importance of state tracking error, control input and terminal penalty; Step 605: Consider system constraints and convert them into penalty terms or hard constraints of the objective function; Step 606: Construct a complete objective function expression.

8. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 7, characterized in that: The step S70 specifically includes: Step 701: Expand the objective function to obtain a quadratic expression for the control sequence; Step 702: construct the Hessian matrix and gradient vector of the quadratic programming problem; Step 703: convert the system constraints into the form of inequality constraints; Step 704: Select an appropriate quadratic programming solution algorithm, such as an interior point method or an effective set method; Step 705: Setting the parameters of the solver, including the maximum number of iterations and the convergence threshold; Step 706: Solve the quadratic programming problem to obtain the optimal control sequence.

9. The method for controlling a continuous stirred reactor based on a Koopman operator according to claim 8, characterized in that: The step S80 specifically includes: Step 801, obtaining the actual state of the continuous stirred tank reactor system at the current moment; Step 802: Based on the current state, use the prediction equation to predict the system behavior in the future prediction time domain; Step 803: construct an optimization problem using the designed objective function; Step 804, solving the optimization problem to obtain a control sequence; Step 805, only the first element of the control sequence is executed as the current control input; Step 806: Wait for the next sampling moment; Step 807, repeating steps 801 to 806 to achieve rolling optimization control of the continuous stirred tank reactor system.

10. The method for controlling a continuous stirred reactor based on Koopman operator according to claim 9, characterized in that: The step S90 specifically includes: Step 901, setting the initial state weight matrix and control weight matrix; Step 902, running a closed-loop control simulation, recording the continuous stirred tank reactor system response and control input; Step 903, calculating performance indicators, including the root mean square value of the tracking error and the energy consumption of the control input; Step 904: Adjust the elements of the state weight matrix and the control weight matrix according to the performance index; Step 905, repeat the simulation and adjustment process until satisfactory control performance is achieved; Step 906, verifying the robustness of the obtained controller under different working conditions; Step 907: Determine the final state weighted matrix and control weighted matrix to achieve high-precision trajectory tracking control of the continuous stirred tank reactor system.

Images