Dynamic constraint multi-objective optimization method for catalytic cracking production indexes

Abstract

The invention provides a dynamic constraint multi-objective optimization method of a catalytic cracking production index, and relates to the technical field of intelligent optimization, a random weight neural network is pre-constructed to establish a relationship among a production index, production data, fuel oil yield and distillation temperature, and based on a fuel oil yield constraint, an equipment production constraint and a constraint boundary of the production index, a dynamic constraint multi-objective optimization method of the catalytic cracking production index is established. Optimizing the set value of the production index, specifically, adopting a dynamic response strategy based on population feasibility after the production constraint condition is changed, and obtaining a better set value of the population under the current constraint condition; obtaining a feasible solution and a non-feasible solution with a relatively small constraint violation degree and a relatively good target value by adopting a target value modification strategy and a parent selection strategy of constraint violation penalty; simulating binary crossover and polynomial mutation to optimize the production index to obtain an optimal set value of the production index; and finally, transmitting the optimized production index set value to a production site.

Description

Technical Field

[0001] The invention relates to the technical field of intelligent optimization, in particular to a dynamic constraint multi-objective optimization method for catalytic cracking production indicators. Background Art

[0002] In the face of global competition and environmental protection requirements, the problem of production index optimization has become increasingly important in industrial production. The ultimate goal of this problem is to determine the optimal set value of industrial production indicators based on the constraints in the actual production process. For production equipment, dynamic environments such as changes in raw material properties and fluctuations in production constraints will cause fluctuations in the set value of the optimal operating variable in the production process. Failure to adjust the set value of the operating variable in a timely manner will reduce product quality, the economic benefits of the enterprise, and even endanger the production safety of the enterprise.

[0003] The petroleum refining industry is one of the pillar industries of the world's energy industry, providing energy for many other important industries (such as the transportation industry). The fluid catalytic cracking unit (FCC) of petroleum is an important secondary processing equipment in the refining industry. It can crack long-chain hydrocarbon petroleum molecules into shorter-chain hydrocarbon products (such as gasoline, diesel, etc.) under the action of catalysts, and then obtain diesel, gasoline and other products. In the catalytic cracking unit, the value of production indicators (such as feed temperature, catalyst temperature, flow rate, etc.) determines the economic benefits of the enterprise and the output of the product, and is also constrained by production demand, unit production capacity, etc. In addition, in actual production, affected by market demand, container capacity, unit operating conditions, etc., the maximum production of different fuel oils changes over time, which in turn affects the economic value of the fuel oil, the energy consumption of the unit, etc. Therefore, it is of great significance to efficiently optimize the operating variables of the catalytic cracking unit in combination with the dynamic conditions and constraints in the production process to reduce production costs and improve the economic benefits of the enterprise.

[0004] In general, the production index optimization problem of the production process of the catalytic cracking unit mainly faces the following challenges:

[0005] 1. Dealing with infeasible solutions

[0006] In a catalytic cracking unit, due to the existence of constraints, even if the decision variable is within the range of its maximum and minimum values, it may violate constraints such as production, which in turn makes the variable infeasible, resulting in discontinuous target values ​​corresponding to the feasible variables, which may eventually cause its true frontier to be discontinuous. Improving some infeasible solutions can easily turn them into feasible solutions or even non-dominated solutions. In addition, when some feasible variables intersect with infeasible operating variables close to the true frontier, it may be easier to guide the feasible solution to move toward the frontier. Therefore, how to deal with infeasible solutions in the optimization process is a challenging problem.

[0007] 2. Design of dynamic response strategy

[0008] In a catalytic cracking unit, the production constraints, container capacity, etc. change over time, which will cause irregular changes in the feasible domain of decision variables and target values. The non-dominated solution at the previous moment may become an infeasible solution. At the same time, the infeasible solution in the previous environment may become a feasible solution or even a non-dominated solution at the current moment. In addition, the real POF changes with the changes in production constraints, showing a discontinuous and unpredictable distribution. However, refineries usually need to determine the best solution for fossil fuels in a short time (within two minutes). Therefore, it is a relatively large challenge to design a fast and effective dynamic response strategy for catalytic cracking units to cope with the irregular changes in the real POF.

[0009] Therefore, it is of great practical application value to study how to timely and effectively optimize the production indicators of each process in the production process of the catalytic cracking unit under a dynamic environment to improve product quality and output, thereby improving the economic benefits of the catalytic cracking enterprise. Summary of the invention

[0010] In view of the shortcomings of the prior art, the purpose of the present invention is to propose a dynamic constraint multi-objective optimization method for catalytic cracking production indicators, comprising:

[0011] Step 1: Population initialization;

[0012] Step 1.1: Obtaining a formula for production constraints of a catalytic cracking unit, wherein the formula for production constraints includes a formula for representing a threshold interval of all data in a production index, a formula for representing a threshold interval of all data in production data, and a formula for representing a constraint relationship between the production index and the production data;

[0013] The production indicators include: wax oil flow FHO, wax oil temperature THO, catalyst flow FC, catalyst temperature TC, lifting steam flow FLS, stripping steam flow FSS, lifting pipe outlet temperature TR, first circuit recovery heat HR1, second circuit recovery heat HR2, third circuit recovery heat HR3; the production data include gasoline output Y1, diesel output Y2, natural gas output Y3, pollutant output YP, and 10% distillation temperature DT of gasoline, diesel and natural gas s，10％ The 90% distillation temperature DT of gasoline, diesel and natural gas s，90％ ;

[0014] Step 1.2: Generate the initial production index based on the threshold range of all data in the production index in the production constraint conditions, which is calculated by the following formula:

[0015]

[0016] Among them, rand is a random number between [0, 1]. is the minimum value of the jth data in the production index, is the maximum value of the jth data in the production index;

[0017] Calculate and obtain N initial production indicators, use the N initial production indicators as N solutions, and the N solutions form an initial population. Set the initial iteration number Gen=1;

[0018] Step 2: Determine whether the production constraint conditions under the current iteration number environment are the same as the production constraint conditions of the previous iteration. If the production constraint conditions under the current iteration number environment are not the same as the production constraint conditions of the previous iteration, execute step 3. If the production constraint conditions under the current iteration number environment are the same as the production constraint conditions of the previous iteration, execute step 4. When the current iteration number is 1, execute step 4.

[0019] Step 3: Calculate the new population based on the current population;

[0020] Step 3.1: Generate a new constraint formula based on the production constraint formula, calculate the constraint violation value cv of each solution based on the new constraint formula, determine whether the solution of the current population is a feasible solution based on the constraint violation value cv of each solution, the constraint violation value cv of the feasible solution is 0, and obtain the number of feasible solutions in the current population;

[0021] Step 3.2: When the number of feasible solutions in the current population is greater than the preset threshold, execute steps 3.3 to 3.5. After step 3.5 is completed, execute step 4. When the number of feasible solutions in the current population is less than or equal to the preset threshold, execute steps 3.6 to 3.9. After step 3.9 is completed, execute step 4.

[0022] Step 3.3: Input the N production indicators corresponding to the N solutions in the current population into the model of the relationship between production indicators and production data to obtain N production data, and then calculate N production targets based on the N production indicators and N production data and the calculation formula of the production target. The production targets include maximizing the production profit f product and minimize production energy consumption f energy , calculate the maximum production profit f product The corresponding target value, and the minimization of production energy consumption f energy The corresponding target value is to maximize the production profit f product The corresponding target value and the minimized production energy consumption f energy The corresponding target values ​​constitute the target value of the penalty term;

[0023] Among them, minimize the production energy consumption f energy The corresponding target value is calculated by the following formula:

[0024]

[0025]

[0026] Where cv' is the normalized constraint violation value corresponding to the jth X(t), r f is the proportion of feasible solutions in the current population, f j '(X(t)) is the minimum production energy consumption f corresponding to the jth X(t) energy The normalized value, f j 1 (X(t)) is the minimum production energy consumption f energy The corresponding target value; X(t) is feasible, which means X(t) is a feasible solution; X(t) is infeasible, which means X(t) is an infeasible solution;

[0027] Among them, maximize the production profit f product The corresponding target value is calculated by the following formula:

[0028]

[0029]

[0030] Among them, f j ”(X(t)) is the maximum production profit f corresponding to the jth X(t) product The normalized value, f j 2 (X(t)) is the maximum production profit f corresponding to the jth X(t) product Normalized value;

[0031] Thus, the target values ​​of the N penalty items are obtained. According to the target values ​​of the N penalty items, the N solutions in the current population are sorted. After the sorting is completed, the N solutions in the current population are arranged in ascending order.

[0032] Step 3.4: Randomly select two solutions from the N solutions in the current population, and select the solution with the smaller sequence number as the parent solution S according to the sequence numbers corresponding to the two solutions. 1 , repeat the above operation to get the parent solution S 2 , solve S for the parent 1 and the parent solution S 2 Perform simulated binary crossover to obtain the offspring solution c 1 and the offspring solution c 2 , solve S for the parent 1 and the parent solution S2 Perform polynomial mutation to obtain offspring solution and offspring solution In the offspring solution c 1 、Offspring solution c 2 、Sub-generation solution and offspring solution When the number reaches N, the offspring population is generated according to the N offspring solutions;

[0033] Step 3.5: Input the production index corresponding to each child solution in the child population into the model of the relationship between the production index and the production data, calculate the production data of each child solution, calculate the production target of each child solution according to the production index corresponding to each child solution and the production data of each child solution, and the calculation formula of the production target, and sort all solutions in the current population and all solutions in the child solution population according to the production target of each solution in the current population and the constraint violation value cv of each solution, as well as the production target corresponding to each child solution in the child solution population and the constraint violation value cv corresponding to each child solution, obtain the top N solutions to generate a new population;

[0034] Step 3.6: Get the non-dominated solution in the previous iteration population of the current iteration number, calculate the average value of all production indicators in the non-dominated connection, and obtain the center point record PS t-1 , recalculate the new initial production index of 50% of N, replace the solution in the current population with the new initial production index, and obtain the first population;

[0035] Step 3.7: Determine the non-dominated solution of the first population and calculate the average value of each production index in the non-dominated solution to obtain the center point record PS t ;

[0036] Step 3.8: According to the constraint violation value cv of each solution in the first population, determine the feasible solutions and infeasible solutions in the first population, and update all solutions in the first population according to the feasible solution ratio rf in the first population, which is specifically implemented by the following formula:

[0037]

[0038] Where Δ represents the change in the center point of the non-dominated solution, that is, PS t and PS t-1 The difference between them, U(0, Δ) is an n-dimensional random vector that obeys a uniform distribution in the interval [0, Δ];

[0039] Step 3.9: Input the updated production index of each solution into the model of the relationship between the production index and the production data, calculate the production data of each solution, and then calculate the production target of each solution according to the production index of each solution, the production data of each solution, and the calculation formula of the production target. According to the production target of each solution and the constraint violation value cv of each solution, sort all the updated solutions and all the solutions in the first population, and obtain the top N solutions to generate a new population;

[0040] Step 4: Generate a new generation population based on the new population or all solutions in the current population. According to the new generation population and the current elite offspring, the current elite offspring contains multiple non-dominated solutions, determine the new elite offspring, and increase the current iteration number by one.

[0041] Step 5: Determine whether the number of iterations after the current number of iterations is increased by one is greater than the preset number of iterations. If the number of iterations after the current number of iterations is increased by one is not greater than the preset number of iterations, use the new generation population as the current population, use the new elite offspring as the current elite offspring, and return to execute step 2. If the number of iterations after the current number of iterations is increased by one is greater than the preset number of iterations, output all solutions in the new generation population, and use all solutions in the new generation population as control inputs for catalytic cracking production.

[0042] Optionally, the formula for the production constraint in step 1.1 includes:

[0043] 1x k ≤x k (t)≤ux k ;

[0044] ikB s,10％ ≤DT s,10％ (t)≤uDT s,10％ ;

[0045] ikB s,90％ ≤DT s,90％ (t)≤uDT s,90％ ;

[0046] oeLh s (t)≤Y s (t)≤uY s (t);

[0047]

[0048] Among them, lx k Indicates the minimum value of the kth data in the production index, ux k Indicates the maximum value of the kth data in the production index X(t), x k (t) represents the kth data in the production index; lDT s，10％Indicates the minimum value of the temperature at which 10% of the fuel oil is distilled, uDT s，10％ Indicates the maximum temperature at which 10% of the fuel oil is distilled, DT s，10％ (t) represents the temperature at which 10% of the fuel oil is distilled under the environment t; lDT s，90％ Indicates the minimum value of the temperature at which 90% of the fuel oil is distilled, uDT s，90％ Indicates the maximum temperature at which 90% of the fuel oil is distilled, DT s，90％ (t) represents the temperature at which 90% of the fuel oil is distilled under the environment t; lY s (t) represents the minimum output of the sth fuel oil under the t environment, uY s (t) represents the maximum output of the sth fuel oil under the t environment, Y s (t) represents the production of the sth fuel oil under environment t; Y1(t) represents the gasoline production under environment t, Y2(t) represents the diesel production under environment t, Y3(t) represents the natural gas production under environment t, YP(t) represents the pollutant production under environment t, and FHO(t) represents the wax oil flow rate under environment t.

[0049] Optionally, step 3.1 includes the following steps:

[0050] The formula of production constraints is rearranged and we get c m (X(t))<0 to obtain g new production constraints, where c m (X(t)) represents the mth production constraint condition. The constraint violation value cv is calculated according to the new production constraint condition. It is specifically expressed by the following formula:

[0051]

[0052] Among them, when cv(X(t)) is equal to 0, it means that X(t) is a feasible solution, which further determines the number of feasible solutions in the current population.

[0053] Optionally, the model of the relationship between the production indicators and the production data is constructed by the following steps:

[0054] Step A1: collecting production indicators and production data in the catalytic cracking unit;

[0055] Step A2: Construct a model of the relationship between production indicators and production data;

[0056] Step A2.1: Determine the input and output of the model of the relationship between production indicators and production data;

[0057] Among them, the production index is used as the model input and expressed by the following formula:

[0058] X=[FHO, THO, FC, TC, FLS, FSS, TR, HR1, HR2, HR3];

[0059] The production data is used as the model output and is expressed by the following formula:

[0060] O=[Y1,Y2,Y3,YP,DT s,10％ ,DT s,90％ ];

[0061] Step A2.2: Based on the input and output, a random weight neural network is used to establish a model of the relationship between production indicators and production data, which is specifically expressed by the following formula:

[0062]

[0063] Where ψ=[ψ1,ψ2,...,ψ P ] is the output weight, P is the number of hidden layer nodes, X represents the input variable, O(X) represents the output variable, W p represents the weight vector between the pth hidden layer neuron and the input, b p represents the threshold of the pth hidden layer neuron, G is the hidden layer activation function, and G can be expressed by the following formula:

[0064]

[0065] Among them, α and β are constants;

[0066] Step A2.3: Calculate the output weights and use the least squares method to find the output weights:

[0067] ψ=(H T H) -1 H T T;

[0068] Among them, H is the hidden layer output matrix, expressed as:

[0069]

[0070] T is the target matrix, expressed as:

[0071] T=[O1,O2,…O S ];

[0072] Where J is the number of samples collected, W1 represents the weight vector between the first hidden layer neuron and the input, b1 represents the threshold of the first hidden layer neuron, X1 and X J All are the input of the model, and O1, O2, and O3 are the output of the model.

[0073] Optionally, the calculation formula of the production target includes:

[0074]

[0075]

[0076] Among them, X(t) represents the production index under environment t, and X(t) is expressed by the following formula:

[0077] X(t)=[FC(t),TC(t),FSS(t),FLS(t),FHO(t),TR(t),THO(t),HR1(t),HR2(t),HR3(t)];

[0078] Wherein, FC(t) represents the catalyst flow rate under the environment of t, TC(t) represents the catalyst temperature under the environment of t, FSS(t) represents the stripping steam flow rate under the environment of t, FLS(t) represents the lifting steam flow rate under the environment of t, FHO(t) represents the wax oil flow rate under the environment of t, TR(t) represents the riser outlet temperature under the environment of t, THO(t) represents the wax oil temperature under the environment of t, HR1(t) represents the first circuit heat recovery under the environment of t, HR2(t) represents the second circuit heat recovery under the environment of t, HR3(t) represents the third circuit heat recovery under the environment of t;

[0079] Among them, PY s is PY1 or PY2 or PY3, PY1 represents the unit price of gasoline, PY2 represents the unit price of diesel, PY3 represents the unit price of natural gas, PYP represents the unit price of pollutant removal, PHO represents the unit price of wax oil, ρ is the energy required to increase 1 mol of catalyst by 1°C, es is the energy required to obtain 1 mol of steam; Y s (t) is Y1(t) or Y2(t) or Y3(t), Y1(t) represents the gasoline production under environment t, Y2(t) represents the diesel production under environment t, Y3(t) represents the natural gas production under environment t, and YP(t) represents the pollutant production under environment t;

[0080] Among them, Y1(t), Y2(t), Y3(t), and YP(t) are calculated through the model of the relationship between production indicators and production data, which are specifically expressed by the following formula:

[0081] Y1(t)=O1(X(t));

[0082] Y2(t)=O2(X(t));

[0083] Y3(t)=O3(X(t));

[0084] YP(t)=O p (X(t));

[0085] Among them, O1(X(t)) represents the first data in the matrix of production data corresponding to the production index under the t environment, O2(X(t)) represents the second data in the matrix of production data corresponding to the production index under the t environment, O3(X(t)) represents the third data in the matrix of production data corresponding to the production index under the t environment, and O p (X(t)) represents the fourth data in the matrix of production data corresponding to the production index under environment t.

[0086] Optionally, step 4 includes the following steps:

[0087] Step 4.1: Take the new population or all solutions in the current population as the second population. For all solutions in the second population, calculate the target value of the penalty item of each solution, and sort all solutions in the second population according to the target value of the penalty item of each solution and the constraint violation value cv corresponding to each solution. After the sorting is completed, the arrangement numbers of all solutions in the second population are from small to large. Randomly select two solutions from the N solutions in the second population. According to the arrangement numbers corresponding to the two solutions, select the solution with the smaller arrangement number as the parent solution. Repeat the above operation to obtain two parent solutions. Simulate binary crossover and polynomial mutation on the two parent solutions until N child solutions are obtained to generate a new child population.

[0088] Step 4.2: For all solutions in the new offspring population, calculate the target value of the penalty item of each solution, and sort all solutions in the second population and all solutions in the new offspring population according to the target value of the penalty item of each solution in the second population and the constraint violation value cv of each solution, as well as the target value of the penalty item of each offspring solution in the new offspring population and the constraint violation value cv corresponding to each offspring solution, obtain the top N solutions to generate a new generation population, compare the solutions in the new generation population with the non-dominated solutions in the elite offspring, determine N non-dominated solutions as the new elite offspring, and increase the current number of iterations by one.

[0089] Optionally, step 3.4 includes the following:

[0090] The simulated binary crossover of the parent solution S1 and the parent solution S2 is realized by the following formula:

[0091]

[0092]

[0093] Where Si1 is the i-th data in the production index corresponding to the parent solution S1, Si2 is the i-th data in the production index corresponding to the parent solution S2, ci1 is the i-th data in the production index corresponding to the child solution c1, ci2 is the i-th data in the production index corresponding to the child solution c2, r is a random number in the range [0, 1]; distribution index η = 5

[0094] Among them, the polynomial mutation of the parent solution S1 and the parent solution S2 is achieved by the following formula:

[0095]

[0096]

[0097] Among them, S i is the i-th data in the production index corresponding to the parent solution, that is, S i Solve S for the parent 1 The i-th data or parent solution S in the corresponding production index 2 The i-th data in the corresponding production index, is the i-th data in the production index corresponding to the sub-generation solution, that is, the sub-generation solution The i-th data and sub-generation solution in the corresponding production index The i-th data in the corresponding production index, λ is the distribution index, where the simulated binary crossover and polynomial mutation of the two parent solutions in step 4.1 are also calculated by the above formula.

[0098] The beneficial effects of adopting the above technical solution are:

[0099] The present invention can optimize the production indicators of multiple catalytic cracking processes according to the current constraints to achieve the goal of maximizing production profits and minimizing production energy consumption in accordance with user needs. At the same time, after the production constraints change, the optimization algorithm will respond dynamically, improve the optimization efficiency, and quickly obtain a set of optimal production indicators under the current environment. The present invention can significantly improve the production efficiency of the catalytic cracking process under a dynamic environment. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 A catalytic cracking production flow chart in an embodiment of the present invention.

[0101] Figure 2 A schematic flow chart of a dynamic constraint multi-objective optimization method for catalytic cracking production indicators in an embodiment of the present invention. DETAILED DESCRIPTION

[0102] The specific implementation of the present invention is further described in detail below in conjunction with the accompanying drawings and examples. The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

[0103] In order to solve the problem that the production constraints such as the fuel oil output of the catalytic cracking unit change over time, and the existing technology cannot quickly track the optimal setting value of the decision variable in a dynamic environment, the present invention proposes a constrained multi-objective optimization method for catalytic cracking production indicators in a dynamic environment. On the one hand, considering that the production constraints change, the change of the optimal solution set is irregular and difficult to predict. The present invention processes the population according to the number of feasible solutions under the new production constraints. When the number of feasible solutions under the new production constraints is large, the parent selection strategy is directly used to cross-mutate the population to generate offspring; when there are fewer feasible solutions, according to the moving direction of the non-dominated solution under the new constraints and the old constraints, the feasible solutions and infeasible solutions are moved to different degrees in the direction of the non-dominated solution. Then, the feasible solution priority principle is used to select N better individuals from the population and offspring. On the other hand, the user hopes that the comprehensive production indicators can be as many and good as possible within the constraints. The present invention uses an adaptive penalty function method to process the constraints, and uses the parent selection strategy to select the parent with better performance to obtain more feasible solutions close to the Pareto frontier. The above method can determine the production decision variables of relevant processes according to the current production status and output constraints, thereby improving production profits and reducing energy consumption in the production process.

[0104] The catalytic cracking production process consists of a series of parallel processes, each of which has its own production tasks. Figure 1 As shown in the figure, the process consists of gasification, lifting, separation, fractionation, etc., as follows:

[0105] After the vacuum wax oil from the atmospheric and vacuum unit is preheated by the heating furnace, it contacts the high-temperature catalyst from the regenerator at the bottom of the riser in the reactor, and the vacuum wax oil is instantly gasified to generate reaction oil and gas. The reaction oil and gas and catalyst particles rise along the riser under the action of the lifting steam. During the rising process, the reaction oil and gas undergo cracking reaction under the action of the catalyst to obtain hydrocarbons and coke. The coke will be deposited on the surface of the catalyst to make the catalyst inactive. When the hydrocarbons and the inactivated catalyst (catalyst to be regenerated) reach the top of the riser, the hydrocarbons and catalyst are separated under the action of the cyclone separator at the outlet of the riser. The hydrocarbons enter the fractionation device through the gas collecting chamber at the top of the settler for product separation, while the catalyst to be regenerated enters the stripping section of the reactor. The catalyst contacts the stripping steam in the stripping section in countercurrent to replace the oil and gas carried by the catalyst to be regenerated. The replaced catalyst to be regenerated enters the coke burning tank of the regenerator through the catalyst pipeline to be regenerated, and the catalyst is regenerated by burning the coke deposited on the surface of the catalyst to be regenerated. The regenerated catalyst (regenerated catalyst) and the flue gas generated during the combustion process move upward under the action of the main wind, and finally separate the catalyst and the flue gas through the separator at the top of the dilute phase tube. The separated catalyst enters the reactor through the regenerated catalyst pipeline to realize the recycling of the catalyst. The flue gas generated by charring is discharged through the smoke exhauster and other devices at the top of the regenerator. After entering the fractionation device, the hydrocarbon compound is introduced into the tower plate, and after passing through multiple tower plates and heat recovery loops, the hydrocarbon compound is fractionated according to different distillation temperatures. The heat recovered from the first loop and the second loop will cause diesel to be fractionated from the hydrocarbon compound, and the heat recovered from the third loop will cause gasoline and natural gas to be fractionated.

[0106] The present invention mainly solves the problem of optimizing the production index of each production process under the dynamic change of production constraints in the above catalytic cracking production process. The optimized production index of the catalytic cracking unit is issued to each process for production, so as to obtain the expected comprehensive production index.

[0107] The performance indicators involved in the catalytic cracking production process include the production indicators of each process, boundary constraints and comprehensive production indicators. Among them, the production indicators are wax oil flow (FHO), wax oil temperature (THO), catalyst flow (FC), catalyst temperature (TC), lifting steam flow (FLS), stripping steam flow (FSS), riser outlet temperature (TR), first circuit recovery heat (HR1), second circuit recovery heat (HR2), third circuit recovery heat (HR3). The comprehensive production indicator is the total production profit (f product ) and total production energy consumption (f energy ).

[0108] In the reactor of the catalytic cracking unit, when the temperature THO(t) of the wax oil (HO) and the catalyst temperature TC(t) increase, the energy consumption required by the unit increases, and the production cost also increases. At the same time, due to the increase in the temperature of the feed oil and the catalyst, the temperature of the riser increases, so the degree of cracking of HO increases. Conversely, when the temperature decreases, the degree of cracking of HO decreases. Generally speaking, when the degree of cracking of HO increases, more gasoline, natural gas, coke and pollutants are produced, while less diesel is produced. Conversely, when the degree of cracking of HO decreases, more diesel is produced but less gasoline and natural gas are produced. When the flow rate FC(t) of the catalyst increases, the reaction rate of cracking is usually increased, the output of gasoline and diesel is increased, and thus the production profit is increased. However, too much catalyst will also increase the content of pollutants in the product, resulting in increased decontamination costs and reduced profits. As a raw material, the flow rate FHO(t) of the wax oil directly affects the output of the final product. The appropriate flow rate helps to increase the output of ideal light hydrocarbon products (such as gasoline) while controlling the ratio of diesel and natural gas. However, too high a flow rate may cause catalyst overload, reduce reaction efficiency, increase pollutants, and reduce production profits. The flow rate FLS(t) of the lifting steam determines the upward driving force of the steam and affects the contact efficiency between the hot catalyst and the wax oil. The larger the FLS(t), the longer the contact time between the catalyst and HO, which can enhance the uniformity of the reaction, increase the yield of light hydrocarbons and medium-light fractions (such as gasoline and diesel), and increase profits. However, too large a flow rate may lead to unnecessary energy consumption. In the stripping section of the reactor, the increase in the flow rate FSS(t) of the stripping steam can more effectively remove the adsorbed light hydrocarbons from the catalyst surface and improve the regeneration effect of the catalyst. However, too much stripping steam may cause a sharp increase in energy consumption, which is not conducive to economic efficiency. Since the cracking reaction is an endothermic process, when the height of the riser increases, the temperature will decrease, thereby affecting the cracking reaction. We use the outlet temperature TR(t) of the riser to represent the intensity of the endothermic cracking inside the reactor. It is affected by other decision variables. Higher temperatures help to generate more light hydrocarbons, but may also increase the generation of pollutants.

[0109] In the regenerator of the catalytic cracking unit, the catalyst (HRC) needs to be heated by burning fuel. The higher the temperature TC(t) of the catalyst, the more energy is consumed.

[0110] In the distillation tower, the temperature of the trays is controlled by three heat recovery devices PA1, PA2 and PA3. The temperature of the trays will affect the production of fuel oil (gasoline, diesel and natural gas). HR1(t), HR2(t) and HR3(t) can recover the distilled energy and recycle the energy, effectively reducing energy consumption, but too much energy recovery will reduce the production of gasoline, diesel and natural gas.

[0111] Therefore, the production process of the catalytic cracking unit is very complex. In order to obtain better comprehensive production indicators, it is necessary to dynamically optimize multiple decision variables at the same time, that is, the production indicators of the catalytic cracking production process.

[0112] In view of the problems existing in the prior art, the present invention provides a dynamic constraint multi-objective optimization method for catalytic cracking production indicators. Figure 2 , which may include the following steps:

[0113] Step 1: Population initialization;

[0114] Step 1.1: Obtaining a formula for production constraints of a catalytic cracking unit, wherein the formula for production constraints includes a formula for representing a threshold interval of all data in a production index, a formula for representing a threshold interval of all data in production data, and a formula for representing a constraint relationship between the production index and the production data;

[0115] The production indicators include: wax oil flow FHO, wax oil temperature THO, catalyst flow FC, catalyst temperature TC, lifting steam flow FLS, stripping steam flow FSS, lifting pipe outlet temperature TR, first circuit recovery heat HR1, second circuit recovery heat HR2, third circuit recovery heat HR3; the production data include gasoline output Y1, diesel output Y2, natural gas output Y3, pollutant output YP, and 10% distillation temperature DT of gasoline, diesel and natural gas s，10％ The 90% distillation temperature DT of gasoline, diesel and natural gas s，90％ ;

[0116] Among them, the production constraints may include: catalytic cracking unit production capacity constraints, fuel oil distillation property constraints, fuel oil output constraints and material balance constraints, which can be specifically reflected by the formula of production constraints. The formula of production constraints includes:

[0117] 1x k ≤x k (t)≤ux k ;

[0118] ikB s,10％ ≤DT s,10％ (t)≤uDT s,10％ ;

[0119] ikB s,90％ ≤DT s,90％ (t)≤uDT s,90％ ;

[0120] oeLh s (t)≤Y s (t)≤uY s (t);

[0121]

[0122] Among them, lx k Indicates the minimum value of the kth data in the production index, ux k Indicates the maximum value of the kth data in the production index X(t), x k (t) represents the kth data in the production index; lDT s，10％ Indicates the minimum value of the temperature at which 10% of the fuel oil is distilled, uDT s，10％ Indicates the maximum temperature at which 10% of the fuel oil is distilled, DT s，10％ (t) represents the temperature at which 10% of the fuel oil is distilled under the environment t; lDT s，90％ Indicates the minimum value of the temperature at which 90% of the fuel oil is distilled, uDT s，90％ Indicates the maximum temperature at which 90% of the fuel oil is distilled, DT s，90％ (t) represents the temperature at which 90% of the fuel oil is distilled under the environment t; lY s (t) represents the minimum output of the sth fuel oil under the t environment, uY s (t) represents the maximum output of the sth fuel oil under the t environment, Y s (t) represents the production of the sth fuel oil under environment t; Y1(t) represents the gasoline production under environment t, Y2(t) represents the diesel production under environment t, Y3(t) represents the natural gas production under environment t, YP(t) represents the pollutant production under environment t, and FHO(t) represents the wax oil flow rate under environment t.

[0123] In the specific implementation, the production constraint conditions include the constraint boundaries of the production indicators and the production conditions, wherein the constraint boundaries of the production indicators are set as shown in Table 1 according to the data collected on site and the user's requirements for the production indicators;

[0124] Table 1 Constraint boundaries of production indicators

[0125] FHO THO(℃) FC TC(℃) FLS Minimum 230825 165.0229574 2020316.079 728.333668 5146.340281 Maximum 240242.66 185 2834843.914 769.2148593 8468.083144 FSS TR(℃) <![CDATA[HR1(J)]]> <![CDATA[HR2(J)]]> <![CDATA[HR3(J)]]> Minimum 3446.700817 550.2005184 27296232.82 26699053.59 18771882.95 Maximum 7034.60063 569.8786939 76650933.65 57153508.84 24616834.63

[0126] For production condition settings, this can include:

[0127] The unit price PYP for pollutant removal was set at 0.872;

[0128] The natural gas unit price PY3 is fixed at 4.3;

[0129] The upper and lower limits of 10% distillation temperature of fuel oil lDT s,10％ 、uDT s,10％ Set to -135℃ and -105℃;

[0130] The upper and lower limits of 90% distillation temperature of fuel oil lDT s,90％ 、uDT s,90％ Set to 270℃ and 300℃;

[0131] The settings of ten production examples of wax oil unit price, gasoline unit price, and diesel unit price are shown in Table 2;

[0132] Table 2 Ten production examples of wax oil unit price, gasoline unit price and diesel unit price

[0133]

[0134]

[0135] Step 1.2: Generate the initial production index based on the threshold range of all data in the production index in the production constraint conditions, which is calculated by the following formula:

[0136]

[0137] Among them, rand is a random number between [0, 1]. is the minimum value of the jth data in the production index, is the maximum value of the jth data in the production index;

[0138] N initial production indicators are calculated and used as N solutions. The N solutions constitute the initial population, and the initial number of iterations Gen is set to 1. In the specific implementation, the size of the population P is set to N = 200, and the maximum number of iterations gen in each dynamic environment is 10 and 20 respectively.

[0139] Step 2: Determine whether the production constraint conditions under the current iteration number environment are the same as the production constraint conditions of the previous iteration. If the production constraint conditions under the current iteration number environment are not the same as the production constraint conditions of the previous iteration, execute step 3. If the production constraint conditions under the current iteration number environment are the same as the production constraint conditions of the previous iteration, execute step 4. When the current iteration number is 1, execute step 4.

[0140] Among them, step 3 can be understood as a dynamic response process.

[0141] Step 3: Calculate the new population based on the current population;

[0142] Step 3.1: Generate a new constraint formula based on the production constraint formula, calculate the constraint violation value cv of each solution based on the new constraint formula, determine whether the solution of the current population is a feasible solution based on the constraint violation value cv of each solution, the constraint violation value cv of the feasible solution is 0, and obtain the number of feasible solutions in the current population;

[0143] Specifically, step 3.1 may include the following:

[0144] The formula of production constraints is rearranged and we get c m (X(t))<0 to obtain g new production constraints, where c m (X(t)) represents the mth production constraint condition. The constraint violation value cv is calculated according to the new production constraint condition. It is specifically expressed by the following formula:

[0145]

[0146] Among them, when cv(X(t)) is equal to 0, it means that X(t) is a feasible solution, which further determines the number of feasible solutions in the current population.

[0147] In the specific implementation, the formula of the production constraint condition is rearranged and the formula in the production constraint condition is: lx k ≤x k (t)≤ux k For example, after moving the items, we can get two new constraints:

[0148] 1x k -x k (t)≤0;

[0149] x k (t)-ux k ≤0;

[0150] Where c1(X(t))=lx k -x k (t), which is the first production constraint, c2(X(t))=x k (t)-ux k , which is the second production constraint. The formulas for other production constraints are the same.

[0151] Step 3.2: When the number of feasible solutions in the current population is greater than the preset threshold, execute steps 3.3 to 3.5. After step 3.5 is completed, execute step 4. When the number of feasible solutions in the current population is less than or equal to the preset threshold, execute steps 3.6 to 3.9. After step 3.9 is completed, execute step 4.

[0152] Step 3.3: Input the N production indicators corresponding to the N solutions in the current population into the model of the relationship between production indicators and production data to obtain N production data, and then calculate N production targets based on the N production indicators and N production data and the calculation formula of the production target. The production targets include maximizing the production profit f product and minimize production energy consumption f energy , calculate the maximum production profit f product The corresponding target value, and the minimization of production energy consumption f energy The corresponding target value is to maximize the production profit f product The corresponding target value and the minimized production energy consumption f energy The corresponding target values ​​constitute the target value of the penalty term;

[0153] The model of the relationship between the production indicators and the production data is constructed by the following steps:

[0154] Step A1: collecting production indicators and production data in the catalytic cracking unit;

[0155] Specifically, the collected data set includes 500 sets of production indicators and production data under different constraints and operating conditions.

[0156] Step A2: Construct a model of the relationship between production indicators and production data;

[0157] Step A2.1: Determine the input and output of the model of the relationship between production indicators and production data;

[0158] Among them, the production index is used as the model input and expressed by the following formula:

[0159] X=[FHO, THO, FC, TC, FLS, FSS, TR, HR1, HR2, HR3];

[0160] The production data is used as the model output and is expressed by the following formula:

[0161] O=[Y1,Y2,Y3,YP,DT s,10％ ,DT s,90％ ]

[0162] Step A2.2: Based on the input and output, a random weight neural network is used to establish a model of the relationship between production indicators and production data, which is specifically expressed by the following formula:

[0163]

[0164] Where ψ=[ψ1,ψ2,...,ψ P] is the output weight, P is the number of hidden layer nodes, X represents the input variable, O(X) represents the output variable, W p represents the weight vector between the pth hidden layer neuron and the input, b p represents the threshold of the pth hidden layer neuron, G is the hidden layer activation function, and G can be expressed by the following formula:

[0165]

[0166] Among them, α and β are constants, α = 1.342, β = 0.658;

[0167] Among them, the parameters of the random weight neural network are set as follows: the number of hidden layer nodes P = 50, the dimension of input data D = 10, and the input weight vector Wp and the hidden layer node threshold bp are randomly generated.

[0168] Step A2.3: Calculate the output weights and use the least squares method to find the output weights:

[0169] ψ=(H T H) -1 H T T;

[0170] Among them, H is the hidden layer output matrix, expressed as:

[0171]

[0172] T is the target matrix, expressed as:

[0173] T=[O1,O2,·O S ];

[0174] Where J is the number of samples collected, W1 represents the weight vector between the first hidden layer neuron and the input, b1 represents the threshold of the first hidden layer neuron, X1 and X J All are the input of the model, and O1, O2, and O3 are the output of the model.

[0175] The calculation formula of the production target includes:

[0176]

[0177]

[0178] Among them, X(t) represents the production index under environment t, and X(t) is expressed by the following formula:

[0179] X(t)=[FC(t),TC(t),FSS(t),FLS(t),FHO(t),TR(t),THO(t),HR1(t),HR2(t),HR3(t)]

[0180] Wherein, FC(t) represents the catalyst flow rate under the environment of t, TC(t) represents the catalyst temperature under the environment of t, FSS(t) represents the stripping steam flow rate under the environment of t, FLS(t) represents the lifting steam flow rate under the environment of t, FHO(t) represents the wax oil flow rate under the environment of t, TR(t) represents the riser outlet temperature under the environment of t, THO(t) represents the wax oil temperature under the environment of t, HR1(t) represents the first circuit heat recovery under the environment of t, HR2(t) represents the second circuit heat recovery under the environment of t, HR3(t) represents the third circuit heat recovery under the environment of t;

[0181] Among them, PY s is PY1 or PY2 or PY3, PY1 represents the unit price of gasoline, PY2 represents the unit price of diesel, PY3 represents the unit price of natural gas, PYP represents the unit price of pollutant removal, PHO represents the unit price of wax oil, ρ is the energy required to increase the temperature of 1 mol of catalyst by 1°C, es is the energy required to obtain 1 mol of steam, in specific implementation, ρ=1.1, es=2260; Y s (t) is Y1(t) or Y2(t) or Y3(t), Y1(t) represents the gasoline production under environment t, Y2(t) represents the diesel production under environment t, Y3(t) represents the natural gas production under environment t, and YP(t) represents the pollutant production under environment t;

[0182] Among them, Y1(t), Y2(t), Y3(t), and YP(t) are calculated through the model of the relationship between production indicators and production data, which are specifically expressed by the following formula:

[0183] Y1(t)=O1(X(t))

[0184] Y2(t)=O2(X(t))

[0185] Y3(t)=O3(X(t))

[0186] YP(t)=O p (X(t))

[0187] Among them, O1(X(t)) represents the first data in the matrix of production data corresponding to the production index under the t environment, O2(X(t)) represents the second data in the matrix of production data corresponding to the production index under the t environment, O3(X(t)) represents the third data in the matrix of production data corresponding to the production index under the t environment, and O p(X(t)) represents the fourth data in the matrix of production data corresponding to the production index under environment t.

[0188] In specific implementation, the following formula can also be used:

[0189] DT s,10％ (t) = O s,10％ (X(t));

[0190] DT s,90％ (t) = O s,90％ (X(t));

[0191] Among them, O s,10％ (X(t)) represents the fifth data in the matrix of production data corresponding to the production index under the t environment, O s,90％ (X(t)) represents the sixth data in the matrix of production data corresponding to the production index under the t environment

[0192] Among them, minimize the production energy consumption f energy The corresponding target value is calculated by the following formula:

[0193]

[0194]

[0195] Where cv' is the normalized constraint violation value corresponding to the jth X(t), r f is the proportion of feasible solutions in the current population, f j '(X(t)) is the minimum production energy consumption f corresponding to the jth X(t) energy The normalized value, f j 1 (X(t)) is the minimum production energy consumption f energy The corresponding target value; X(t) is feasible, which means X(t) is a feasible solution; X(t) is infeasible, which means X(t) is an infeasible solution;

[0196] Among them, maximize the production profit f product The corresponding target value is calculated by the following formula:

[0197]

[0198]

[0199] Among them, f j ”(X(t)) is the maximum production profit f corresponding to the jth X(t) product The normalized value, f j 2(X(t)) is the maximum production profit f corresponding to the jth X(t) product Normalized value;

[0200] Thus, the target values ​​of the N penalty items are obtained. According to the target values ​​of the N penalty items, the N solutions in the current population are sorted. After the sorting is completed, the arrangement sequence of the N solutions in the current population is from small to large. In the specific implementation, the solutions in the current population can be sorted by the NSGA-II method.

[0201] Step 3.4: Randomly select two solutions from the N solutions in the current population, and select the solution with the smaller sequence number as the parent solution S according to the sequence numbers corresponding to the two solutions. 1 , repeat the above operation to get the parent solution S 2 , solve S for the parent 1 and the parent solution S 2 Perform simulated binary crossover to obtain the offspring solution c 1 and the offspring solution c 2 , solve S for the parent 1 and the parent solution S 2 Perform polynomial mutation to obtain offspring solution and offspring solution In the offspring solution c 1 、Offspring solution c 2 、Sub-generation solution and offspring solution When the number reaches N, the offspring population is generated according to the N offspring solutions;

