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 PDF

Info

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
Application number
CN201910408795.2A
Other languages
Chinese (zh)
Other versions
CN111914382A (en
Inventor
薛锋
史旭华
唐俊苗
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ningbo University
Original Assignee
Ningbo University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Ningbo University filed Critical Ningbo University
Priority to CN201910408795.2A priority Critical patent/CN111914382B/en
Publication of CN111914382A publication Critical patent/CN111914382A/en
Application granted granted Critical
Publication of CN111914382B publication Critical patent/CN111914382B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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

Constraint evolution optimization method of atmospheric and vacuum device based on agent model
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 ：
Figure BSA0000183217290000021
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,
Figure BSA0000183217290000022
when i+.j, +.>
Figure BSA0000183217290000023
Next, a maximum likelihood function L is calculated according to the following formula b ：
Figure BSA0000183217290000024
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 }:
Figure BSA0000183217290000031
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:
Figure BSA0000183217290000032
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 ：
Figure BSA0000183217290000033
In the above, k i Is the i-th element in the kernel covariance vector k.
Next, according to the formula
Figure BSA0000183217290000034
Calculating the output predictive value +.>
Figure BSA0000183217290000035
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):
Figure BSA0000183217290000036
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,
Figure BSA00001832172900000311
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>
Figure BSA0000183217290000037
And->
Figure BSA0000183217290000038
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:
Figure BSA0000183217290000039
in the above-mentioned method, the step of,
Figure BSA00001832172900000310
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:
Figure BSA0000183217290000041
step (4) firstly changing J (x) into:
Figure BSA0000183217290000042
wherein->
Figure BSA0000183217290000043
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 vectors
Figure BSA0000183217290000044
Wherein->
Figure BSA0000183217290000045
Represents the D-th particle vector, d=1, 2, …, D.
Step (4.3): calculating the respective particle vectors according to the above formula (6)
Figure BSA0000183217290000046
Corresponding objective function value
Figure BSA0000183217290000047
Wherein 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 be
Figure BSA0000183217290000048
The 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 +.>
Figure BSA0000183217290000049
And updating the running speed vector of each particle according to the following formula>
Figure BSA00001832172900000410
Figure BSA00001832172900000411
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 pair
Figure BSA00001832172900000417
The correction is carried out on each element: if->
Figure BSA00001832172900000412
If the element is greater than 1, modifying the element to be 1; if->
Figure BSA00001832172900000413
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 formula
Figure BSA00001832172900000414
Updating individual particle vectors +.>
Figure BSA00001832172900000415
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 that
Figure BSA00001832172900000416
And (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:
Figure BSA0000183217290000051
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:
Figure BSA0000183217290000052
wherein->
Figure BSA0000183217290000053
Representing the nonlinear transformation process of the b-th GPR model,/->
Figure BSA0000183217290000054
Is Gaussian noise->
Figure BSA0000183217290000055
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:
Figure BSA00001832172900000521
wherein the method comprises the steps of
Figure BSA0000183217290000056
Namely, 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 vector
Figure BSA0000183217290000057
And from matrix Z b Find the AND vector from the column vectors of (1)>
Figure BSA0000183217290000058
Nearest N b The column vectors, thus forming the b-th local data matrix +.>
Figure BSA0000183217290000059
Step (11) will
Figure BSA00001832172900000510
The row vectors of the 1 st to the alpha st row form a local input matrix +.>
Figure BSA00001832172900000511
Will->
Figure BSA00001832172900000512
The row vectors of the last row of the rows constitute a local output vector +.>
Figure BSA00001832172900000513
And establish->
Figure BSA00001832172900000514
And->
Figure BSA00001832172900000515
RBF neural network model in between:
Figure BSA00001832172900000516
wherein (1)>
Figure BSA00001832172900000517
Nonlinear transformation procedure representing the b-th RBF neural network model, < >>
Figure BSA00001832172900000518
For error vectors, the specific implementation procedure is as follows:
step (11.1) after setting the hidden layer node number as xi, randomly slave matrix
Figure BSA00001832172900000519
And selecting xi column vectors as initial center point vectors of the clusters respectively.
Step (11.2) calculating a matrix
Figure BSA00001832172900000520
And 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.
Step (11.5) calculating RBF parameters according to the following formula
Figure BSA00001832172900000612
Figure BSA0000183217290000061
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 ：
Figure BSA0000183217290000062
Then the output vector of z transformed by all ζ hidden layer neuron nodes is u= [ u ] 1 ，u 2 ，…，u ξ ] T 。
Step (11.7) matrix
Figure BSA0000183217290000063
Each 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)>
Figure BSA00001832172900000610
Obtaining hidden layer output matrix
Figure BSA00001832172900000611
Step (11.8) according to the formula
Figure BSA0000183217290000064
Calculating hidden layer output matrix U b To the output layer->
Figure BSA0000183217290000065
Regression coefficient vector B between b 。
Step (11.9) calculating an estimated value of the RBF neural network model pair output
Figure BSA0000183217290000066
And error matrix
Figure BSA0000183217290000067
Step (12) updates the objective function J (x) as follows:
Figure BSA0000183217290000069
wherein->
Figure BSA0000183217290000068
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 determined
Figure BSA0000183217290000074
The 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. />
Figure BSA0000183217290000071
Step (4): the objective function J (x) is first changed to:
Figure BSA0000183217290000072
wherein the method comprises the steps of
Figure BSA0000183217290000073
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 β 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:
Figure BSA0000183217290000081
step (8): the objective function J (x) is updated first as:
Figure BSA00001832172900000814
wherein the method comprises the steps of
Figure BSA0000183217290000082
That 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 vector
Figure BSA0000183217290000083
And from matrix Z b Find the AND vector from the column vectors of (1)>
Figure BSA0000183217290000084
Nearest N b The column vectors, thus forming the b-th local data matrix +.>
Figure BSA0000183217290000085
Step (11): will be
Figure BSA0000183217290000086
Line vectors of 1 st line to alpha st line in (a)Composing a local input matrix->
Figure BSA0000183217290000087
Will->
Figure BSA0000183217290000088
The row vectors of the last row of the rows constitute a local output vector +.>
Figure BSA0000183217290000089
And establish->
Figure BSA00001832172900000810
And->
Figure BSA00001832172900000811
RBF neural network model in between:
Figure BSA00001832172900000812
step (12): the objective function J (x) is updated first as:
Figure BSA00001832172900000815
wherein the method comprises the steps of
Figure BSA00001832172900000813
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 ，…，ψ β 。/>
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
Figure BSA0000183217290000091
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):
Figure FSA0000183217280000011
Figure FSA0000183217280000012
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,
Figure FSA0000183217280000013
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>
Figure FSA0000183217280000014
And->
Figure FSA0000183217280000015
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:
Figure FSA0000183217280000016
wherein->
Figure FSA0000183217280000017
And taking the following formula (2) as an optimization target:
Figure FSA0000183217280000018
in the above-mentioned method, the step of,
Figure FSA0000183217280000019
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:
Figure FSA00001832172800000110
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 vectors
Figure FSA0000183217280000021
Wherein->
Figure FSA0000183217280000022
Represents the D-th particle vector, d=1, 2, …, D;
step (4.3): calculating respective particle vectors according to the above formula (2)
Figure FSA0000183217280000023
Corresponding objective function value
Figure FSA0000183217280000024
Wherein 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 be
Figure FSA0000183217280000025
The 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 as
Figure FSA0000183217280000026
And updating the running speed of each particle according to the following formulaDegree vector->
Figure FSA0000183217280000027
Figure FSA0000183217280000028
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 pair
Figure FSA0000183217280000029
The correction is carried out on each element: if->
Figure FSA00001832172800000210
If the element is greater than 1, modifying the element to be 1; if->
Figure FSA00001832172800000211
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 formula
Figure FSA00001832172800000212
Updating individual particle vectors +.>
Figure FSA00001832172800000213
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 that
Figure FSA00001832172800000214
The 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:
Figure FSA00001832172800000215
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:
Figure FSA0000183217280000031
wherein->
Figure FSA0000183217280000032
Representing the nonlinear transformation process of the b-th GPR model,/->
Figure FSA0000183217280000033
Is Gaussian noise->
Figure FSA0000183217280000034
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:
Figure FSA0000183217280000035
wherein->
Figure FSA0000183217280000036
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 vector
Figure FSA0000183217280000037
And from matrix N b Find the AND vector from the column vectors of (1)>
Figure FSA0000183217280000038
Nearest N b The column vectors, thus forming the b-th local data matrix +.>
Figure FSA0000183217280000039
Step (11): will be
Figure FSA00001832172800000310
The row vectors of the 1 st to the alpha st row form a local input matrix +.>
Figure FSA00001832172800000311
Will->
Figure FSA00001832172800000312
The row vectors of the last row of the rows constitute a local output vector +.>
Figure FSA00001832172800000313
And establish->
Figure FSA00001832172800000314
And->
Figure FSA00001832172800000315
RBF neural network model in between: />
Figure FSA00001832172800000316
Wherein (1)>
Figure FSA00001832172800000317
Nonlinear transformation procedure representing the b-th RBF neural network model, < >>
Figure FSA00001832172800000318
Is an error vector;
step (12): first, J (x) is updated as follows:
Figure FSA00001832172800000319
wherein->
Figure FSA00001832172800000320
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 ：
Figure FSA00001832172800000321
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,
Figure FSA00001832172800000322
when i+.j, +.>
Figure FSA00001832172800000323
i＝1，2，…，N，j＝1，2，…，N；
Next, a maximum likelihood function L is calculated according to the following formula b ：
Figure FSA0000183217280000041
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 }:
Figure FSA0000183217280000042
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:
Figure FSA0000183217280000043
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 ：
Figure FSA0000183217280000044
In the above, k i Is the i-th element in the kernel covariance vector k;
next, according to the formula
Figure FSA0000183217280000045
Calculating the output predictive value +.>
Figure FSA0000183217280000046
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 。
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 established
Figure FSA0000183217280000047
And->
Figure FSA0000183217280000048
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 matrix
Figure FSA0000183217280000049
Selecting xi column vectors as initial center point vectors of each cluster respectively;
step (11.2) calculating a matrix
Figure FSA00001832172800000410
The 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 υ ：
Figure FSA00001832172800000411
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 υ ：
Figure FSA0000183217280000051
Then the output vector of z transformed by all ζ hidden layer neuron nodes is u= [ u ] 1 ，u 2 ，…，u ξ ] T ；
Step (11.7) matrix
Figure FSA0000183217280000052
Each 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)>
Figure FSA0000183217280000053
Obtaining hidden layer output matrix->
Figure FSA0000183217280000054
Step (11.8) according to the formula
Figure FSA0000183217280000055
Calculating hidden layer output matrix U b To the output layer->
Figure FSA0000183217280000056
Regression coefficient vector B between b ；
Step (11.9) calculating an estimated value of the RBF neural network model pair output
Figure FSA0000183217280000057
And error matrix
Figure FSA0000183217280000058
/>
CN201910408795.2A 2019-05-07 2019-05-07 Constraint evolution optimization method of atmospheric and vacuum device based on agent model Active CN111914382B (en)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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

