CN111914382B - Constraint evolution optimization method of atmospheric and vacuum device based on agent model - Google Patents
Constraint evolution optimization method of atmospheric and vacuum device based on agent model Download PDFInfo
- Publication number
- CN111914382B CN111914382B CN201910408795.2A CN201910408795A CN111914382B CN 111914382 B CN111914382 B CN 111914382B CN 201910408795 A CN201910408795 A CN 201910408795A CN 111914382 B CN111914382 B CN 111914382B
- Authority
- CN
- China
- Prior art keywords
- vector
- matrix
- output
- vectors
- particle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/30—Computing systems specially adapted for manufacturing
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Feedback Control In General (AREA)
Abstract
The invention discloses a constraint evolution optimization method of an atmospheric and vacuum device based on a proxy model, which aims at solving the problems of time-consuming calculation in the atmospheric and vacuum rectification process, and equation and inequality constraint optimization for guaranteeing stable working conditions and safety requirements, and aims at optimizing aiming at the maximum profit of production. Specifically, the method adopts a point adding strategy in optimization iteration, updates and perfects global and local proxy models, and performs optimization search on the built global and local proxy models so as to obtain an optimization solution. Compared with the traditional method, the method has the greatest advantage that the aimed optimization problem is more in line with the actual running condition, namely the constraint problems of equality and inequality related to the stability and safety requirements of the working condition are considered. In addition, the method of the invention searches global and local training data through newly defined optimization objective functions in a reciprocating and optimized way, thereby better ensuring the diversity of the training data and further continuously improving the accuracy of the proxy model.
Description
Technical Field
The invention relates to an optimization method for maximum profit of atmospheric and vacuum device production in the petrochemical field, in particular to a constraint evolution optimization method for atmospheric and vacuum device based on a proxy model.
Background
With the development of global economy, it is becoming increasingly important that the chemical production process be optimized in real time. Real-time optimization of the process is often expressed as achieving the aims of reducing energy consumption, reducing industrial emission or further improving economic benefit and product quality under the condition of ensuring product quality and operation condition safety in a short time as much as possible. As a faucet device for petroleum processing, atmospheric and vacuum rectification has a role in the oil refining industry, bears the primary separation of crude oil, and provides the weight of raw materials for subsequent oil refining chemical devices. The rectification process is a constraint, multivariable, nonlinear and strongly coupled complex system, and structurally, the rectification device can generally consist of dozens to hundreds of trays, complex mass and heat transfer operations are completed on the trays, and according to different separation requirements, different column scales can reach hundreds or even thousands of model equations, and the variables are complicated.
The complexity of the process model presents challenges to existing optimization methods. The optimization method of the rectification process generally comprises traditional optimization, commercial platform optimization, evolution optimization and the like. The traditional optimization method mainly comprises a sequence linear programming, a penalty function, an augmented Lagrangian function, a generalized reduced gradient, an interior point method, a sequence quadratic programming method (Sequential Quadartic Programming, SQP) and the like, wherein the SQP method is the main flow of the traditional optimization because of fast convergence and strong nonlinear processing, but the method has higher requirements on a mathematical model of an optimization proposition, and the iteration process can increase the difficulty for solving based on derivative information of the optimization proposition. In recent years, various foreign suppliers and high-tech companies have also introduced respective commercial optimization software, such as the RT-Opt of Aspen tech company, the ROMeo of Simulation Science company, the Profix Max of Honeywell company, the RTO+ of Emerson company, etc., however, the core optimization algorithm of these software is basically SQP, and the core technology is blocked strictly, and the optimization method cannot meet the needs of some special devices.
The evolutionary algorithm has the advantages which are incomparable with the traditional optimization algorithm, such as good robustness, no need of strict object structure and derivative information, etc. In recent years, application of evolutionary algorithms to optimization of rectification operations has become a hotspot of research. The genetic algorithm is used for cost optimization of a simple rectification process, such as by optimizing the initial population, resulting in satisfactory operating parameters. However, the evolutionary algorithm always has the contradiction of more times of calling the target function, low efficiency and difficulty in meeting real-time optimization. For the above problems, replacing the process mechanism model with a proxy model is an effective solution. For example, an Artificial Neural Network (ANN) can be used for simulating an atmospheric and vacuum rectifying tower and optimizing the operation conditions so as to improve the overall economic benefit; the evolutionary design of the rectification process can be realized by means of a neural network, and the optimization of the operation parameters of the rectification process can be realized by combining an expert system; there are also documents in which an atmospheric tower model established by ASPEN PLUS is proxied with a support vector machine to complete the optimum design of the atmospheric tower; aiming at the fact that the agent model is not updated in a point adding mode in the optimizing process, the initial sample number and the selection of sample points influence the optimizing effect, literature proposes a waste heat recovery process of an atmospheric tower based on Gaussian process agent model optimization, the agent model is updated in a point adding mode in the iteration process, and the result shows that the time consumption model calling times and optimizing time are obviously reduced while the optimizing target is achieved.
However, the optimization method technology related in the existing scientific research literature and patent materials is mostly limited to the boundary constraint problem of decision variables of the atmospheric and vacuum process, and the equality and inequality constraints for ensuring working conditions and safety requirements are not considered in the optimization process.
Disclosure of Invention
The main technical problems to be solved by the invention are as follows: aiming at the problems of time-consuming calculation and equation and inequality constraint optimization for guaranteeing working condition stability and safety requirements in the atmospheric and vacuum rectification process, how to implement optimization aiming at the maximum profit. Specifically, the method adopts a point adding strategy in optimization iteration, updates and perfects global and local agent models, performs optimization search on the built global and local agent models, thereby obtaining an optimized solution, and finally applies the optimized solution to a constraint evolution optimization process of profit maximization in an atmospheric and vacuum process.
The technical scheme adopted for solving the technical problems is as follows: a constraint evolution optimization method of an atmospheric and vacuum device based on a proxy model comprises the following steps:
generating N input and output sampling data by using a mechanism model of an atmospheric and vacuum rectification process, and recording an input data matrix as X epsilon R α×N The output data matrix is recorded as Y epsilon R β×N Wherein the α input measurement variables in the input data matrix X include: m is m 1 Feed flow, m of seed crude oil 2 The feed temperature, m, of the individual rectification columns 3 Reflux temperature and m 3 Reflux ratio, m 4 Mixing ratio of individual crude oil mixers, α=m 1 +m 2 +2m 3 +m 4 The beta output measurement variables in the output data matrix Y are the flow rates of beta product oils respectively, R α×N Representing a real matrix of dimension alpha x N, R β×N Representing a real matrix in the β×n dimension.
Step (2) establishes a GPR model between the input data matrix X and the beta output measurement variables, respectively, using a Gaussian process regression (Gaussian Process Regression, abbreviation: GPR) algorithm: y is Y b =f b (X)+E b Wherein f b Representing the nonlinear transformation process of the b-th GPR model, E b Is all that isValue zero, variance sigma b Gaussian noise of Y b To output the b-th row vector in the data matrix Y, b=1, 2, …, β.
The principle of the GPR algorithm is described as follows:
model parameters that the GPR algorithm needs to determine include: variance sigma b Nuclear parameter c b And eta b 。
First, a kernel covariance matrix C ε R is calculated according to the following formula N×N The ith row and the jth column of element C ij :
In the above, x i And x j Column vectors for the ith and jth columns, respectively, in the input data matrix X, when i=j,when i+.j, +.>
Next, a maximum likelihood function L is calculated according to the following formula b :
In the above expression, |c| represents a determinant of the calculation matrix C.
Then, the maximum likelihood function is recalculated with respect to the model parameter set Θ b ={σ b ,c b ,η b Partial derivative of }:
finally, the optimal solution obtained by carrying out maximized solution on the partial derivative in the formula (3) by using a conjugate gradient method is the model parameter set theta of the GPR algorithm b 。
Determining model parameter set theta of GPR algorithm b ={σ b ,c b ,η b After } the arbitrary input data vector z E R can be obtained α×1 Predicting the corresponding output value:wherein f b The represented nonlinear transformation process is as follows:
first, a kernel covariance vector k ε R is calculated according to the following formula N×1 :
In the above, k i Is the i-th element in the kernel covariance vector k.
Next, according to the formulaCalculating the output predictive value +.>At the same time, the GPR algorithm also gives the variance cov (z) =c of the prediction error z -k T C -1 k, wherein C z =η b +σ b 。
Step (3) determining a production profit maximization objective function and corresponding constraint conditions thereof as defined in the following formula (5):
in the above formula, x is R α×1 The data composition of alpha input measured variables, J (x) represents the production profit corresponding to x as input,Is recorded by m in x 1 Feed of seed crude oilFlow data composition, F (x) ε R β×1 Expressed as x ε R α×1 Beta product oil flows are input and calculated according to a mechanism model of the atmospheric and vacuum rectification process, and the beta product oil flows are ∈10>And->Unit price, x of crude oil and product oil respectively min And x max Respectively represent the upper limit and the lower limit of the input data vector, h p (x) And g is equal to q (x) The P-th and Q-th equality constraints, the number of equality constraints for P, the number of inequality constraints for Q, and the abbreviation for Subject To for s.t. are expressed respectively, meaning of the constraints.
Since equation (5) defines an optimization problem with multiple constraint conditions, if the particle swarm optimization algorithm is used for optimization solution, an equivalent unconstrained optimization problem is defined according to equation (5), namely:
in the above-mentioned method, the step of,for new objective function, γ=1, 2, …, (p+q), J 0 =max{J(x 1 ),J(x 2 ),…,J(x N ) The expression J (x) 1 ),J(x 2 ),…,J(x N ) Maximum value of (A), lambda is penalty factor, V γ (x) Representing a gamma penalty function, defined as follows:
step (4) firstly changing J (x) into:wherein->Then using the formula (6) as an optimization target, and solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors s 1 ,s 2 ,…,s β Wherein s is b ∈R α×1 The specific implementation process is as follows:
step (4.1): after initializing b=1, setting parameters of a particle swarm optimization algorithm, specifically including: total number of particles D, maximum number of iterations Im, acceleration factor ε 1 And epsilon 2 Generally, im=600, ε is preferable 1 =ε 2 =2。
Step (4.2): setting the iteration number iter=1, and in the interval x min ,x max ]Up-randomly generating D particle vectorsWherein->Represents the D-th particle vector, d=1, 2, …, D.
Step (4.3): calculating the respective particle vectors according to the above formula (6)Corresponding objective function valueWherein the penalty factor λ is according to the formula λ= |j 0 |·χ d Calculated by/(P+Q), χ d The number of times the d-th particle vector violates the constraint in equation (5) is expressed.
Step (4.4): will beThe particle vector corresponding to the maximum objective function value in the (a) is recorded as a vector t, and the maximum value obtained by each particle vector in the whole iterative process is obtainedParticle vectors corresponding to the objective function values are recorded as +.>And updating the running speed vector of each particle according to the following formula>
In the above, vector v d ∈R α×1 Each element is a section [ -1,1 []Random number, rand on 1 And rand 2 Is in interval [0,1 ]]Any random number within.
Step (4.5): vector pairThe correction is carried out on each element: if->If the element is greater than 1, modifying the element to be 1; if->If the element in the element is smaller than-1, modifying the element to be-1; in other cases, no modifications are made to the elements.
Step (4.6): according to the formulaUpdating individual particle vectors +.>And correcting the elements in each particle vector to be in the interval x min ,x max ]Inner part
Step (4.7): judging whether the condition is satisfied: item < Im? If yes, returning to the step (4.3) after setting item=item+1; if not, then the mostThe resulting optimal particle vector s b Is thatAnd (4.8) executing the step (4.8) according to the particle vector corresponding to the maximum value.
Step (4.8): judging whether the condition is satisfied: b < β? If yes, setting b=b+1, and returning to the step (4.2); if not, beta optimal particle vectors s are obtained 1 ,s 2 ,…,s β 。
Step (5) is respectively carried out by s 1 ,s 2 ,…,s β To input data vector and calculate F(s) according to the mechanism model of the atmospheric and vacuum device 1 ),F(s 2 ),…,F(s β ) Then, s is again 1 ,s 2 ,…,s β Added to the input data matrix, i.e. x= [ X, s 1 ,s 2 ,…,s β ]And correspondingly F(s) 1 ),F(s 2 ),…,F(s β ) Is added into the output data matrix, i.e. y= [ Y, F(s) 1 ),F(s 2 ),…,F(s β )]。
Step (6) calculates a crowding distance dist (i) of each column vector in the input data matrix X according to the following formula:
in the above formula, i=1, 2, …, N, min { } represents the minimum value of all elements in { }, the symbol is represented by the length of the vector is calculated and, N represents the total number of column vectors in the input data matrix X, N being changed continuously with the addition of data vectors.
Step (7) all column vectors satisfying the condition dist (i) > delta in the input data matrix X are formed into a global input data matrix X G And correspondingly forms the column vectors of the same columns in the output data matrix into a global output data matrix Y G And respectively establish X G And Y is equal to G GPR model between the β output measurement variables:wherein->Representing the nonlinear transformation process of the b-th GPR model,/->Is Gaussian noise->For outputting data matrix Y G Line b of (b), b=1, 2, …, β, the reference symbol G being the initial of the word Global, means Global.
Step (8) updates the objective function J (x) as follows:wherein the method comprises the steps ofNamely, the GPR model in the step (7) is used for substituting the mechanism model to calculate x epsilon R α×1 For the beta product oil flows corresponding to the input data, solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors phi 1 ,φ 2 ,…,φ β The specific implementation process is the same as that of the steps (4.1) to (4.8).
Step (9) is respectively carried out by phi 1 ,φ 2 ,…,φ β To input data vector and calculate F (phi) according to the mechanism model of the atmospheric and vacuum device 1 ),F(φ 2 ),…,F(φ β ) Then phi is again arranged 1 ,φ 2 ,…,φ β Added to the input data matrix, i.e. x= [ X, phi 1 ,φ 2 ,…,φ β ]And correspondingly F (phi) 1 ),F(φ 2 ),…,F(φ β ) Added to the output data matrix, i.e. y= [ Y, F (phi) 1 ),F(φ 2 ),…,F(φ β )]。
Step (10) X and Y b Is combined into a matrix Z b =[X T ,Y b T ] T ∈R (α+1)×N Then, the data vector phi is again b And F (phi) b ) The b-th element of the list is combined into a column vectorAnd from matrix Z b Find the AND vector from the column vectors of (1)>Nearest N b The column vectors, thus forming the b-th local data matrix +.>
Step (11) willThe row vectors of the 1 st to the alpha st row form a local input matrix +.>Will->The row vectors of the last row of the rows constitute a local output vector +.>And establish->And->RBF neural network model in between:wherein (1)>Nonlinear transformation procedure representing the b-th RBF neural network model, < >>For error vectors, the specific implementation procedure is as follows:
step (11.1) after setting the hidden layer node number as xi, randomly slave matrixAnd selecting xi column vectors as initial center point vectors of the clusters respectively.
Step (11.2) calculating a matrixAnd dividing the column vectors into corresponding clusters according to the minimum distance value.
And (11.3) calculating the mean value vector of all the home column vectors in each cluster, wherein the vector is the new center point vector of each cluster.
Step (11.4) determining whether each center point vector converges? If not, returning to the step (11.2); if yes, recording the converged central point vector as O 1 ,O 2 ,…,O ζ And step (11.5) is performed.
Where v=1, 2, …, ζ.
Step (11.6) calculating the input vector z εR according to the following formula α×1 Output u transformed by hidden layer v-th neuron node v :
Then the output vector of z transformed by all ζ hidden layer neuron nodes is u= [ u ] 1 ,u 2 ,…,u ξ ] T 。
Step (11.7) matrixEach row vector in (2) is taken as an input vector, and the output vector of each input vector converted by hidden layer neuron nodes is calculated according to the step (11.6)>Obtaining hidden layer output matrix
Step (11.8) according to the formulaCalculating hidden layer output matrix U b To the output layer->Regression coefficient vector B between b 。
Step (11.9) calculating an estimated value of the RBF neural network model pair outputAnd error matrix
Step (12) updates the objective function J (x) as follows:wherein->Namely, the RBF neural network model in the step (11) is used for replacing the mechanism model to calculate x epsilon R α×1 For the beta product oil flows corresponding to the input data, solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors psi 1 ,ψ 2 ,…,ψ β The specific implementation process is the same as that of the steps (4.1) to (4.8).
Step (13): respectively at psi 1 ,ψ 2 ,…,ψ β To input data vector and calculate F (psi) according to the mechanism model of atmospheric and vacuum rectification process 1 ),F(ψ 2 ),…,F(ψ β ) Then, psi is added again 1 ,ψ 2 ,…,ψ β Added to the input data matrix X, i.e. x= [ X, phi ] 1 ,φ 2 ,…,φ β ]And correspondingly F (phi) 1 ),F(φ 2 ),…,F(φ β ) Added to the output data matrix Y, i.e. y= [ Y, F (phi) 1 ),F(φ 2 ),…,F(φ β )]。
Step (14): and (3) repeating the steps (6) to (13) for n times, and ending the whole optimization process.
Compared with the traditional method, the method has the advantages that:
first, the method of the present invention considers the constraint problems of equations and inequalities in the optimization problem, as shown in equation (5). The constraint conditions of the equations and the inequality are determined according to the stability and the safety requirements of working conditions, so that the method has the greatest advantage over the traditional method in that the aimed optimization problem is more consistent with the actual running condition. In addition, the method of the invention searches global and local training data through newly defined optimization objective functions in a reciprocating and optimized way, thereby better ensuring the diversity of the training data and further continuously improving the accuracy of the proxy model.
Drawings
Figure 1 shows three crude oil rectification systems in a refinery.
FIG. 2 is a comparison of the results of the method of the present invention with the implementation of an optimization based on a pure mechanism model.
Detailed Description
The process according to the invention is described in detail below with reference to the drawings and to the specific examples. Fig. 1 shows a rectification system of three sets of crude oil in a refinery, which consists of a crude oil blending device, an atmospheric furnace AF, an atmospheric rectification tower ADT, a vacuum tower VDT, a heat exchange network NET and the like, wherein the sets of crude oil rectification devices CDU1 and CDU2 share the heat exchange network NET1. As shown in FIG. 1, 5 paths of desalted and dehydrated crude oil F1-F5 are mixed by a mixer according to a mixing ratio c1-c5 and then redistributed, heat exchange is carried out between the desalted and dehydrated crude oil F1-F5 and a distillation product with higher temperature through a heat exchange network, a tube type heating furnace is heated to a certain temperature, at the moment, part of the crude oil is gasified, and oil vapor and unvaporized oil enter an atmospheric tower and a vacuum tower in turn through an oil transfer line and are divided into naphtha, straight run distillate (SR-Jet), straight run distillate Diesel (SR-Diesel), atmospheric Gas Oil (AGO) and atmospheric residuum (ATMRES).
Thus, in the present embodiment the number of measured variables α=19 is entered, including m 1 Feed flow of =5 crude oils, m 2 Feed temperature of=3 rectifying columns, m 3 =3 reflux temperatures, m 3 =3 reflux ratios, m 4 Mixing ratio δ of =5 crude oil mixers 1 ,δ 2 ,…,δ 5 Then α=m 1 +m 2 +2m 3 +m 4 The number of output measurement variables β=5. The invention discloses a constraint evolution optimization method of an atmospheric and vacuum device based on a proxy model, and the specific implementation mode is as follows.
Step (1): generating N=90 input/output sampling data by using a mechanism model of an atmospheric and vacuum rectification process, and recording an input data matrix as X epsilon R α×N The output data matrix is recorded as Y epsilon R β×N 。
Step (2): a Gaussian Process Regression (GPR) algorithm is used to build a GPR model between an input data matrix X and β output measurement variables, respectively: y is Y b =f b (X)+E b 。
Step (3): the following formula is determinedThe production profit maximization objective function and the corresponding constraint condition, wherein gamma 1 ,γ 2 ,γ 3 The limiting parameters of the three mixers in fig. 1, respectively. />
Step (4): the objective function J (x) is first changed to:wherein the method comprises the steps ofThen using the formula (6) as an optimization target, and solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors s 1 ,s 2 ,…,s β The parameters of the particle swarm optimization algorithm are set as follows: population number d=80, maximum iteration number im=600, acceleration factor epsilon 1 =ε 2 =2。
Step (5): respectively in s 1 ,s 2 ,…,s β F(s) is calculated according to a mechanism model of the atmospheric and vacuum rectification process and is input with a data vector 1 ),F(s 2 ),…,F(s β ) Then, s is again 1 ,s 2 ,…,s β Added to the input data matrix, i.e. x= [ X, s 1 ,s 2 ,…,s β ]。
Step (6): the crowding distance dist (i) of each column vector in the input data matrix X is calculated.
Step (7): all column vectors meeting the condition dist (i) > delta in the input data matrix X are formed into a global input data matrix X G And correspondingly forms the corresponding column vectors in the output data matrix into a global output data matrix Y G And respectively establish X G And Y is equal to G GPR model between the β output measurement variables:
step (8): the objective function J (x) is updated first as:wherein the method comprises the steps ofThat is, the Gaussian process regression model in the step (7) is used for substituting the mechanism model to calculate x epsilon R α×1 For the beta product oil flows corresponding to the input data, solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors phi 1 ,φ 2 ,…,φ β 。
Step (9): respectively phi is marked by 1 ,φ 2 ,…,φ β F (phi) is calculated according to a mechanism model of the atmospheric and vacuum rectification process for inputting data vectors 1 ),F(φ 2 ),…,F(φ β ) Then phi is again arranged 1 ,φ 2 ,…,φ β Added to the input data matrix X and F (phi) is correspondingly applied 1 ),F(φ 2 ),…,F(φ β ) Added to the output data matrix Y.
Step (10): x and Y are combined b Is combined into a matrix Z b =[X T ,Y b T ] T ∈R (α+1)×N Then, the data vectors phi are respectively processed b And F (phi) b ) The b-th element of the list is combined into a column vectorAnd from matrix Z b Find the AND vector from the column vectors of (1)>Nearest N b The column vectors, thus forming the b-th local data matrix +.>
Step (11): will beLine vectors of 1 st line to alpha st line in (a)Composing a local input matrix->Will->The row vectors of the last row of the rows constitute a local output vector +.>And establish->And->RBF neural network model in between:
step (12): the objective function J (x) is updated first as:wherein the method comprises the steps ofNamely, the RBF neural network model in the step (11) is used for replacing the mechanism model to calculate x epsilon R α×1 For the beta product oil flows corresponding to the input data, solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors psi 1 ,ψ 2 ,…,ψ β 。/>
Step (13): respectively at psi 1 ,ψ 2 ,…,ψ β To input data vector and calculate F (psi) according to the mechanism model of atmospheric and vacuum rectification process 1 ),F(ψ 2 ),…,F(ψ β ) Then, psi is added again 1 ,ψ 2 ,…,ψ β Added to the input data matrix X and F (phi) is correspondingly applied 1 ),F(φ 2 ),…,F(φ β ) Added to the output data matrix Y.
Step (14): after repeating steps (6) to (13) n=400 times, the whole optimization process is ended.
Finally, the optimization results of the implementation steps are shown in table 1, and the obtained profit and the call number comparison of the mechanism model are shown in fig. 2. Table 1 lists the optimal values of the input measured variables before optimization, based on pure mechanism model optimization, and given by the method of the present invention. As can be seen from FIG. 2, the unit profit of the first three crude towers is optimized to be only about 3X 10 5 Through the algorithm of the invention, the economic profit of 5 crude oils is improved by times after the operation parameters of the device are optimized, and the unit profit can reach 6.007 multiplied by 10 5 If pure mechanism constraint optimization is adopted, the unit profit obtained after optimization is 6.188 multiplied by 10 5 The distribution of main operation parameters of the optimized algorithm is relatively close; from the aspect of the calling times of the mechanism model, the calling times of the optimization algorithm of the invention to the mechanism model are controlled to be 400 times, and 9263 times are needed for obtaining the equivalent economic benefit of the algorithm by adopting pure mechanism constraint optimization. Therefore, the multi-agent constraint optimization algorithm provided by the invention can obviously reduce the evaluation times of the time-consuming calculation model under the condition of obtaining a better optimization effect.
TABLE 1 crude oil blending coefficient and main operating parameter optimization results comparison
The above embodiments are merely illustrative of specific implementations of the invention and are not intended to limit the invention. Any modification made to the present invention that comes within the spirit of the present invention and the scope of the appended claims falls within the scope of the present invention.
Claims (3)
1. The constrained evolution optimization method of the atmospheric and vacuum device based on the agent model is characterized by comprising the following steps of:
step (1): generating N input/output sampling numbers by using mechanism model of atmospheric and vacuum deviceAccording to the data, the input data matrix is recorded as X epsilon R α×N The output data matrix is recorded as Y epsilon R β×N Wherein the α input measurement variables in the input data matrix X include: m is m 1 Feed flow, m of seed crude oil 2 The feed temperature, m, of the individual rectification columns 3 Reflux temperature and m 3 Reflux ratio, m 4 Mixing ratio of individual crude oil mixers, α=m 1 +m 2 +2m 3 +m 4 The beta output measurement variables in the output data matrix Y are the flow rates of beta product oils respectively, R α×N A real matrix representing the dimension α×n;
step (2): a GPR model between an input data matrix X and beta output measurement variables is established by using a GPR algorithm: y is Y b =f b (X)+E b Wherein f b Representing the nonlinear transformation process of the b-th GPR model, E b Mean value is zero and variance is sigma b Gaussian noise of Y b For the b-th row vector in the output data matrix Y, b=1, 2, …, β;
step (3): determining a production profit maximization objective function and corresponding constraint conditions defined in the following formula (1):
in the above formula, x is R α×1 The data composition of alpha input measured variables, J (x) represents the production profit corresponding to x as input,Is recorded by m in x 1 Feed flow data composition, F (x), ε R, of crude oil β×1 Expressed as x ε R α×1 Beta product oil flows are input and calculated according to a mechanism model of the atmospheric and vacuum rectification process, and the beta product oil flows are ∈10>And->Unit price, x of crude oil and product oil respectively min And x max Respectively represent the upper limit and the lower limit of the input data vector, h p (x) And g is equal to q (x) The P-th and Q-th equality constraints, the number of equality constraints for P, the number of inequality constraints for Q, and the abbreviation for Subject To for s.t. are expressed respectively, meaning of the constraints;
step (4): first, J (x) is changed into:wherein->And taking the following formula (2) as an optimization target:
in the above-mentioned method, the step of,for new objective function, γ=1, 2, …, (p+q), J 0 =max{J(x 1 ),J(x 2 ),…,J(x N ) The expression J (x) 1 ),J(x 2 ),…,J(x N ) Maximum value of (A), lambda is penalty factor, V γ (x) Representing a gamma penalty function, defined as follows:
and optimizing by using particle swarmAlgorithm solution to obtain beta optimal particle vectors s 1 ,s 2 ,…,s β Wherein s is b ∈R α×1 The specific implementation process is as follows:
step (4.1): after initializing b=1, setting parameters of a particle swarm optimization algorithm, specifically including: total number of particles D, maximum number of iterations Im, acceleration factor ε 1 And epsilon 2 ;
Step (4.2): setting the iteration number iter=1, and in the interval x min ,x max ]Up-randomly generating D particle vectorsWherein->Represents the D-th particle vector, d=1, 2, …, D;
step (4.3): calculating respective particle vectors according to the above formula (2)Corresponding objective function valueWherein the penalty factor λ is according to the formula λ= |j 0 |·χ d Calculated by/(P+Q), χ d Representing the number of times the d-th particle vector does not satisfy the constraint condition in the formula (1);
step (4.4): will beThe particle vector corresponding to the maximum objective function value in the iteration process is recorded as a vector t, and the particle vector corresponding to the maximum objective function value obtained by each particle vector in the whole iteration process is recorded asAnd updating the running speed of each particle according to the following formulaDegree vector->
In the above, vector v d ∈R α×1 Each element is a section [ -1,1 []Random number, rand on 1 And rand 2 Is in interval [0,1 ]]Random numbers within;
step (4.5): vector pairThe correction is carried out on each element: if->If the element is greater than 1, modifying the element to be 1; if->If the element in the element is smaller than-1, modifying the element to be-1; in other cases, no modifications are made to the elements;
step (4.6): according to the formulaUpdating individual particle vectors +.>And correcting the elements in each particle vector to be in the interval x min ,x max ]An inner part;
step (4.7): judging whether the condition is satisfied: item < Im? If yes, returning to the step (4.3) after setting item=item+1; if not, the optimal particle vector s is finally obtained b Is thatThe particle vector corresponding to the maximum value in (2) and executing the step (4.8);
step (4.8): judging whether the condition is satisfied: b < β? If yes, setting b=b+1, and returning to the step (4.2); if not, beta optimal particle vectors s are obtained 1 ,s 2 ,…,s β ;
Step (5): respectively in s 1 ,s 2 ,…,s β To input data vector and calculate F(s) according to the mechanism model of the atmospheric and vacuum device 1 ),F(s 2 ),…,F(s β ) Then, s is again 1 ,s 2 ,…,s β Added to the input data matrix X, i.e. x= [ X, s 1 ,s 2 ,…,s β ]And correspondingly F(s) 1 ),F(s 2 ),…,F(s β ) Added to the output data matrix Y, i.e. y= [ Y, F(s) 1 ),F(s 2 ),…,F(s β )];
Step (6): the crowding distance dist (i) for each column vector in the input data matrix X is calculated according to the formula:
in the above formula, i=1, 2, …, N, min { } represents the minimum value of all elements in { }, the symbol is represented by the length of the vector is calculated and, N represents the total number of column vectors in the input data matrix X, and N is changed continuously along with the operation of adding the data vectors;
step (7): all column vectors meeting the condition dist (i) > delta in the input data matrix X are formed into a global input data matrix X G And correspondingly forms the column vectors of the same column in the output data matrix Y into a global output data matrix Y G And respectively establish X G And Y is equal to G GPR model between the β output measurement variables:wherein->Representing the nonlinear transformation process of the b-th GPR model,/->Is Gaussian noise->For outputting data matrix Y G B=1, 2, …, β, δ is the distance threshold, G is the initial of the word Global, meaning globally;
step (8): first, J (x) is updated as follows:wherein->Then solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors phi 1 ,φ 2 ,…,φ β The specific implementation process is the same as that of the steps (4.1) to (4.8);
step (9): respectively phi is marked by 1 ,φ 2 ,…,φ β For inputting data vector, and respectively calculating to obtain vector F (phi) according to mechanism model of atmospheric and vacuum device 1 ),F(φ 2 ),…,F(φ β ) Then phi is again arranged 1 ,φ 2 ,…,φ β Added to the input data matrix X, i.e. x= [ X, phi ] 1 ,φ 2 ,…,φ β ]And correspondingly F (phi) 1 ),F(φ 2 ),…,F(φ β ) Added to the output data matrix Y, i.e. y= [ Y, F (phi) 1 ),F(φ 2 ),…,F(φ β )];
Step (10): x and Y are combined b Is combined into a matrix Z b =[X T ,Y b T ] T ∈R (α+1)×N Then, the data vector phi is again b Vector F (phi) b ) The b-th element of the list is combined into a column vectorAnd from matrix N b Find the AND vector from the column vectors of (1)>Nearest N b The column vectors, thus forming the b-th local data matrix +.>
Step (11): will beThe row vectors of the 1 st to the alpha st row form a local input matrix +.>Will->The row vectors of the last row of the rows constitute a local output vector +.>And establish->And->RBF neural network model in between: />Wherein (1)>Nonlinear transformation procedure representing the b-th RBF neural network model, < >>Is an error vector;
step (12): first, J (x) is updated as follows:wherein->Then solving by using a particle swarm optimization algorithm to obtain beta optimal particle vectors psi 1 ,ψ 2 ,…,ψ β The specific implementation process is the same as that of the steps (4.1) to (4.8);
step (13): respectively at psi 1 ,ψ 2 ,…,ψ β To input data vector and calculate F (psi) according to the mechanism model of the atmospheric and vacuum device 1 ),F(ψ 2 ),…,F(ψ β ) Then, psi is added again 1 ,ψ 2 ,…,ψ β Added to X, i.e. X= [ X, phi ] 1 ,φ 2 ,…,φ β ]And correspondingly F (phi) 1 ),F(φ 2 ),…,F(φ β ) Added to the output data matrix Y, i.e. y= [ Y, F (phi) 1 ),F(φ 2 ),…,F(φ β )];
Step (14): and (3) repeating the steps (6) to (13) for n times, and ending the whole optimization process.
2. The constrained evolution optimization method of the atmospheric and vacuum device based on the proxy model according to claim 1, wherein the detailed implementation process of the GPR model between the input data matrix X and the β output measurement variables in the step (2) is specifically as follows:
first, a kernel covariance matrix C ε R is calculated according to the following formula N×N The ith row and the jth column of element C ij :
In the above, x i And x j Column vectors for the ith and jth columns, respectively, in the input data matrix X, when i=j,when i+.j, +.>i=1,2,…,N,j=1,2,…,N;
Next, a maximum likelihood function L is calculated according to the following formula b :
In the above expression, |c| represents a determinant of the calculation matrix C;
then, the maximum likelihood function is recalculated with respect to the model parameter set Θ b ={σ b ,c b ,η b Partial derivative of }:
finally, the optimal solution obtained by carrying out maximized solution on the partial derivative in the formula (3) by using a conjugate gradient method is the model parameter set theta of the GPR algorithm b ;
Determining model parameter set theta of GPR algorithm b ={σ b ,c b ,η b After } the arbitrary input data vector z E R can be obtained α×1 Predicting the corresponding output value:wherein f b Represented byThe nonlinear transformation process of (2) is as follows:
first, a kernel covariance vector k ε R is calculated according to the following formula N×1 :
In the above, k i Is the i-th element in the kernel covariance vector k;
3. The constrained evolution optimization method of an atmospheric and vacuum apparatus based on a proxy model according to claim 1, wherein the step (11) is establishedAnd->The detailed implementation process of the RBF neural network model comprises the following steps:
step (11.1) after setting the hidden layer node number as xi, randomly slave matrixSelecting xi column vectors as initial center point vectors of each cluster respectively;
step (11.2) calculating a matrixThe distance between each column vector and the xi central point vectors, and dividing the column vector into corresponding clustering clusters according to the minimum value of the distance;
step (11.3) calculating the mean value vector of all the home column vectors in each cluster, wherein the vector is the new central point vector of each cluster;
step (11.4) determining whether each center point vector converges? If not, returning to the step (11.2); if yes, recording the converged central point vector as O 1 ,O 2 ,…,O ξ And performing step (11.5);
step (11.5) calculating RBF parameter θ according to the following formula υ :
Wherein v=1, 2, …, ζ;
step (11.6) calculating the input vector z εR according to the following formula α×1 Output u transformed by hidden layer upsilon neuron nodes υ :
Then the output vector of z transformed by all ζ hidden layer neuron nodes is u= [ u ] 1 ,u 2 ,…,u ξ ] T ;
Step (11.7) matrixEach row vector in (2) is taken as an input vector, and the output vector of each input vector converted by hidden layer neuron nodes is calculated according to the step (11.6)>Obtaining hidden layer output matrix->
Step (11.8) according to the formulaCalculating hidden layer output matrix U b To the output layer->Regression coefficient vector B between b ;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910408795.2A CN111914382B (en) | 2019-05-07 | 2019-05-07 | Constraint evolution optimization method of atmospheric and vacuum device based on agent model |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910408795.2A CN111914382B (en) | 2019-05-07 | 2019-05-07 | Constraint evolution optimization method of atmospheric and vacuum device based on agent model |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111914382A CN111914382A (en) | 2020-11-10 |
CN111914382B true CN111914382B (en) | 2023-04-21 |
Family
ID=73242984
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910408795.2A Active CN111914382B (en) | 2019-05-07 | 2019-05-07 | Constraint evolution optimization method of atmospheric and vacuum device based on agent model |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111914382B (en) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112465198B (en) * | 2020-11-16 | 2023-07-28 | 湖南大学 | Double-agent assisted search energy optimization method for small sample data of park |
CN113110060B (en) * | 2021-04-29 | 2023-01-06 | 中海石油炼化有限责任公司 | Real-time optimization method of reforming device |
CN113408194B (en) * | 2021-06-10 | 2022-12-02 | 北京宜能高科科技有限公司 | General disc optimization method of atmospheric and vacuum distillation unit |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109543263A (en) * | 2018-11-01 | 2019-03-29 | 宁波大学 | A kind of method for building up of integrated atmospheric distillation process agent model |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8670960B2 (en) * | 2010-03-16 | 2014-03-11 | Schlumberger Technology Corporation | Proxy methods for expensive function optimization with expensive nonlinear constraints |
US9733629B2 (en) * | 2014-07-21 | 2017-08-15 | Honeywell International Inc. | Cascaded model predictive control (MPC) approach for plantwide control and optimization |
-
2019
- 2019-05-07 CN CN201910408795.2A patent/CN111914382B/en active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109543263A (en) * | 2018-11-01 | 2019-03-29 | 宁波大学 | A kind of method for building up of integrated atmospheric distillation process agent model |
Also Published As
Publication number | Publication date |
---|---|
CN111914382A (en) | 2020-11-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111914382B (en) | Constraint evolution optimization method of atmospheric and vacuum device based on agent model | |
Yu et al. | A combined genetic algorithm/simulated annealing algorithm for large scale system energy integration | |
Ryberg et al. | Metamodel-based multidisciplinary design optimization for automotive applications | |
Shelokar et al. | An ant colony classifier system: application to some process engineering problems | |
Yuzgec et al. | Dynamic neural-network-based model-predictive control of an industrial baker's yeast drying process | |
CN107070704A (en) | A kind of Trusted Web services combined optimization method based on QoS | |
Bi et al. | Proximal alternating-direction-method-of-multipliers-incorporated nonnegative latent factor analysis | |
CN107045569B (en) | Gear reducer optimization design method based on clustering multi-target distribution estimation algorithm | |
CN112489733A (en) | Octane number loss prediction method based on particle swarm algorithm and neural network | |
CN106600001B (en) | Glass furnace Study of Temperature Forecasting method based on Gaussian mixtures relational learning machine | |
CN113128124B (en) | Multi-grade C-Mn steel mechanical property prediction method based on improved neural network | |
CN114239400A (en) | Multi-working-condition process self-adaptive soft measurement modeling method based on local double-weighted probability hidden variable regression model | |
CN111914381B (en) | Operation optimization method of atmospheric and vacuum device based on KPLSR model | |
CN109493921B (en) | Multi-agent model-based normal pressure rectification process modeling method | |
Liu et al. | Distributed optimization subject to inseparable coupled constraints: a case study on plant-wide ethylene process | |
CN110728031B (en) | Multi-objective optimization method for balancing complex petrochemical process production energy based on ANN modeling | |
CN109543263B (en) | Method for establishing integrated atmospheric distillation process proxy model | |
Liang et al. | Optimal data fusion based on information quality function | |
Kim et al. | Comparison of Derivative‐Free Optimization: Energy Optimization of Steam Methane Reforming Process | |
Nian et al. | Strategy of changing cracking furnace feedstock based on improved group search optimization | |
CN114187975A (en) | Method and system for optimizing and automatically updating data model of diesel hydrogenation device | |
Liu et al. | An improved quantum particle swarm algorithm for solving multi-objective fuzzy flexible job shop scheduling problem | |
CN109388062B (en) | Global coordination distributed predictive control algorithm based on system decomposition indexes | |
Xue et al. | Crude oil distillation optimization using surrogate-aided constrained evolutionary optimization | |
Li et al. | A Nonlinear Proportional Integral Derivative-Incorporated Stochastic Gradient Descent-based Latent Factor Model |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |