Double-tower system process optimization method based on LU decomposition

Abstract

The invention discloses a double-tower system process optimization method based on LU decomposition, and the method comprises the steps: taking the total value of a liquid phase product produced by a double-tower system formed by a depropanizing tower and a debutanizing tower as an optimization target, and constructing an optimization model by combining a constraint equation and a normalization equation; and performing iterative solution on the optimization model based on LU decomposition and a quadratic sequence programming method. According to the method, by judging the correlation of the constraint equation and optimizing the condition number, the influence of ill-conditioned problems can be effectively avoided, and the numerical stability in the solving process is remarkably improved. According to the method, the redundant constraint equation is removed, and the condition number is improved through LU decomposition basis changing, so that unnecessary calculation can be reduced, and the convergence speed of iteration is improved.

Description

Technical Field

[0001] The invention belongs to the technical field of large-scale chemical process optimization and relates to a double-tower system process optimization method based on LU decomposition. Background Art

[0002] Depropanizer and debutanizer are two important distillation operation units in the ethylene production process. The two towers exchange materials and energy with each other. In order to keep the entire system in the optimal operating state, it is necessary to combine the two distillation processes of depropanizer and debutanizer and study them from a system perspective. Since the depropanizer and debutanizer have strong coupling before and after, the optimization solution is more complicated than that of a single tower. Even if the single tower is optimized very well, the combined tower system may not be able to be solved. The reason is that the combined models are interrelated and affect each other, which can easily lead to related or contradictory constraint equations, and the problem may become unsolvable at this time. In addition, the dimensions of the combined distillation towers reach thousands, and there are many constraint equations. The condition number of the constraint matrix directly affects the well-posedness of the problem. If the condition number is poor, it means that the derivative matrix is ​​close to singularity, the solution will become unstable, the numerical error will be amplified, and then the accuracy of the final solution will be affected, and it may even be impossible to solve. In addition, the large scale of the problem also leads to complex calculations. The optimization of the entire process involves a large number of mass balance and energy balance equations. In addition, it is also necessary to consider the phase equilibrium equations in the tower, such as Raoult's law or Gibbs free energy model, and the dynamic behavior of the tower, such as reflux ratio, product purity, etc. All these physical and chemical processes and optimization constraints can be described by a set of complex nonlinear equations, and the variables and equation dimensions involved are large, but compared with a single tower, the model structure of the dual tower is more obvious and the sparsity is stronger.

[0003] LU decomposition is a method of decomposing a matrix into the product of a lower triangular matrix and an upper triangular matrix. Specifically, if an n*n square matrix A can be decomposed into a lower triangular matrix L and an upper triangular matrix U, that is, P*A=L*U, then this decomposition of A is called LU decomposition, and the P matrix is ​​called a permutation matrix. LU decomposition mainly uses the Doolittle algorithm and the Crout algorithm.

[0004] The Doolittle algorithm constructs the L matrix and the U matrix as follows: Let matrix A be an m*n matrix, L be an m*m unit lower triangular matrix with the main diagonal elements being 1, and U be an m*n upper triangular matrix, so LU is also called a triangular decomposition of the A matrix, where the three matrices A, L, and U are in the form of:

[0005]

[0006] P*A=L*U

[0007] From the above formula, comparing the terms in matrix A and LU, we can get the general term:

[0008]

[0009] Therefore, the expression of the upper triangular matrix U is obtained as:

[0010]

[0011] When u jj = 0, directly set l ij =0(i>j).

[0012] The Crout algorithm and the Doolittle algorithm differ in how to construct the lower triangular matrix L and the upper triangular matrix U. The diagonal of the L matrix constructed by the Crout algorithm is not 1, but the diagonal elements of the upper triangular matrix U constructed by the Crout algorithm are 1. The difference in the calculation formula from the Doolittle algorithm is:

[0013]

[0014] LU decomposition can retain and utilize the sparse structure of the matrix during the decomposition process, thereby reducing the amount of calculation and storage requirements, significantly improving the efficiency of processing sparse matrices; it is particularly suitable for processing large-scale sparse matrices. Therefore, the present invention proposes a dual-tower system process optimization method based on LU decomposition to achieve coupled optimization of the depropanizer and debutanizer. Summary of the invention

[0015] The purpose of the present invention is to provide a double-tower system process optimization method based on LU decomposition.

[0016] The invention specifically comprises the following steps:

[0017] Step 1, establish the equilibrium equations of the depropanizer and debutanizer, and normalize the equations;

[0018] The depropanizer and debutanizer are both distillation towers, and the steady-state process model of the distillation tower satisfies the material balance, energy balance, and vapor-liquid phase reaction balance equations, while the molar fraction composition of the vapor-liquid phase satisfies the normalization equation.

[0019] Assuming that the number of plates in the distillation tower is M, the condenser and reboiler can be regarded as two upper and lower plates in the modeling process, then the actual number of plates is N = M + 2. There is a feed, a vapor phase side discharge and a liquid phase side discharge on each plate. Under steady-state operation, because there is no material accumulation in the plate, the amount of material entering a plate is equal to the amount of material leaving the plate. For each plate, the sum of the total molar flow of the liquid phase and the gas phase entering the plate is equal to the sum of the total molar flow of the liquid phase and the gas phase leaving the plate.

[0020] The material balance equation for the first plate is obtained as follows:

[0021] V 2 y i,2 -V 1 y i,1 -(S 1 +L 1 )x i,1 =0;

[0022] Among them, V 2 is the gas phase flow rate leaving the first tray, y i,2 is the mole fraction of component i in the gas phase leaving the first tray, V 1 is the gas phase flow rate entering the first tray, y i,1 is the mole fraction of the i-th component in the gas phase entering the first tray, S 1 is the liquid phase flow rate entering the first tray, L 1 is the liquid phase flow rate leaving the first tray, x i,1 is the mole fraction of the i-th component in the liquid phase leaving the first tray.

[0023] The material balance equation for the jth tray is:

[0024] F j z i,j +L j-1 x i,j-1 +V j+1 y i,j+1 -(L j +S j )x i,j -(V j +G j )y i,j =0,(2≤j≤N-1);

[0025] Among them, F j is the liquid phase flow rate entering the jth tray, z i,j The mole fraction of the i-th component in the liquid phase entering the j-th tray, L j-1 Liquid flow rate leaving from the j-1th tray, x i,j-1 The mole fraction of the i-th component in the liquid phase leaving the j-1-th plate, V j+1 The gas phase flow rate entering the jth tower plate. i,j+1 The mole fraction of the i-th component in the gas phase entering the j+1-th tray, L j The liquid phase flow rate leaving the jth tray, S j The side liquid phase flow rate leaving from the jth tray, G j Side gas flow rate leaving from the jth tray, y i,j is the mole fraction of the i-th component in the gas phase entering the j-th tower plate.

[0026] The material balance equation for the Nth plate is obtained as follows:

[0027] L N-1 x i,N-1 -(V N +G N )y i,N -L N x i,N =0;

[0028] Among them, L N-1 is the liquid phase flow rate leaving from the N-1th tray, x i,N-1 is the mole fraction of the i-th component in the liquid phase leaving the N-1th plate, V N is the gas phase flow rate leaving from the Nth plate, G N is the side gas flow rate leaving from the Nth tray, y i,N is the mole fraction of the i-th component in the gas phase leaving the N-th plate. N is the liquid phase flow rate leaving from the Nth plate, x i,N is the mole fraction of component i in the liquid phase leaving the Nth tray.

[0029] The energy balance equation for the first tray is obtained as follows:

[0030]

[0031] Among them, V 2 represents the gas phase flow rate leaving from the second tray, represents the specific enthalpy of the gas phase leaving the second tray, V 1 represents the gas phase flow rate entering the first tray 1, represents the specific enthalpy of the gas phase entering the first tray, L 1 represents the liquid phase flow rate leaving from the first tray, S 1 represents the liquid phase flow entering the first tray, represents the specific enthalpy of the liquid phase leaving the first tray, Q 1 Represents the heat added or removed on the first stage.

[0032] The energy balance equation for the jth plate is obtained as follows:

[0033]

[0034] Among them, F j represents the liquid phase flow rate entering the jth tray, represents the specific enthalpy of the liquid phase entering the jth tray, V j+1 represents the gas phase flow rate entering the j+1th tray, represents the specific enthalpy of the gas phase entering the j+1th tray, L j-1 represents the liquid phase flow rate leaving from the j-1th tray, represents the specific enthalpy of the liquid phase leaving the j-1th plate, Q j represents the heat added or removed at stage j, V j represents the gas phase flow rate leaving from plate j, G j represents the side gas flow rate leaving tower j, S j Side stream liquid phase leaving from column j. represents the specific enthalpy of the gas phase leaving the jth plate, represents the specific enthalpy of the liquid phase leaving from the j-th tray.

[0035] The energy balance equation for the Nth plate is obtained as follows:

[0036]

[0037] Among them, L N-1 represents the liquid phase flow rate leaving from the N-1th tray, is the specific enthalpy of the liquid phase leaving the N-1th plate, Q N is the heat added or removed at the Nth stage, V N is the gas phase flow rate leaving from the Nth plate, G N is the side gas flow rate leaving from the Nth tray, is the specific enthalpy of the gas phase leaving the Nth tray.

[0038] The vapor-liquid phase reaction equilibrium equation is obtained as follows: i,j =K i,j x i,j (1≤i≤C,1≤j≤N);

[0039] Among them, y i,j is the mole fraction of the i-th component in the gas phase of the j-th tray, x i,j The mole fraction of the i-th component in the liquid phase on the j-th tray; K i,j is the equilibrium constant, which represents the relative concentration relationship between the gas phase and the liquid phase of the i-th component on the tower plate j. It is usually a constant related to temperature and pressure. C is the total number of components in the system, and N is the total number of tower plates.

[0040] The molecular normalization equation is:

[0041] Among them, y i,j is the mole fraction of the i-th component in the gas phase of the j-th tray, x i,j is the mole fraction of the i-th component in the liquid phase of the j-th tower plate, and C is the total number of components in the system.

[0042] Step 2: Establish an optimization model;

[0043] The total value of the liquid products produced by the depropanizer and debutanizer is taken as the optimization target. The total value of the products is calculated based on the market prices of C3 and C4 products and the outputs of the depropanizer and debutanizer. This is used as a basis for optimizing the production process to establish the optimization objective function F. The objective function F, constraint equations, and boundary constraints are combined to establish an optimization model. The constraint equations are the three equilibrium equations and the normalization equation in step one.

[0044] Step 3: Solve the optimization model established in step 2;

[0045] We choose to use the quadratic sequence programming method SQP to solve the problem. In each iteration, we construct a quadratic programming subproblem to approximate the original nonlinear problem, and then solve this subproblem to obtain the search direction and step size; then we combine LU decomposition to process the constraint matrix, solve the correlation, and perform infeasible recovery. The specific solution steps are as follows:

[0046] Step 3.1: Preprocess the constraint equations:

[0047] The original model of the problem is converted into a standard form suitable for processing by a numerical solver, including determining the initial values ​​of the optimization variables and scaling the variables with large unit range differences.

[0048] Step 3.2, derivative of constraint equation:

[0049] Each tray in the double tower system has a corresponding equilibrium equation, which combines the material balance equation, energy balance equation, vapor-liquid phase reaction balance equation, and the normalized equation of the vapor-liquid phase molar fraction group involved. Each constraint calculates the partial derivative of each variable separately, and the calculated derivative matrix is ​​called the Jacobian matrix ceqGRAD(x):

[0050]

[0051] Among them, x represents all the variables involved in the twin tower optimization process: x 1 ,x 2 ,...,x n , c 1 ,c 2 ,...,c m is the set of constraint equations, is the partial derivative of the ith equilibrium equation with respect to the jth variable.

[0052] Step 3.3: Remove redundant constraints by determining whether the constraint equations are linearly related and consistent:

[0053] In each iteration, the rank of the Jacobian matrix of the constraint is indirectly determined by LU decomposition, and the correlation of the constraint is determined by determining the rank of the derivative matrix of the constraint equation, that is, the Jacobian matrix. After LU decomposition of the Jacobian matrix, the L matrix and the U matrix are obtained, L*U=lu(ceqGRAD), and then the rank of the matrix is ​​quickly determined according to the structure of the U matrix, thereby determining the independence of the constraint equation. If the constraint equations are linearly related and consistent, these redundant constraint equations can be removed; if they are related but inconsistent, the system will immediately determine that the problem has no solution.

[0054] Step 3.4, iterative solution, change basis to deal with singular matrix;

[0055] The SQP algorithm is called to iteratively solve the model, converting the complex nonlinear optimization problem into a series of quadratic programming sub-problems. In each iteration, the algorithm constructs a local quadratic programming sub-problem based on the current iteration point x. k The objective function F and the constraint c(x) are quadratically approximated; at the iteration point x k The QP subproblem is expressed as follows:

[0056]

[0057] Where, d k is the search direction vector, g(x k ) is the derivative of the optimization objective function F, J(x k ) T is the constraint c(x k ), the derivative Jacobian matrix of W(x k ) is the second derivative matrix of the Lagrangian function, which is in g(x k ) and J(x k ) T The Lagrangian function is a function that integrates the constraint and the objective function, and its form is: L(x,λ,v,π)=F(x)+λ T c i (x)+v T (xx U )-π T (xx L ), λ, v, π are the corresponding coefficients. L represents the lower bound of the variable x, x U Represents the upper bound of the variable x.

[0058] Call the existing algorithm to solve the equation and get the search direction d k , used to update the current iteration point x k , the update formula is: k+1 =x k+d k After updating the iteration point, check whether the optimization termination conditions are met, such as the improvement of the objective function, the excessive number of iterations, or the violation of the equilibrium equation constraints. k The solution depends on the Jacobian matrix of the constraints, so once the matrix becomes singular, it is impossible to continue solving.

[0059] For the above singular cases, the constrained Jacobian matrix J is decomposed by LU:

[0060]

[0061] Among them J B The Jacobian matrix corresponding to the basic variables is called the basic matrix, J N Represents the Jacobian matrix of the dependent variable. The corresponding basic variable is: x 1 ,x 2 ,x 3 ,,,x m After LU decomposition, some basic variables and dependent variables are replaced by permutation matrix to replace the basic variables that cause the matrix singularity and restore the solvability of the matrix. The re-partitioned basic variables are

[0062] Step 3.5: Determine whether the optimization conditions are met:

[0063] Determine whether the optimization conditions are met. If not, introduce a new x k+1 , repeat steps 3.2-3.5 to continue the next iteration solution; if satisfied, terminate the calculation and output the optimization results: the maximum benefit value of the dual-tower system and the values ​​of each operating variable when the maximum benefit is obtained.

[0064] Furthermore, the specific steps of removing redundant constraints in step 3.3 are:

[0065] Step A1, setting the mark error topDep to 1e-10;

[0066] Step A2: Convert the Jacobian matrix ceqGRAD obtained in step 3.2 into a sparse form and perform LU decomposition. ceqGRAD contains many zero elements. Sparse LU decomposition decomposes the sparse matrix into a lower triangular matrix L and an upper triangular matrix U, as well as the related permutation matrix P.