Patent Citations (1)

* Cited by examiner, † Cited by third party
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
Hosseini et al. Hybrid imperialist competitive algorithm, variable neighborhood search, and simulated annealing for dynamic facility layout problem
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
Zhang et al. Multi-objective optimization for materials design with improved NSGA-II
Lin et al. Novel adaptive hybrid rule network based on TS fuzzy rules using an improved quantum-behaved particle swarm optimization
CN112489733A (en) Octane number loss prediction method based on particle swarm algorithm and neural network
CN113128124B (en) Multi-grade C-Mn steel mechanical property prediction method based on improved neural network
CN106600001B (en) Glass furnace Study of Temperature Forecasting method based on Gaussian mixtures relational learning machine
Bi et al. Proximal alternating-direction-method-of-multipliers-incorporated nonnegative latent factor analysis
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
Tian et al. Dynamic operation optimization based on improved dynamic multi-objective dragonfly algorithm in continuous annealing process.
Liang et al. Optimal data fusion based on information quality function
CN116415177A (en) Classifier parameter identification method based on extreme learning machine
Liu et al. Distributed optimization subject to inseparable coupled constraints: a case study on plant-wide ethylene process
CN114187975A (en) Method and system for optimizing and automatically updating data model of diesel hydrogenation device
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
Xue et al. Crude oil distillation optimization using surrogate-aided constrained evolutionary optimization
Nian et al. Strategy of changing cracking furnace feedstock based on improved group search optimization
CN111639797A (en) Gumbel-softmax technology-based combined optimization method
Liu et al. An improved quantum particle swarm algorithm for solving multi-objective fuzzy flexible job shop scheduling problem

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