[0202] Among them, select the parent solution S 1 and the parent solution S 2 This can be achieved through the parent selection operator of dCMOEA; the crossover probability and distribution index of the simulated binary crossover are set to 0.9 and 5, respectively, and the mutation probability and distribution index of the polynomial mutation are set to 0.1 and 20, respectively.

[0203] Among them, step 3.4 specifically includes the following contents:

[0204] Solve S for the parent 1 and the parent solution S 2 The simulated binary crossover is implemented by the following formula:

[0205]

[0206]

[0207] Among them, S i 1 Solve S for the parent 1 The i-th data in the corresponding production index, S i2 To show the parent generation solution S 2 The i-th data in the corresponding production index, c i 1 Solve c for the child 1 The i-th data in the corresponding production index, c i 2 Represents the sub-generation solution c 2 The i-th data in the corresponding production index, r is a random number in the range [0, 1], and the distribution index η = 5;

[0208] Among them, the solution S for the parent 1 and the parent solution S 2 Polynomial mutation is achieved by the following formula:

[0209]

[0210]

[0211] Among them, S i is the i-th data in the production index corresponding to the parent solution, that is, S i Solve S for the parent 1 The i-th data or parent solution S in the corresponding production index 2 The i-th data in the corresponding production index, is the i-th data in the production index corresponding to the sub-generation solution, that is, the sub-generation solution The i-th data and sub-generation solution in the corresponding production index The corresponding production index is the i-th data, λ is the distribution index, and λ=20 in the specific implementation.

[0212] Step 3.5: Input the production index corresponding to each child solution in the child population into the model of the relationship between the production index and the production data, calculate the production data of each child solution, calculate the production target of each child solution according to the production index corresponding to each child solution and the production data of each child solution, and the calculation formula of the production target, and sort all solutions in the current population and all solutions in the child solution population according to the production target of each solution in the current population and the constraint violation value cv of each solution, as well as the production target corresponding to each child solution in the child solution population and the constraint violation value cv corresponding to each child solution, and obtain the top N solutions to generate a new population.

[0213] Among them, all solutions in the current population and all solutions in the child solution population can be sorted by the NSGA-II method, the production data of each child solution can be calculated according to the formula of the model of the relationship between the production index and the production data in step 3.3, the production target of each child solution can be calculated according to the calculation formula of the production target in step 3.3, and the constraint violation value cv of each solution can be obtained according to the formula of the constraint violation value cv in step 3.1.

[0214] Step 3.6: Get the non-dominated solution in the previous iteration population of the current iteration number, calculate the average value of all production indicators in the non-dominated connection, and obtain the center point record PS t-1 , recalculate the new initial production index of 50% of N, replace the solution in the current population with the new initial production index, and obtain the first population;

[0215] It should be noted that when replacing the solution in the current population with a new initial production indicator, it is necessary to ensure that there is at least one feasible solution in the population.

[0216] Step 3.7: Determine the non-dominated solution of the first population and calculate the average value of each production index in the non-dominated solution to obtain the center point record PS t ;

[0217] Step 3.8: According to the constraint violation value cv of each solution in the first population, determine the feasible solutions and infeasible solutions in the first population, and update all solutions in the first population according to the feasible solution ratio rf in the first population. The updating process can be understood as the moving direction of the center point of the non-dominated solution to make certain adjustments to the feasible solutions and infeasible solutions respectively, so as to realize the process of moving the feasible solutions and infeasible solutions to the non-dominated solutions, which is specifically realized by the following formula:

[0218]

[0219] Where Δ represents the change in the center point of the non-dominated solution, that is, PS t and PS t-1 The difference between them, U(0, Δ) is an n-dimensional random vector that obeys a uniform distribution in the interval [0, Δ];

[0220] Step 3.9: Input the updated production index of each solution into the model of the relationship between the production index and the production data, calculate the production data of each solution, and then calculate the production target of each solution according to the production index of each solution, the production data of each solution, and the calculation formula of the production target. According to the production target of each solution and the constraint violation value cv of each solution, sort all the updated solutions and all the solutions in the first population, and obtain the top N solutions to generate a new population;

[0221] Step 4: Generate a new generation population based on the new population or all solutions in the current population. According to the new generation population and the current elite offspring, the current elite offspring contains multiple non-dominated solutions, determine the new elite offspring, and increase the current iteration number by one.

[0222] Step 4.1: Take the new population or all solutions in the current population as the second population. For all solutions in the second population, calculate the target value of the penalty item of each solution, and sort all solutions in the second population according to the target value of the penalty item of each solution and the constraint violation value cv corresponding to each solution. After the sorting is completed, the arrangement numbers of all solutions in the second population are from small to large. Randomly select two solutions from the N solutions in the second population. According to the arrangement numbers corresponding to the two solutions, select the solution with the smaller arrangement number as the parent solution. Repeat the above operation to obtain two parent solutions. Simulate binary crossover and polynomial mutation on the two parent solutions until N child solutions are obtained to generate a new child population.

[0223] Among them, for all solutions in the second population, the target value of the penalty item of each solution is calculated. Specifically, the production data of each solution is calculated. The specific calculation refers to the formula of the model of the relationship between the production index and the production data in step 3.3. Then, according to the production index and the production data, the production target is calculated. The specific calculation refers to the calculation formula of the production target in step 3.3. Then, the target value of the penalty item corresponding to the production target is calculated. The specific calculation refers to the minimization of production energy consumption f in step 3.3. energy The corresponding target value calculation formula and the maximum production profit f product The calculation formula of the corresponding target value.

[0224] Among them, sorting all solutions in the second population can be achieved by the NSGA-II method, the two parent solutions can be obtained by the parent selection operator of dCMOEA, and the two parent solutions are simulated binary crossover and polynomial mutation. Please refer to the formula in step 3.4.

[0225] Step 4.2: For all solutions in the new offspring population, calculate the target value of the penalty item of each solution, and sort all solutions in the second population and all solutions in the new offspring population according to the target value of the penalty item of each solution in the second population and the constraint violation value cv of each solution, as well as the target value of the penalty item of each offspring solution in the new offspring population and the constraint violation value cv corresponding to each offspring solution, obtain the top N solutions to generate a new generation population, compare the solutions in the new generation population with the non-dominated solutions in the elite offspring, determine N non-dominated solutions as the new elite offspring, and increase the current number of iterations by one.

[0226] Among them, for all solutions in the new offspring population, the target value of the penalty item of each solution is calculated. The specific process refers to the calculation process of the target value of the penalty item of each solution for all solutions in the second population.

[0227] It should be noted that the elite offspring with iteration number 1 is the new generation population with iteration number 1.

[0228] Step 5: Determine whether the number of iterations after the current number of iterations is increased by one is greater than the preset number of iterations. If the number of iterations after the current number of iterations is increased by one is not greater than the preset number of iterations, use the new generation population as the current population, use the new elite offspring as the current elite offspring, and return to execute step 2. If the number of iterations after the current number of iterations is increased by one is greater than the preset number of iterations, output all solutions in the new generation population, and use all solutions in the new generation population as control inputs for catalytic cracking production.

[0229] Among them, the optimized catalytic cracking production indicators will be distributed to each process through the industrial control system.

[0230] The optimized values ​​of decision variables such as temperature and flow rate are used as control inputs for catalytic cracking production.

[0231] In the above example, we used random weight neural network to establish the relationship between catalytic cracking production index and fuel oil output, used constraint processing method to deal with infeasible solutions, and used population feasibility to dynamically respond to changes in production constraints, thereby improving the optimization efficiency of the algorithm.

[0232] In order to verify that the proposed dynamic constraint multi-objective optimization method for catalytic cracking production indicators has an advanced level, the hypervolume (HV) index is used as the evaluation index. The larger the HV, the better the optimization effect of the algorithm.

[0233]

[0234] Where VOL is the Lebesgue measure used to measure volume, and Z is the set of feasible non-dominated solutions obtained by optimization. i represents the hypervolume parameter consisting of the i-th solution and the reference point in the set Z. According to our preliminary experiments, the two target values ​​are in the hundreds of thousands. In order to calculate the HV of the non-dominated solution, this paper proposes to product and f energy The objective values ​​of are divided by 100000 and 100000000 respectively, and (0,10) is used as the reference point for calculating the non-dominated solution HV.

[0235] The proposed method is compared with other proposed dynamic constraint multi-objective optimization algorithms dCMOEA, DC-MOEA, DC-NSGA-II-A, and TDCEA. In the experiment, we selected 10 test instances, in which nt=20 is set in each test instance, that is, there are 20 different production constraints, and the change frequency τt of each production constraint is set to 10 generations and 20 generations. At the same time, in order to reduce the impact of static optimization, the algorithm runs 40 times under the first production constraint, and the above experiment is repeated 30 times to obtain the average value and variance of HV, where the values ​​outside the brackets are the average values ​​of the 30 experimental results, and the values ​​inside the brackets are the variances of the 30 experimental results. The specific experimental results are shown in Table 3.

[0236] Table 3 Experimental results

[0237]

[0238] It can be seen from Table 3 that the method of the present invention has an advanced level in solving the dynamic constraint multi-objective optimization problem of catalytic cracking production indicators.

[0239] The above description is only a preferred embodiment of the present disclosure and an explanation of the technical principles used. Those skilled in the art should understand that the scope of the invention involved in the embodiments of the present disclosure is not limited to the technical solutions formed by a specific combination of the above-mentioned technical features, but should also cover other technical solutions formed by any combination of the above-mentioned technical features or their equivalent features without departing from the above-mentioned inventive concept. For example, the above-mentioned features are replaced with the technical features with similar functions disclosed in the embodiments of the present disclosure (but not limited to) to form a technical solution.

Claims

1. A dynamic constrained multi-objective optimization method for catalytic cracking production indicators, characterized in that: include: Step 1: Population initialization; Step 1.1: Obtaining a formula for production constraints of a catalytic cracking unit, wherein the formula for production constraints includes a formula for representing a threshold interval of all data in a production index, a formula for representing a threshold interval of all data in production data, and a formula for representing a constraint relationship between the production index and the production data; The production indicators include: wax oil flow FHO, wax oil temperature THO, catalyst flow FC, catalyst temperature TC, lifting steam flow FLS, stripping steam flow FSS, lifting pipe outlet temperature TR, first circuit recovery heat HR1, second circuit recovery heat HR2, third circuit recovery heat HR3; the production data include gasoline output Y1, diesel output Y2, natural gas output Y3, pollutant output YP, and 10% distillation temperature DT of gasoline, diesel and natural gas s，10％ The 90% distillation temperature DT of gasoline, diesel and natural gas s，90％ ; Step 1.2: Generate the initial production index based on the threshold range of all data in the production index in the production constraint conditions, which is calculated by the following formula: Among them, rand is a random number between [0, 1]. is the minimum value of the jth data in the production index, is the maximum value of the jth data in the production index; Calculate and obtain N initial production indicators, use the N initial production indicators as N solutions, and the N solutions form an initial population. Set the initial iteration number Gen=1; Step 2: Determine whether the production constraint conditions under the current iteration number environment are the same as the production constraint conditions of the previous iteration. If the production constraint conditions under the current iteration number environment are not the same as the production constraint conditions of the previous iteration, execute step 3. If the production constraint conditions under the current iteration number environment are the same as the production constraint conditions of the previous iteration, execute step 4. When the current iteration number is 1, execute step 4. Step 3: Calculate the new population based on the current population; Step 3.1: Generate a new constraint formula based on the production constraint formula, calculate the constraint violation value cv of each solution based on the new constraint formula, determine whether the solution of the current population is a feasible solution based on the constraint violation value cv of each solution, the constraint violation value cv of the feasible solution is 0, and obtain the number of feasible solutions in the current population; Step 3.2: When the number of feasible solutions in the current population is greater than the preset threshold, execute steps 3.3 to 3.

5. After step 3.5 is completed, execute step 4. When the number of feasible solutions in the current population is less than or equal to the preset threshold, execute steps 3.6 to 3.

9. After step 3.9 is completed, execute step 4. Step 3.3: Input the N production indicators corresponding to the N solutions in the current population into the model of the relationship between production indicators and production data to obtain N production data, and then calculate N production targets based on the N production indicators and N production data and the calculation formula of the production target. The production targets include maximizing the production profit f product and minimize production energy consumption f energy , calculate the maximum production profit f product The corresponding target value, and the minimization of production energy consumption f energy The corresponding target value is to maximize the production profit f product The corresponding target value and the minimized production energy consumption f energy The corresponding target values ​​constitute the target value of the penalty term; Among them, minimize the production energy consumption f energy The corresponding target value is calculated by the following formula: Where cv' is the normalized constraint violation value corresponding to the jth X(t), r f is the proportion of feasible solutions in the current population, f j '(X(t)) is the minimum production energy consumption f corresponding to the jth X(t) energy The normalized value, f j 1 (X(t)) is the minimum production energy consumption f energy The corresponding target value; X(t) is feasible, which means X(t) is a feasible solution; X(t) is infeasible, which means X(t) is an infeasible solution; Among them, maximize the production profit f product The corresponding target value is calculated by the following formula: Among them, f″ j (X(t)) is the maximum production profit f corresponding to the jth X(t) product The normalized value, f j 2 (X(t)) is the maximum production profit f corresponding to the jth X(t) product Normalized value; Thus, the target values ​​of the N penalty items are obtained. According to the target values ​​of the N penalty items, the N solutions in the current population are sorted. After the sorting is completed, the N solutions in the current population are arranged in ascending order. Step 3.4: Randomly select two solutions from the N solutions in the current population, and select the solution with the smaller sequence number as the parent solution S according to the sequence numbers corresponding to the two solutions. 1 , repeat the above operation to get the parent solution S 2 , solve S for the parent 1 and the parent solution S 2 Perform simulated binary crossover to obtain the offspring solution c 1 and the offspring solution c 2 , solve S for the parent 1 and the parent solution S 2 Perform polynomial mutation to obtain offspring solution and offspring solution In the offspring solution c 1 、Offspring solution c 2 、Sub-generation solution and offspring solution When the number reaches N, the offspring population is generated according to the N offspring solutions; Step 3.5: Input the production index corresponding to each child solution in the child population into the model of the relationship between the production index and the production data, calculate the production data of each child solution, calculate the production target of each child solution according to the production index corresponding to each child solution and the production data of each child solution, and the calculation formula of the production target, and sort all solutions in the current population and all solutions in the child solution population according to the production target of each solution in the current population and the constraint violation value cv of each solution, as well as the production target corresponding to each child solution in the child solution population and the constraint violation value cv corresponding to each child solution, obtain the top N solutions to generate a new population; Step 3.6: Get the non-dominated solution in the previous iteration population of the current iteration number, calculate the average value of all production indicators in the non-dominated connection, and obtain the center point record PS t-1 , recalculate the new initial production index of 50% of N, replace the solution in the current population with the new initial production index, and obtain the first population; Step 3.7: Determine the non-dominated solution of the first population and calculate the average value of each production index in the non-dominated solution to obtain the center point record PS t ; Step 3.8: According to the constraint violation value cv of each solution in the first population, determine the feasible solutions and infeasible solutions in the first population, and update all solutions in the first population according to the feasible solution ratio rf in the first population, which is specifically implemented by the following formula: Where Δ represents the change in the center point of the non-dominated solution, that is, PS t and PS t-1 The difference between them, U(0, Δ) is an n-dimensional random vector that obeys a uniform distribution in the interval [0, Δ]; Step 3.9: Input the updated production index of each solution into the model of the relationship between the production index and the production data, calculate the production data of each solution, and then calculate the production target of each solution according to the production index of each solution, the production data of each solution, and the calculation formula of the production target. According to the production target of each solution and the constraint violation value cv of each solution, sort all the updated solutions and all the solutions in the first population, and obtain the top N solutions to generate a new population; Step 4: Generate a new generation population based on the new population or all solutions in the current population. According to the new generation population and the current elite offspring, the current elite offspring contains multiple non-dominated solutions, determine the new elite offspring, and increase the current iteration number by one. Step 5: Determine whether the number of iterations after the current number of iterations is increased by one is greater than the preset number of iterations. If the number of iterations after the current number of iterations is increased by one is not greater than the preset number of iterations, use the new generation population as the current population, use the new elite offspring as the current elite offspring, and return to execute step 2. If the number of iterations after the current number of iterations is increased by one is greater than the preset number of iterations, output all solutions in the new generation population, and use all solutions in the new generation population as control inputs for catalytic cracking production.