[0067] Step A3, determine the redundancy of the constraints; determine whether the constraints are relevant by determining whether the U matrix obtained by LU decomposition of the Jacobian matrix of the constraints has all-zero rows at the current iteration point. If a row of the U matrix is ​​all-zero at the current iteration point, it means that the constraints of the corresponding row are relevant, that is, the redundancy needs to be removed, and the consistency is determined in the next step. If there is no all-zero row, it means that there are no redundant constraints, then skip steps A4 and A5, and continue with the subsequent solution.

[0068] Step A4, determine the consistency of the constraints; determine whether the product of the L matrix and the constraint c matrix is ​​zero in the row corresponding to the redundant constraint. If the row corresponding to the redundant constraint is zero, it means that the redundant constraints are consistent. Delete this row of constraints directly to remove redundancy and continue to solve; if the product of the L matrix and the constraint c matrix in the row corresponding to the redundant constraint is not zero, it means that the redundant constraints are inconsistent, the constraint equations are contradictory, and can no longer be solved. The optimization process ends and a prompt indicating that there is no solution is output.

[0069] Step A5, remove redundant constraints and update the constraint system; find the index of the redundant constraints through the elementary permutation matrix P obtained by LU decomposition, empty the rows where the redundant constraints are located in the constraint matrix and Jacobian matrix, update the constraint system, remove redundant equations, retain independent and consistent equations, and continue the subsequent optimization solution.

[0070] Furthermore, the specific process of processing the singular matrix basis change described in step 3.4 is:

[0071] Step B1, initializing the condition number;

[0072] The condition number can be regarded as a measure of the degree to which a matrix is ​​close to being singular. When the condition number is large, it means that the determinant of the matrix tends to zero and the matrix tends to be singular. Let the condition number κ(J B ) is initialized to 1. Then in each iteration, the current iteration point x k The basic matrix J B Calculate its condition number. The calculation formula of the condition number is: κ(J B )=||J B ||·||J B -1 ||, where ||J B || indicates J B 2-norm of the matrix.

[0073] Step B2, determine whether the condition number is acceptable;

[0074] The judgment condition for changing the basis is: the calculated condition number κ(J B ) and the pre-set threshold κ threshlod For comparison, if the current matrix J B The condition number κ'(J B) is greater than the threshold and relative to the previous condition number κ(J B ) increases by the set multiple.

[0075] Step B3, LU decomposition for base change;

[0076] When the condition number is poor, the basis change operation is performed through LU decomposition, and the LU decomposition is performed on the constrained Jacobian matrix J. The relationship between the basic matrix and the Jacobian matrix is: Among them J B represents the basic matrix, J N Represents the Jacobian matrix of the dependent variable. The corresponding basis variable is: x 1 ,x 2 ,x 3 ,,,x m Then select the basis vector with better condition number to replace the original basis vector, and re-divide the basis variables according to the elementary transformation matrix of LU decomposition into Substitute the new basis variables into the updated Jacobian matrix. The relationship between the Jacobian matrix and the basic matrix in the second variable partition form is The new Jacobian matrix has the following relationship with the previous Jacobian matrix: P is the elementary permutation matrix obtained by LU decomposition, which is used to perform simple row permutation operations on the original matrix. P must satisfy P T P=PP T =I. Re-decompose the basis matrix after the basis change

[0077] Step B4, verify the new condition number;

[0078] After the base change operation is completed, the condition number of the updated basic matrix is ​​recalculated and verified to see if it meets the preset suitability requirements. If the new condition number meets the requirements, the variables after the base change are replaced; otherwise, the basis variables are adjusted until a suitable basis is found.

[0079] Compared with the existing technology for solving the optimization of the double-tower simultaneous process, the present invention has the following significant technical advantages:

[0080] 1. Improve the stability of the solution: By judging the correlation of the constraint equations and optimizing the condition number, the present invention can effectively avoid the influence of ill-conditioned problems and significantly improve the numerical stability in the solution process.

[0081] 2. Speed ​​up the solution: Removing redundant constraint equations and improving the condition number by LU decomposition can reduce unnecessary calculations and improve the convergence speed of iteration. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1This is a double tower simulation diagram of the depropanizer and debutanizer;

[0083] Figure 2 It is a schematic diagram of the distillation tower model;

[0084] Figure 3 is the inflow and outflow model of the jth tower plate;

[0085] Figure 4 Solve the overall flow chart for the LU_SQP algorithm for double-tower optimization.

[0086] Figure 5 This is a flowchart for handling redundancy of constraint equations in twin-tower optimization; DETAILED DESCRIPTION

[0087] The present invention will be further described below in conjunction with the accompanying drawings.

[0088] like Figure 1 As shown, the double tower system consists of two distillation towers, namely a depropanizer (E-DA-404) and a debutanizer (E-DA-405), which are mainly used to separate different components in the raw materials.

[0089] The function of the depropanizer (E-DA-404) is to separate propane from other heavier components in the raw materials: in this process, propane is collected as a lighter component in liquid form at the top of the tower and becomes a product; the heavier components flow out from the bottom of the tower. The feed of the depropanizer (E-DA-404) is divided into two parts: one part is the bottom liquid from the condensate stripping tower and the bottom liquid from the C3 liquefied gas redistillation tower, which are introduced into the 34th tower plate of the depropanizer (stream S502). The other part is the bottom liquid from the deethanizer, which is introduced into the 18th tower plate of the depropanizer (stream S538).

[0090] The condensate stripping tower is a device used to separate liquid hydrocarbons and gases. The bottom liquid contains more heavy components, such as C4 and above hydrocarbons. The C3 liquefied gas redistillation tower is a device used to further separate and purify propane. The bottom liquid may contain heavier components that are not completely separated. C3 and C4 represent carbon three and carbon four, respectively. They are hydrocarbons composed of only carbon and hydrogen, also called hydrocarbons. Carbon three contains 3 carbon atoms and carbon four has 4 carbon atoms.

[0091] like Figure 2As shown, the top distillate of the depropanizer (E-DA-404) is condensed by propylene refrigerant, and the condensed liquid enters the reflux tank. Part of the condensed liquid is pumped out as the reflux L1 of the tower, and the other part (stream S511) is sent to the propylene dryer as the top liquid product S1 (mainly propane) of the depropanizer (E-DA-404). In order to reduce the heat load of the bottom reboiler and save low-pressure steam, an intermediate reboiler is provided in the middle of the depropanizer, using quenching water as the heat source. The bottom product LN of the depropanizer (E-DA-404) is a bottom liquid (stream S503) containing C4 and heavier components, and stream S503 flows into the feed port of the debutanizer (E-DA-405).

[0092] In the debutanizer (E-DA-405), the input is the bottom liquid containing C4 and heavier components output from the depropanizer (E-DA-404). The C4 component is then separated at the top of the tower (the top vapor phase product V1 of the debutanizer (E-DA-405)) and collected as a product, while the heavy components C5 and above flow out from the bottom of the tower (the bottom product LN of the debutanizer (E-DA-405)). C5 is a hydrocarbon compound containing 5 carbon atoms. The top distillate of the debutanizer (E-DA-405) is condensed with water and enters the reflux tank, where a portion of the liquid is pumped out as the tower reflux L1, and the other portion (stream S592) is sent out as a mixed C4 product (the top liquid phase product S1 of the debutanizer (E-DA-405)). The bottom reboiler uses low-pressure steam as a heat source. The bottom product of the debutanizer (E-DA-405) is light gasoline (stream S581), which is mixed with gasoline from the gasoline stripping tower, cooled, and sent to the pyrolysis gasoline hydrogenation unit for further processing.

[0093] A dual-tower system process optimization method based on LU decomposition specifically comprises the following steps:

[0094] Step 1, establish the equilibrium equations of the depropanizer and debutanizer, and normalize the equations;

[0095] The depropanizer and debutanizer described in the present invention are both distillation towers, and the steady-state process model of the distillation tower satisfies the material balance, energy balance, and vapor-liquid phase reaction balance equations, and the molar fraction composition of the vapor-liquid phase satisfies the normalization equation. The tray in the distillation tower is a device used to promote mass transfer and heat transfer between the gas and liquid phases; under steady-state operating conditions, the material balance, phase balance, and energy balance on the tray are key concepts for describing the behavior of the tray.

[0096] Assuming that the number of plates in the distillation tower is M, the condenser and the reboiler can be regarded as two upper and lower plates during the modeling process, then the actual number of plates is N=M+2.

[0097] Each tray has a feed, a vapor phase side discharge and a liquid phase side discharge. Under steady-state operation, the amount of material entering a tray is equal to the amount of material leaving the tray because there is no material accumulation in the tray. For each tray, the sum of the total molar flow of the liquid and vapor phases entering the tray is equal to the sum of the total molar flow of the liquid and vapor phases leaving the tray.

[0098] Since the total number of moles of the i-th component on the first tray is conserved during the reaction, the number of moles of the i-th component in the gas phase entering the first tray is V 1 y i,1 Add the number of moles of the i-th component in the liquid phase leaving through the feed and the first tray (S 1 +L 1 )x i,1 Equal to the number of moles of the i-th component in the gas phase leaving the first tray to the second tray V 2 y i,2 ; The material balance equation for the first plate is obtained as follows:

[0099] V 2 y i,2 -V 1 y i,1 -(S 1 +L 1 )x i,1 =0;

[0100] Among them, V 2 is the gas phase flow rate leaving the first tray, y i,2 is the mole fraction of component i in the gas phase leaving the first tray, V 1 is the gas phase flow rate entering the first tray, y i,1 is the mole fraction of the i-th component in the gas phase entering the first tray, S 1 is the liquid phase flow rate entering the first tray, L 1 is the liquid phase flow rate leaving the first tray, x i,1 is the mole fraction of the i-th component in the liquid phase leaving the first tray.

[0101] like Figure 3 As shown in the figure, since the total molar number of the i-th component on the j-th plate is conserved during the reaction, the molar number of the i-th component in the liquid phase and the gas phase entering the j-th plate (F j z i,j +L j-1 x i,j-1 +V j+1 y i,j+1 ) is equal to the number of moles of the i-th component in the liquid and gas phases leaving the tower plate j (L j +S j )x i,j and (V j+G j )y i,j ); The material balance equation for the jth tray is:

[0102] F j z i,j +L j-1 x i,j-1 +V j+1 y i,j+1 -(L j +S j )x i,j -(V j +G j )y i,j =0,(2≤j≤N-1);

[0103] Among them, F j is the liquid phase flow rate entering the jth tray, z i,j The mole fraction of the i-th component in the liquid phase entering the j-th tray, L j-1 Liquid flow rate leaving from the j-1th tray, x i,j-1 The mole fraction of the i-th component in the liquid phase leaving the j-1-th plate, V j+1 The gas phase flow rate entering the jth tower plate. i,j+1 The mole fraction of the i-th component in the gas phase entering the j+1-th tray, L j The liquid phase flow rate leaving the jth tray, S j The side liquid phase flow rate leaving from the jth tray, G j Side gas flow rate leaving from the jth tray, y i,j is the mole fraction of the i-th component in the gas phase entering the j-th tower plate.

[0104] Since the total molar number of the i-th component on the N-th plate is conserved during the reaction, the molar number of the i-th component in the liquid phase leaving the N-1-th plate is L N-1 x i,N-1 Subtract the moles of the i-th component in the gas and liquid phases leaving the N-th plate (V N +G N )y i,N and L N x i,N is equal to zero; the material balance equation for the Nth plate is obtained as follows:

[0105] L N-1 x i,N-1 -(V N +G N )y i,N -L N x i,N =0;

[0106] Among them, L N-1is the liquid phase flow rate leaving from the N-1th tray, x i,N-1 is the mole fraction of the i-th component in the liquid phase leaving the N-1th plate, V N is the gas phase flow rate leaving from the Nth plate, G N is the side gas flow rate leaving from the Nth tray, y i,N is the mole fraction of the i-th component in the gas phase leaving the N-th plate. N is the liquid phase flow rate leaving from the Nth plate, x i,N is the mole fraction of component i in the liquid phase leaving the Nth tray.

[0107] For the first tray condenser, the total enthalpy of the gas phase entering the condenser plus the added heat equals the total enthalpy of the gas phase leaving the condenser plus the total enthalpy of the liquid phase leaving the condenser. The energy balance equation for the first tray is:

[0108]

[0109] Among them, V 2 represents the gas phase flow rate leaving from the second tray, represents the specific enthalpy of the gas phase leaving the second tray, V 1 represents the gas phase flow rate entering the first tray 1, represents the specific enthalpy of the gas phase entering the first tray, L 1 represents the liquid phase flow rate leaving from the first tray, S 1 represents the liquid phase flow rate entering the first tray, represents the specific enthalpy of the liquid phase leaving the first tray, Q 1 Represents the heat added or removed on the first stage.

[0110] According to the law of conservation of energy, the total enthalpy of the gas and liquid phases entering the jth plate is Plus the added heat Q j and the enthalpy of the feed = The total volume of gas and liquid leaving the jth plate and The energy balance equation for the jth plate is obtained:

[0111]

[0112] Among them, F j represents the liquid phase flow rate entering the jth tray, represents the specific enthalpy of the liquid phase entering the jth tray, V j+1 represents the gas phase flow rate entering the j+1th tray, represents the specific enthalpy of the gas phase entering the j+1th tray, L j-1 represents the liquid phase flow rate leaving from the j-1th tray, represents the specific enthalpy of the liquid phase leaving the j-1th plate, Q j represents the heat added or removed at stage j, V j represents the gas phase flow rate leaving from plate j, G j represents the side gas flow rate leaving tower j, S j Side stream liquid phase leaving from column j. represents the specific enthalpy of the gas phase leaving the jth plate, represents the specific enthalpy of the liquid phase leaving from the j-th tray.

[0113] The total enthalpy of the liquid phase entering the reboiler on the last tray Plus the added heat Q N Equal to the total enthalpy of the vapor and liquid phases leaving the reboiler The energy balance equation for the Nth plate is obtained:

[0114]

[0115] Among them, L N-1 represents the liquid phase flow rate leaving from the N-1th tray, is the specific enthalpy of the liquid phase leaving the N-1th plate, Q N is the heat added or removed at the Nth stage, V N is the gas phase flow rate leaving from the Nth plate, G N is the side gas flow rate leaving from the Nth tray, is the specific enthalpy of the gas phase leaving the Nth tray.

[0116] Vapor-liquid phase reaction equilibrium equation: i,j =K i,j x i,j (1≤i≤C,1≤j≤N);

[0117] Among them, y i,j is the mole fraction of the i-th component in the gas phase of the j-th tray, x i,j The mole fraction of the i-th component in the liquid phase on the j-th tray; K i,j is the equilibrium constant, which represents the relative concentration relationship between the gas phase and the liquid phase of the i-th component on the tower plate j. It is usually a constant related to temperature and pressure. C is the total number of components in the system, and N is the total number of tower plates.