2. The dynamic constraint multi-objective optimization method for catalytic cracking production indicators according to claim 1, characterized in that: The formula for the production constraints in step 1.1 includes: lx k ≤x k (t)≤ux k ; lDT s,10％ ≤DT s,10％ (t)≤uDT s,10％ ; lDT s,90％ ≤DT s,90％ (t)≤uDT s,90％ ; glass s (t)≤Y s (t)≤uY s (t); Among them, lx k Indicates the minimum value of the kth data in the production index, ux k Indicates the maximum value of the kth data in the production index X(t), x k (t) represents the kth data in the production index; lDT s，10％ Indicates the minimum value of the temperature at which 10% of the fuel oil is distilled, uDT s，10％ Indicates the maximum temperature at which 10% of the fuel oil is distilled, DT s，10％ (t) represents the temperature at which 10% of the fuel oil is distilled under the environment t; lDT s，90％ Indicates the minimum value of the temperature at which 90% of the fuel oil is distilled, uDT s，90％ Indicates the maximum temperature at which 90% of the fuel oil is distilled, DT s，90％ (t) represents the temperature at which 90% of the fuel oil is distilled under the environment t; lY s (t) represents the minimum output of the sth fuel oil under the t environment, uY s (t) represents the maximum output of the sth fuel oil under the t environment, Y s (t) represents the production of the sth fuel oil under environment t; Y1(t) represents the gasoline production under environment t, Y2(t) represents the diesel production under environment t, Y3(t) represents the natural gas production under environment t, YP(t) represents the pollutant production under environment t, and FHO(t) represents the wax oil flow rate under environment t.

3. The dynamic constraint multi-objective optimization method for catalytic cracking production indicators according to claim 1, characterized in that: Step 3.1 includes the following steps: The formula of production constraints is rearranged and we get c m (X(t))<0 to obtain g new production constraints, where c m (X(t)) represents the mth production constraint condition. The constraint violation value cv is calculated according to the new production constraint condition. It is specifically expressed by the following formula: Among them, when cv(X(t)) is equal to 0, it means that X(t) is a feasible solution, which further determines the number of feasible solutions in the current population.

4. The dynamic constraint multi-objective optimization method for catalytic cracking production indicators according to claim 1, characterized in that: The model of the relationship between the production indicators and the production data is constructed by the following steps: Step A1: collecting production indicators and production data in the catalytic cracking unit; Step A2: Construct a model of the relationship between production indicators and production data; Step A2.1: Determine the input and output of the model of the relationship between production indicators and production data; Among them, the production index is used as the model input and expressed by the following formula: X=[FHO, THO, FC, TC, FLS, FSS, TR, HR1, HR2, HR3]; The production data is used as the model output and is expressed by the following formula: O=[Y1,Y2,Y3,YP,DT s,10％ ,DT s,90％ ]; Step A2.2: Based on the input and output, a random weight neural network is used to establish a model of the relationship between production indicators and production data, which is specifically expressed by the following formula: Where ψ=[ψ1,ψ2,...,ψ P ] is the output weight, P is the number of hidden layer nodes, X represents the input variable, O(X) represents the output variable, W p represents the weight vector between the pth hidden layer neuron and the input, b p represents the threshold of the pth hidden layer neuron, G is the hidden layer activation function, and G can be expressed by the following formula: Among them, α and β are constants; Step A2.3: Calculate the output weights and use the least squares method to find the output weights: ψ=(H T H) -1 H T T; Among them, H is the hidden layer output matrix, expressed as: T is the target matrix, expressed as: T=[O1,O2,…O S ]; Where J is the number of samples collected, W1 represents the weight vector between the first hidden layer neuron and the input, b1 represents the threshold of the first hidden layer neuron, X1 and X J All are the input of the model, and O1, O2, and O3 are the output of the model.

5. The dynamic constraint multi-objective optimization method for catalytic cracking production indicators according to claim 3, characterized in that: The calculation formula of the production target includes: Among them, X(t) represents the production index under environment t, and X(t) is expressed by the following formula: X(t)=[FC(t),TC(t),FSS(t),FLS(t),FHO(t),TR(t),THO(t),HR1(t),HR2(t),HR3(t)]; Wherein, FC(t) represents the catalyst flow rate under the environment of t, TC(t) represents the catalyst temperature under the environment of t, FSS(t) represents the stripping steam flow rate under the environment of t, FLS(t) represents the lifting steam flow rate under the environment of t, FHO(t) represents the wax oil flow rate under the environment of t, TR(t) represents the riser outlet temperature under the environment of t, THO(t) represents the wax oil temperature under the environment of t, HR1(t) represents the first circuit heat recovery under the environment of t, HR2(t) represents the second circuit heat recovery under the environment of t, HR3(t) represents the third circuit heat recovery under the environment of t; Among them, PY s is PY1 or PY2 or PY3, PY1 represents the unit price of gasoline, PY2 represents the unit price of diesel, PY3 represents the unit price of natural gas, PYP represents the unit price of pollutant removal, PHO represents the unit price of wax oil, ρ is the energy required to increase 1 mol of catalyst by 1°C, es is the energy required to obtain 1 mol of steam; Y s (t) is Y1(t) or Y2(t) or Y3(t), Y1(t) represents the gasoline production under environment t, Y2(t) represents the diesel production under environment t, Y3(t) represents the natural gas production under environment t, and YP(t) represents the pollutant production under environment t; Among them, Y1(t), Y2(t), Y3(t), and YP(t) are calculated through the model of the relationship between production indicators and production data, which are specifically expressed by the following formula: Y1(t)=O1(X(t)); Y2(t)=O2(X(t)); Y3(t)=O3(X(t)); YP(t)=O p (X(t)); Among them, O1(X(t)) represents the first data in the matrix of production data corresponding to the production index under the t environment, O2(X(t)) represents the second data in the matrix of production data corresponding to the production index under the t environment, O3(X(t)) represents the third data in the matrix of production data corresponding to the production index under the t environment, and O p (X(t)) represents the fourth data in the matrix of production data corresponding to the production index under environment t.

6. The dynamic constraint multi-objective optimization method for catalytic cracking production indicators according to claim 1, characterized in that: Step 4 includes the following steps: Step 4.1: Take the new population or all solutions in the current population as the second population. For all solutions in the second population, calculate the target value of the penalty item of each solution, and sort all solutions in the second population according to the target value of the penalty item of each solution and the constraint violation value cv corresponding to each solution. After the sorting is completed, the arrangement numbers of all solutions in the second population are from small to large. Randomly select two solutions from the N solutions in the second population. According to the arrangement numbers corresponding to the two solutions, select the solution with the smaller arrangement number as the parent solution. Repeat the above operation to obtain two parent solutions. Simulate binary crossover and polynomial mutation on the two parent solutions until N child solutions are obtained to generate a new child population. Step 4.2: For all solutions in the new offspring population, calculate the target value of the penalty item of each solution, and sort all solutions in the second population and all solutions in the new offspring population according to the target value of the penalty item of each solution in the second population and the constraint violation value cv of each solution, as well as the target value of the penalty item of each offspring solution in the new offspring population and the constraint violation value cv corresponding to each offspring solution, obtain the top N solutions to generate a new generation population, compare the solutions in the new generation population with the non-dominated solutions in the elite offspring, determine N non-dominated solutions as the new elite offspring, and increase the current number of iterations by one.

7. The dynamic constraint multi-objective optimization method for catalytic cracking production indicators according to claim 6, characterized in that: Step 3.4 includes the following: Solve S for the parent 1 and the parent solution S 2 The simulated binary crossover is implemented by the following formula: Among them, S i 1 Solve S for the parent 1 The i-th data in the corresponding production index, S i 2 To show the parent generation solution S 2 The i-th data in the corresponding production index, c i 1 Solve c for the child 1 The i-th data in the corresponding production index, c i 2 Represents the sub-generation solution c 2 The i-th data in the corresponding production index, r is a random number in the range [0, 1]; the distribution index η = 5 Among them, the solution S for the parent 1 and the parent solution S 2 Polynomial mutation is achieved by the following formula: Among them, S i is the i-th data in the production index corresponding to the parent solution, that is, S i Solve S for the parent 1 The i-th data or parent solution S in the corresponding production index 2 The i-th data in the corresponding production index, is the i-th data in the production index corresponding to the sub-generation solution, that is, the sub-generation solution The i-th data and sub-generation solution in the corresponding production index The i-th data in the corresponding production index, λ is the distribution index, where the simulated binary crossover and polynomial mutation of the two parent solutions in step 4.1 are also calculated by the above formula.

Images