[0118] In any given mixture, the sum of the mole fractions of all components must equal 1. That is, the total amount of all components in the mixture must account for 100% of the entire mixture. In the context of gas and liquid phases, the normalization equation can be applied to the gas phase and the liquid phase separately. The molecular normalization equation is:

[0119] Among them, y i,jis the mole fraction of the i-th component in the gas phase of the j-th tray, x i,j is the mole fraction of the i-th component in the liquid phase of the j-th tower plate, and C is the total number of components in the system.

[0120] Table 1 shows the model variables and their correlation coefficients of the middle distillation tower, and the symbols, names and units of all variables in this embodiment.

[0121] Table 1

[0122]

[0123] Step 2: Establish an optimization model;

[0124] Enterprises always aim to pursue maximum economic benefits, which can usually be reflected in the quantity, quality and price of products produced. In the ethylene production process, the depropanizer and debutanizer produce C3 and C4 liquid products respectively, which are the benefits that enterprises can directly see. The C3 product extracted from the first tray of the depropanizer and the C4 product extracted from the first tray of the debutanizer are the main sources of value. In contrast, the output from the bottom of the debutanizer is usually of lower value and is generally used as raw materials or fuel for other production processes, so it is not the main consideration here.

[0125] In order to optimize the production process, the total value of the liquid products produced by the depropanizer and debutanizer is taken as the optimization target, while taking into account the impact of product prices.

[0126] Based on the market price of C3 and C4 products and the output of the depropanizer and debutanizer, the total value of the products can be calculated. At present, the market price of propylene is about $610 / ton, while the price of butane is about $620 / ton. As a basis for optimizing the production process, the following optimization objective function is established:

[0127] Max:F=(610*S1_bw+620*S1_dw) / 10000;

[0128] Wherein, F represents the optimization objective function, and its unit is 10,000 yuan / hour; S1_bw and S1_dw represent the liquid product flow rates at the top of the depropanizer and debutanizer, respectively, and their units are Kmol / hr.

[0129] In the distillation process of the distillation tower, increasing the output of the product will affect the quality composition of the product. In order to ensure that the product meets the quality requirements, it is necessary to limit the key components of the product. In the depropanizer of this article, the main component of the liquid product is C3 product, with a content of more than 99%, among which C3H6 is the largest, which is 96.7% at the working point. C3H6 in the product can be used as the main indicator of the purity of C3 product. For the debutanizer, the liquid product is mainly C4 component, and C5 and components above C5 are slag in the product. For this reason, the content of C5 and components above C5 must be controlled. Here, the content of C5 is set as an indicator to examine the purity of the debutanizer C4 product.

[0130] The optimization model is established by combining the objective function F, constraint equations (three major equilibrium equations and normalization equation of the distillation tower) and boundary constraints.

[0131] In this embodiment, the boundary constraints are:

[0132] Depropanizer top product: C3H6≥0.93;

[0133] Debutanizer top product: C5≤0.0095;

[0134] It means that the proportion of C3H6 in the top product of the depropanizer tower is not less than 0.93, and the proportion of C5 in the top product of the debutanizer tower is not greater than 0.0095. For this optimization objective function and constraint equation, the dimension of all variables involved in the double tower model in this embodiment is 4190, and the dimension of fixed variables is 103.

[0135] Step 3: Solve the optimization model established in step 2;

[0136] Based on the large scale and high sparsity of the matrix of the twin tower problem, the quadratic sequence programming (SQP) method is selected to solve it. SQP is particularly suitable for solving large-scale and sparse nonlinear programming problems. The core idea of ​​the method is to approximate the original nonlinear problem by constructing a quadratic programming subproblem in each iteration, and then solve this subproblem to obtain the search direction and step size. It transforms the original complex problem into a series of easier-to-handle quadratic programming problems. Combined with LU decomposition to process the constraint matrix, solve the correlation, and perform infeasible recovery, the solution efficiency and accuracy are improved.

[0137] like Figure 4 As shown, the specific solution steps of the LU_SQP algorithm are as follows:

[0138] Step 3.1: Preprocess the constraint equations:

[0139] The purpose of preprocessing is to convert the original model of the problem into a standard form suitable for numerical solver processing, ensuring that the subsequent algorithm can be effectively calculated. The double tower optimization problem is preprocessed, including determining the feasible initial values ​​of the optimization variables and scaling the variables with large unit range differences. That is, by adjusting the variables to similar scales, it is avoided that certain variables have too much influence on the results during the optimization process. Since the general variables of the distillation tower have their clear physical meanings, they can be scaled according to the physical meaning and size of the variables. The scaled variables optimized for the distillation tower in this embodiment are shown in Table 2:

[0140] Table 2

[0141]

[0142]

[0143] Step 3.2, derivative of constraint equation:

[0144] Each plate in the double tower system has a corresponding equilibrium equation, which combines the material balance equation, energy balance equation, vapor-liquid phase reaction balance equation and the normalized equation of the vapor-liquid phase molar fraction group involved, with a total of 4293 variables and 3390 constraints. The equation is denoted as c i (x), i = 1....3390, the variable is denoted as x i , i=1...4293. Each constraint calculates partial derivatives for each variable. Since there are too many equations and variables, they are not listed here. The calculated derivative matrix is ​​called the Jacobian matrix, denoted by ceqGRAD(x). ceqGRAD(x) has the following form:

[0145]

[0146] Among them, x represents all the variables involved in the twin tower optimization process: x 1 ,x 2 ,...,x n , c 1 ,c 2 ,...,c m is the set of constraint equations, is the partial derivative of the ith equilibrium equation with respect to the jth variable.

[0147] Step 3.3: Remove redundant constraints by determining whether the constraint equations are linearly related and consistent:

[0148] In each iteration, the rank of the Jacobian matrix of the constraint is indirectly determined through LU decomposition. Because the twin-tower optimization involves nonlinear equations, it is difficult to directly extract a clear constraint coefficient matrix like linear programming. Therefore, the correlation of the constraints is determined by judging the rank of the derivative matrix of the constraint equation, that is, the Jacobian matrix. After performing LU decomposition on the Jacobian matrix, we get the L matrix and the U matrix, L*U=lu(ceqGRAD), and then according to the structure of the U matrix, we can quickly determine the rank of the matrix and thus determine the independence of the constraint equation. If the constraint equations are linearly related and consistent, these redundant constraint equations can be removed to reduce the occupancy of system resources during the calculation process and accelerate convergence. If they are related but inconsistent, the system will immediately determine that the problem has no solution, thereby avoiding invalid calculations in subsequent iterations. Figure 5 As shown, the specific steps of LU decomposition to remove redundant constraints are:

[0149] Step A1, setting the mark error topDep to 1e-10;

[0150] Computers use floating-point numbers to represent real numbers, but the representation precision of floating-point numbers is limited and cannot accurately represent all real numbers. In particular, after a series of operations, a number that should theoretically be zero may become a small value that is very close to zero but not zero. This error comes from rounding or truncation of limited precision. For example, in some calculations, the theoretical result may be 0, but the actual result may be 1.23*e-10. If you directly judge whether this result is equal to zero, you will get an incorrect judgment. Therefore, it is necessary to judge whether the result is close to zero in a very small range, rather than equal to zero.

[0151] Step A2: LU decomposition;

[0152] Convert the Jacobian matrix ceqGRAD obtained in step 3.2 into a sparse form and perform LU decomposition. ceqGRAD contains many zero elements. Storing this matrix in the form of a sparse matrix can avoid storing and calculating a large number of zero elements. Sparse LU decomposition decomposes the sparse matrix into a lower triangular matrix L and an upper triangular matrix U, as well as the related permutation matrix P. The calculation will skip the zero elements in the sparse matrix, which greatly saves computing resources and time.

[0153] Step A3: determine the relevant redundancy of constraints;

[0154] Whether the constraints are relevant is determined by judging whether the U matrix obtained by LU decomposition of the Jacobian matrix of the constraints has all-zero rows at the current iteration point. If a row of the U matrix is ​​all-zero at the current iteration point, it means that the constraints of the corresponding row are relevant, that is, the redundancy needs to be removed. Because computers use floating-point numbers to represent real numbers, the judgment of all zeros here is actually compared with a very small number tolDep. The results show that the 67th and 2301st rows of the U matrix have all-zero rows, which means that the 67th and 2301st constraints are redundant, and the consistency is determined in the next step. If there is no all-zero row, it means that there are no redundant constraints, so continue with the subsequent solution and skip steps A4 and A5.

[0155] Step A4: Determine the consistency of the constraints;

[0156] Determine whether the product of the L matrix and the constraint c matrix is ​​zero in the row corresponding to the redundant constraint. If the row corresponding to the redundant constraint is zero, it means that the redundant constraints are consistent, and the redundancy can be removed by directly deleting this row of constraints. If the product of the L matrix and the constraint c matrix in the row corresponding to the redundant constraint is not zero, it means that the redundant constraints are inconsistent, the constraint equations are contradictory, and can no longer be solved. The optimization process can be terminated and a no solution prompt is output. In this embodiment, after calculation, the corresponding rows of rows 67 and 2301 are equal to 0, indicating that the two constraints of rows 67 and 2301 of the combined tower are related and consistent, and one of the rows of constraints can be eliminated through step A6 to continue solving.

[0157] Step A5: remove redundant constraints and update the constraint system;

[0158] The elementary permutation matrix P obtained through LU decomposition finds the index 67 and 2301 of the redundant constraints, records the 2301th row as the constraint index to be removed, marks it as remove, and sets the remove row to be empty in the constraint matrix and Jacobian matrix, updates the constraint system, and now one constraint is reduced, and the subsequent optimization solution continues. Because the redundant equations have been removed, only independent and consistent equations are retained, and the updated constraint equation set will be simpler, thereby reducing the computational burden and speeding up the solution. At the same time, after removing the redundant constraints, the independence of the constraints is enhanced, making the Jacobian matrix more favorable, reducing the possibility of ill-conditioned matrices, and improving numerical stability, laying a good foundation for the next step of optimizing the condition number of the matrix.

[0159] Step 3.4, iterative solution, change basis to deal with singular matrix;

[0160] The SQP algorithm is called to iteratively solve the model, converting the complex nonlinear optimization problem into a series of quadratic programming sub-problems. In each iteration, the algorithm constructs a local quadratic programming sub-problem based on the current iteration point x. k The quadratic approximation of the objective function F and the constraint c(x) can approximate the local characteristics of the original problem.k The QP subproblem is expressed as follows:

[0161]

[0162] Where, d k is the search direction vector, g(x k ) is the derivative of the optimization objective function F, J(x k ) T is the constraint c(x k ), the derivative Jacobian matrix of W(x k ) is the second derivative matrix of the Lagrangian function, which is in g(x k ) and J(x k ) T The Lagrangian function is a function that integrates the constraint and the objective function, and its form is: L(x,λ,v,π)=F(x)+λ T c i (x)+v T (xx U )-π T (xx L ), λ, v, π are the corresponding coefficients. L represents the lower bound of the variable x, x U Represents the upper bound of the variable x.

[0163] The solution to the quadratic programming problem is relatively easy. We can directly call the existing algorithm to solve the equation and obtain the search direction d k , used to update the current iteration point x k , the update formula is: k+1 =x k +d k After updating the iteration point, check whether the optimization termination conditions are met, such as the improvement of the objective function, the excessive number of iterations, or the violation of the equilibrium equation constraints. k The solution depends on the Jacobian matrix of the constraints, so once the matrix becomes singular, it is impossible to continue solving.

[0164] For the above singular cases, the constrained Jacobian matrix J is decomposed by LU:

[0165]

[0166] Among them J B The Jacobian matrix corresponding to the basic variables is called the basic matrix, J N Represents the Jacobian matrix of the dependent variable. The corresponding basic variable is: x 1 ,x 2 ,x 3 ,,,xm After LU decomposition, some basic variables and dependent variables are replaced by permutation matrix to replace the basic variables that cause the matrix singularity and restore the solvability of the matrix. The re-partitioned basic variables are The specific process of processing singular matrix basis change is:

[0167] Step B1, initializing the condition number;

[0168] The condition number can be regarded as a measure of the degree to which a matrix is ​​close to being singular. When the condition number is large, it means that the determinant of the matrix tends to zero and the matrix tends to be singular. Let the condition number κ(J B ) is initialized to 1. Then in each iteration, the current iteration point x k The basic matrix J B Calculate its condition number. The calculation formula of the condition number is: κ(J B )=||J B ||·||J B -1 ||, where ||J B || indicates J B 2-norm of the matrix.

[0169] Step B2, determine whether the condition number is acceptable;

[0170] The judgment condition for changing the basis is: the calculated condition number κ(J B ) and the pre-set threshold κ threshlod For comparison, if the current matrix J B The condition number κ'(J B ) is greater than the threshold and relative to the previous condition number κ(J B ) has increased by 20 times. In this embodiment, the threshold κ is taken threshlod is 1000, and in the first 7 iterations, κ(J B ) is less than the threshold, indicating that the matrix is ​​well-posed. In the eighth calculation, the current matrix J B The condition number reaches 1.2*e13, which is 20 times greater than the previous condition number, which means that the basic matrix is ​​singular and the calculation cannot be performed. A basis change operation is required to recover.

[0171] Step B3, LU decomposition for base change;

[0172] When the condition number is poor, the basis change operation is performed through LU decomposition, and the LU decomposition is performed on the constrained Jacobian matrix J. The relationship between the basic matrix and the Jacobian matrix is: Among them J B represents the basic matrix, J N Represents the Jacobian matrix of the dependent variable. The corresponding basis variable is: x 1 ,x2 ,x 3 ,,,x m Then select the basis vector with better condition number to replace the original basis vector, and re-divide the basis variables according to the elementary transformation matrix of LU decomposition into Substitute the new basis variables into the updated Jacobian matrix. The relationship between the Jacobian matrix and the basic matrix in the second variable partition form is The new Jacobian matrix has the following relationship with the previous Jacobian matrix: P is the elementary permutation matrix obtained by LU decomposition, which is used to perform simple row permutation operations on the original matrix. P must satisfy P T P=PP T =I. Re-decompose the basis matrix after the basis change

[0173] Step B4, verify the new condition number;

[0174] After the base change operation is completed, the condition number of the updated basic matrix is ​​recalculated and verified to see whether it meets the preset suitability requirements. If the new condition number meets the requirements, the variables after the base change are replaced; otherwise, the base variables are adjusted until a suitable base is found. In the embodiment, the condition number judgment condition is met after the base change once, and the subsequent solution process is entered.

[0175] Step 3.5: Determine whether the optimization conditions are met:

[0176] Determine whether the optimization conditions are met. If not, introduce a new x k+1 , repeat steps 3.2-3.5 to continue the next iteration solution; if satisfied, terminate the calculation and output the optimization results: the maximum benefit value of the dual-tower system and the values ​​of each operating variable when the maximum benefit is obtained.

[0177] While ensuring that the operating conditions are feasible, the method implemented in this patent is used to optimize the solution of the dual-tower system. The software used is Matlab, version R2022b, and the operating environment is a 64-bit Windows 11 laptop. The solution results are compared with the existing solver SNOPT. The calculation results show that under the same given conditions such as initial values ​​and convergence accuracy, both the patented method and SNOPT can correctly converge to the optimal value. However, the solution efficiency varies greatly, and the calculation efficiency of the patented method is much higher than that of SNOPT. Table 3 shows the calculation conditions and optimization results of the two methods.

[0178] Table 3 The solution results of the optimization problem using the patented method and SNOPT

[0179]

[0180]

[0181] As can be seen from Table 3, the calculation time of SNOPT is much longer than that of the method of this patent, and the number of iterations of the method of this patent is also much smaller than that of SNOPT. The main reason for this result is that the convergence speed of SNOPT in solving the optimization problem is too slow. It may be that the model is not ideal and the constraints are redundant, which causes SNOPT to add a lot of unnecessary calculations. In addition, from the calculation results of the objective function value, it can be seen that the objective function solved by the method of this patent has higher accuracy.

[0182] In addition to being used in the optimization of a double-tower distillation tower, the method of the present invention can also be extended to the optimization problem of multiple towers, laying a good foundation for the simulation and optimization of the entire large-scale ethylene production process. In the future, it can be applied to various chemical process optimization problems, especially in application scenarios involving complex constraints and pathological problems. The effect is particularly obvious, providing an efficient and stable optimization solution for related industries.

Claims

1. A dual-tower system process optimization method based on LU decomposition, specifically comprising the following steps: Step 1, establish the equilibrium equations of the depropanizer and debutanizer, and normalize the equations; The depropanizer and debutanizer are both distillation towers, and the steady-state process model of the distillation tower satisfies the material balance, energy balance, and vapor-liquid phase reaction balance equations, and the molar fraction composition of the vapor-liquid phase satisfies the normalization equation; Assuming that the number of plates in the distillation tower is M, in the modeling process, the condenser and the reboiler can be regarded as two upper and lower plates, then the actual number of plates is N = M + 2; there is a feed, a vapor phase side discharge and a liquid phase side discharge on each plate. Under steady-state operation, because there is no material accumulation in the plate, the amount of material entering a plate is equal to the amount of material leaving the plate; for each plate, the sum of the total molar flow of the liquid phase and the gas phase entering the plate is equal to the sum of the total molar flow of the liquid phase and the gas phase leaving the plate; The material balance equation for the first plate is obtained as follows: V2y i,2 -V1y i,1 -(S1+L1)x i,1 =0; in, V2 is the gas phase flow rate leaving the first tray, y i,2 is the mole fraction of the i component in the gas phase leaving the first tray, V1 is the gas phase flow rate entering the first tray, y i,1 is the mole fraction of the i-th component in the gas phase entering the first tray, S1 is the liquid phase flow rate entering the first tray, L1 is the liquid phase flow rate leaving the first tray, x i,1 is the mole fraction of the i-th component in the liquid phase leaving the first tray; The material balance equation for the jth plate is obtained as follows: F j z i,j +L j-1 x i,j-1 +V j+1 y i,j+1 -(L j +S j )x i,j -(V j +G j )y i,j =0,(2≤j≤N-1); Among them, F j is the liquid phase flow rate entering the jth tray, z i,j The mole fraction of the i-th component in the liquid phase entering the j-th tray, L j-1 Liquid flow rate leaving from the j-1th tray, x i,j-1 The mole fraction of the i-th component in the liquid phase leaving the j-1-th plate, V j+1 The gas phase flow rate entering the jth tray; y i,j+1 The mole fraction of the i-th component in the gas phase entering the j+1-th tray, L j The liquid phase flow rate leaving the jth tray, S j The side liquid phase flow rate leaving from the jth tray, G j The side gas flow rate leaving the jth tray, y i,j is the mole fraction of the i-th component in the gas phase entering the j-th tray; The material balance equation for the Nth plate is obtained as follows: L N-1 x i,N-1 -(V N +G N )y i,N -L N x i,N =0; Among them, L N-1 is the liquid phase flow rate leaving from the N-1th tray, x i,N-1 is the mole fraction of the i-th component in the liquid phase leaving the N-1th plate, V N is the gas phase flow rate leaving from the Nth plate, G N is the side gas flow rate leaving from the Nth tray, y i,N is the mole fraction of the i-th component in the gas phase leaving the N-th plate; L N is the liquid phase flow rate leaving from the Nth plate, x i,N is the mole fraction of component i in the liquid phase leaving the Nth tray; The energy balance equation for the first tray is obtained as follows: Wherein, V2 represents the gas phase flow rate leaving from the second tray, represents the specific enthalpy of the gas phase leaving the second tray, V1 represents the gas phase flow rate entering the first tray 1, represents the specific enthalpy of the gas phase entering the first tray, L1 represents the liquid phase flow rate leaving the first tray, S1 represents the liquid phase flow rate entering the first tray, represents the specific enthalpy of the liquid phase leaving the first tray, Q1 represents the heat added or removed at the first tray; The energy balance equation for the jth plate is obtained as follows: Among them, F j represents the liquid phase flow rate entering the jth tray, represents the specific enthalpy of the liquid phase entering the jth tray, V j+1 represents the gas phase flow rate entering the j+1th tray, represents the specific enthalpy of the gas phase entering the j+1th tray, L j-1 represents the liquid phase flow rate leaving from the j-1th tray, represents the specific enthalpy of the liquid phase leaving the j-1th plate, Q j represents the heat added or removed at stage j, V j represents the gas phase flow rate leaving from plate j, G j represents the side gas flow rate leaving tower j, S j The side stream liquid phase flow rate leaving from column j; represents the specific enthalpy of the gas phase leaving the jth plate, represents the specific enthalpy of the liquid phase leaving the jth tray; The energy balance equation for the Nth plate is obtained as follows: Among them, L N-1 represents the liquid phase flow rate leaving from the N-1th tray, is the specific enthalpy of the liquid phase leaving the N-1th plate, Q N is the heat added or removed at the Nth stage, V N is the gas phase flow rate leaving from the Nth plate, G N is the side gas flow rate leaving from the Nth tray, is the specific enthalpy of the gas phase leaving the Nth tray; The vapor-liquid phase reaction equilibrium equation is obtained as follows: i,j =K i,j x i,j (1≤i≤C,1≤j≤N); Among them, y i,j is the mole fraction of the i-th component in the gas phase of the j-th tray, x i,j The mole fraction of the i-th component in the liquid phase on the j-th tray; K i,j is the equilibrium constant, which indicates the relative concentration relationship between the gas phase and the liquid phase of the i-th component on the tower plate j, and is usually a constant related to temperature and pressure; C is the total number of components in the system, and N is the total number of tower plates; The molecular normalization equation is: Among them, y i,j is the mole fraction of the i-th component in the gas phase of the j-th tray, x i,j is the mole fraction of the i-th component in the liquid phase of the j-th tray, and C is the total number of components in the system; Step 2: Establish an optimization model; The total value of the liquid products produced by the depropanizer and debutanizer is taken as the optimization target. Based on the market prices of C3 and C4 products and the outputs of the depropanizer and debutanizer, the total value of the products is calculated. This is used as a basis for optimizing the production process to establish the optimization objective function F. The objective function F, constraint equations, and boundary constraints are combined to establish an optimization model. The constraint equations are the three equilibrium equations and the normalization equation in step one. Step 3: Solve the optimization model established in step 2; The quadratic sequence programming method SQP is chosen for solution. In each iteration, a quadratic programming subproblem is constructed to approximate the original nonlinear problem, and then the subproblem is solved to obtain the search direction and step size; the constraint matrix is ​​processed by combining LU decomposition, the correlation is solved, and infeasible recovery is performed; the specific solution steps are as follows: Step 3.1: Preprocess the constraint equations: Convert the original model of the problem into a standard form suitable for processing by the numerical solver, including determining the initial values ​​of the optimization variables and scaling the variables with large unit range differences; Step 3.2, derivative of constraint equation: Each plate in the double tower system has a corresponding equilibrium equation, which combines the material balance equation, energy balance equation, vapor-liquid phase reaction balance equation and the normalized equation of the vapor-liquid phase molar fraction group involved; each constraint calculates the partial derivative of each variable separately, and the calculated derivative matrix is ​​called the Jacobian matrix ceqGRAD(x): Among them, x represents all variables involved in the twin tower optimization process: x1, x2, ..., x n ,c1,c2,...,c m is the set of constraint equations, is the partial derivative of the ith equilibrium equation with respect to the jth variable; Step 3.3: Remove redundant constraints by determining whether the constraint equations are linearly related and consistent: In each iteration, the rank of the Jacobian matrix of the constraint is indirectly determined by LU decomposition, and the correlation of the constraint is determined by determining the rank of the derivative matrix of the constraint equation, that is, the Jacobian matrix. After LU decomposition of the Jacobian matrix, the L matrix and the U matrix are obtained, L*U=lu(ceqGRAD), and then the rank of the matrix is ​​quickly determined according to the structure of the U matrix, thereby determining the independence of the constraint equation. If the constraint equations are linearly related and consistent, these redundant constraint equations can be removed; if they are related but inconsistent, the system will immediately determine that the problem has no solution. Step 3.4, iterative solution, change basis to deal with singular matrix; The SQP algorithm is called to iteratively solve the model, converting the complex nonlinear optimization problem into a series of quadratic programming sub-problems; in each iteration, the algorithm constructs a local quadratic programming sub-problem based on the current iteration point x k The objective function F and the constraint c(x) are quadratically approximated; at the iteration point x k The QP subproblem is expressed as follows: Where, d k is the search direction vector, g(x k ) is the derivative of the optimization objective function F, J(x k ) T is the constraint c(x k ), the derivative Jacobian matrix of W(x k ) is the second derivative matrix of the Lagrangian function, which is in g(x k ) and J(x k ) T On the basis of , we can continue to find the first derivative to get; the Lagrangian function is a function that integrates the constraint and the objective function, and the form is: L(x,λ,v,π)=F(x)+λ T c i (x)+v T (xx U )-π T (xx L ), λ, v, π are the corresponding coefficients; x L represents the lower bound of the variable x, x U Represents the upper bound of the variable x; Call the existing algorithm to solve the equation and get the search direction d k , used to update the current iteration point x k , the update formula is: k+1 =x k +d k After updating the iteration point, check whether the optimization termination conditions are met, such as the degree of improvement of the objective function, the excessive number of iterations, or the degree of violation of the equilibrium equation constraints; because d k The solution depends on the Jacobian matrix of the constraints, so once the matrix becomes singular, it is impossible to continue solving; For the above singular cases, the constrained Jacobian matrix J is decomposed by LU: Among them J B The Jacobian matrix corresponding to the basic variables is called the basic matrix, J N Represents the Jacobian matrix of dependent variables. The corresponding basic variables are: x1,x2,x3,,,x m ; After LU decomposition, some basic variables and dependent variables are replaced by permutation matrix to replace the basic variables that cause the matrix singularity and restore the solvability of the matrix. The re-divided basic variables are Step 3.5: Determine whether the optimization conditions are met: Determine whether the optimization conditions are met. If not, introduce a new x k+1 , repeat steps 3.2-3.5 to continue the next iteration solution; if satisfied, terminate the calculation and output the optimization results: the maximum benefit value of the dual-tower system and the values ​​of each operating variable when the maximum benefit is obtained.

2. The dual-tower system process optimization method based on LU decomposition as claimed in claim 1, characterized in that: The specific steps for removing redundant constraints described in step 3.3 are: Step A1, setting the mark error topDep to 1e-10; Step A2, converting the Jacobian matrix ceqGRAD obtained in step 3.2 into a sparse form and performing LU decomposition; ceqGRAD contains many zero elements; sparse LU decomposition decomposes the sparse matrix into a lower triangular matrix L and an upper triangular matrix U, and a related permutation matrix P; Step A3: determine the relevant redundancy of constraints; Whether the constraints are relevant is determined by judging whether the U matrix obtained by LU decomposition of the Jacobian matrix of the constraints has all zero rows at the current iteration point; If the U matrix has a row of all zeros at the current iteration point, it means that the constraints of the corresponding row are related, that is, the redundancy needs to be removed, and the consistency is judged in the next step; if there is no all-zero row, it means that there is no redundant constraint, so skip steps A4 and A5 and continue with the subsequent solution; Step A4, determine the consistency of the constraints; determine whether the product of the L matrix and the constraint c matrix is ​​zero in the row corresponding to the redundant constraint; if the row corresponding to the redundant constraint is zero, it means that the redundant constraints are consistent, and the row of constraints is directly deleted to remove the redundancy, and the solution is continued; if the product of the L matrix and the constraint c matrix in the row corresponding to the redundant constraint is not zero, it means that the redundant constraints are inconsistent, the constraint equations are contradictory, and the solution cannot be found, the optimization process is terminated, and a prompt of no solution is output; Step A5, remove redundant constraints and update the constraint system; find the index of the redundant constraints through the elementary permutation matrix P obtained by LU decomposition, empty the rows where the redundant constraints are located in the constraint matrix and Jacobian matrix, update the constraint system, remove redundant equations, retain independent and consistent equations, and continue the subsequent optimization solution.

3. The dual-tower system process optimization method based on LU decomposition as claimed in claim 1, characterized in that: The specific process of processing the singular matrix basis change described in step 3.4 is: Step B1, initializing the condition number; The condition number can be seen as a measure of how close a matrix is ​​to being singular; when the condition number is large, it means that the determinant of the matrix tends to zero and the matrix tends to be singular; Let the condition number κ of the Jacobian matrix (J B ) is initialized to 1; Then in each iteration, the current iteration point x k The basic matrix J B Calculate its condition number; The calculation formula of condition number is: κ(J B )=||J B ||·||J B -1 ||, where ||J B || indicates J B 2-norm of the matrix; Step B2, determine whether the condition number is acceptable; The judgment condition for changing the basis is: the calculated condition number κ(J B ) and the pre-set threshold κ threshlod For comparison, if the current matrix J B The condition number κ'(J B ) is greater than the threshold and relative to the previous condition number κ(J B ) increases by a set multiple; Step B3, LU decomposition for base change; When the condition number is poor, the basis change operation is performed through LU decomposition, and the LU decomposition is performed on the constrained Jacobian matrix J. The relationship between the basic matrix and the Jacobian matrix is: Among them J B represents the basic matrix, J N Represents the Jacobian matrix of the dependent variables. The corresponding basis variables are: x1, x2, x3,,,x m ; Then select the basis vector with better condition number to replace the original basis vector, and re-divide the basis variables according to the elementary transformation matrix of LU decomposition into Substitute the new basis variables into the updated Jacobian matrix. The relationship between the Jacobian matrix and the basic matrix in the second variable partition form is The new Jacobian matrix has the following relationship with the previous Jacobian matrix: P is the elementary permutation matrix obtained by LU decomposition, which is used to perform simple row permutation operations on the original matrix; P must satisfy P T P=PP T =I; Re-decompose the basis matrix after the basis change Step B4, verify the new condition number; After the basis change operation is completed, the condition number of the updated basic matrix is ​​recalculated and verified to see whether it meets the preset suitability requirements; if the new condition number meets the requirements, the variables after the basis change are replaced; otherwise, the basis variables are continued to be adjusted until a suitable basis is found.

Images