Nonlinear modeling and interior point algorithm for the material flow optimization in petroleum refinery

: This paper established a mathematical model with nonconvex bilinear terms. It formulated the complex material flow in the petroleum refinery scenario based on the concept of the “P model”. The mathematical model described the nonlinear constraints such as linear and nonlinear mass and volume intersection flow blending of crude and middle material physical properties. Additionally, it described the complex inflow and outflow in secondary devices as nonlinear constraints such as delta-base structure and physical property transfer. It is highly difficult to determine the direction and quantity of each material in the network of refineries. An improved interior point algorithm with an initial point strategy was proposed to find a high-quality feasible solution in a short time. The real instances from the petroleum refinery were employed to compare and analyze the solutions from the improved algorithm and commercial solver. The experimental results show that the proposed algorithm framework can balance the solution quality and computational efficiency and perform well in different scenarios of refinery material flow networks.


Introduction
The pooling problem, also known as the mixed-flow problem, is a classic revenue maximization problem involving supply, storage, demand, and specified flow mixture constraints and usually appears in petroleum refinery and wastewater treatment scenarios [1].Crude oils flowing into the petroleum refinery network are refined and processed to produce gasoline, kerosene, diesel, and other liquid fuels, hundreds of kinds of lubricating oil, paraffin, asphalt, and other refinery products, and tens of thousands of synthetic resin, synthetic rubber, synthetic fiber, and other chemical products.Various physical and chemical processes, including crude oil separation, light oil, and gas processing, should be performed for the crude oil processing of various products.Each of these processes is relatively independent in constructing a refinery production unit network, such as the atmospheric-vacuum distillation unit and catalytic cracking, catalytic reforming, hydrocracking, coking, and other secondary devices.Refinery optimization refers to the overall optimization of the whole process of refinery enterprises, including raw material procurement, production, processing, and product delivery, for maximizing the economic benefits of enterprises while satisfying various production constraints.In terms of the pooling problem, the mixture of quality attributes often leads to nonconvex bilinear constraints, causing significant decision-making difficulties.Therefore, the research mainly focuses on the formulation and solution of nonlinear models.
There are three classical modeling methods for the pooling problem: the "P(proportional) model" proposed by Haverly [2], the "Q(quantity) model" improved by Ben-Tal [3], and the "PQ(proportional and quatity) model" further modified by Quesada and Grossmann [4].The "Q model" takes actual material inflow and product outflow of the processing system as the decision variables.The "P model" replaces the amount of material inflow and outflow by the transformation proportion relationship between all processes to obtain the equivalent form of the "Q model" with fewer variables.Based on the "Q model", the "PQ model" adds constraints of proportional variables and product output to obtain a tighter convex envelope relaxation model [5].
The global optimization algorithms for solving the pooling problem in academia and industry are listed as follows: the nonlinear programming method based on Lagrange relaxation proposed by Ben-Tal et al. [3,6], the branch-and-bound method based on convex relaxation proposed by Zamora and Grossmann [7], the branch-and-bound method for convex relaxation with enhanced lower bound based on discretization proposed by Ruiz and Grossmann [8], reformulation-linearization-technique(RLT) nonlinear method proposed by Meyer and Floudas [9], and external approximation nonlinear method proposed by Bergamini [10].Although these global optimization algorithms can theoretically guarantee the optimality of the results, their solution's efficiency cannot meet the practical operation requirements for complex problems.
Compared with the global optimization algorithms, several local optimization algorithms have been proposed to solve the pooling problem.These algorithms include the continuous linear programming solution (SLP) designed for refinery scenarios [11][12][13][14], the nonmonotonic continuous linear programming (NMSLP) method [15], the nonlinear term discretization method [16][17][18], the Benders decomposition method proposed by Floudas and Aggarwal [19], the distributed recursive method for solving pooling problems in commercial software [20][21][22], and an interior-point method for solving nonlinear problems in commercial software [23,24].Among them, the Benders decomposition and interior point methods fail to converge even for a long time.The continuous linear programming and distributed recursive methods easily fall into the local optimum.
This paper proposes a formulation of a nonlinear model suitable for refinery mixed flow based on the modeling of the "P model" in the pooling problem combined with complex real business scenarios.Furthermore, an improved interior point method algorithm is presented, which can find a high-quality feasible solution quickly.This algorithm has contributed to the generalized mixed-flow processing problems considering the complex nonlinear constraints in the refinery process and the generality of the interior point method.
The main contributions of this paper are: 1) A nonconvex and nonlinear model is established considering complex processes, such as mass linear, mass nonlinear, volume linear, and volume nonlinear intersection flow blending of crude and middle material physical properties with delta-base structure and physical property transfer.
2) An improved interior point method is presented to solve the nonconvex and nonlinear model in the allowable time.
3) Numerical experiment shows that the result of the improved interior point method is superior to that of the commercial solver.

Material flow in crude oil refinery
The main difficulty of production optimization in the crude oil refinery mainly lies in the significant amount of different raw materials and final products and the complex mutual supply relationship between various devices.Also, to meet different requirements for different final products, it is necessary to perform secondary processing of raw materials and intermediate components for selecting different processing options.It is challenging to construct a detailed optimization model for the above system.Therefore, this paper simplifies the entire industrial production process into five parts: raw material procurement, atmospheric and vacuum unit processing, secondary unit processing, blending, and product sales.By refining the transformation relationships within and between these five parts, the entire system can be well described through a relatively simple and neat model.It is noticeable that no necessary procedure is omitted.The performed simplifications include a summary and refinement of the actual process.The material quantity of inflow and outflow of each device is the main decision variable in the optimization model.These decision variables should meet specific equality or inequality relations caused by production requirements.Meanwhile, there are limited materials processed in each device and system due to the device capacity, production technology, and product delivery indicators.In addition, the inflow and outflow of devices should meet specific constraints caused by processing capacity and physical indicators.Finally, raw materials and final products should meet the business operation constraints, such as purchasing, sales, and warehouse storage.Even after some simplification, the great number of raw materials and final products, the complex processing and transformation relations, and the embedded scheme selection will still lead to a nonconvex and nonlinear model, resulting in significant difficulties for the model's formulation and solution.This paper will give the model formulation in detail in Sections 2.2.1-2.2.6 and the solution algorithm of the optimization model in Section 2.3.

X
The amount of the utility system g of link s flowing into device e is the utility system consumption of device e of process s (104 tons), { , , }

X
The amount of the utility system g of link s outflowing into device e is the utility system consumption of device e of process s (104 tons), { , , } In the model, the material set is M , the utility system set is G , and the device set is E .The superscript { , , , , } s gm cy ec th xs  can identify the corresponding production parts, representing raw material procurement, atmospheric-vacuum distillation unit processing, secondary device processing, blending, and product sales, respectively.For example, s M represents the materials that should be processed in production part s and s E represents the devices employed in production part s.In this study, the consumption pairs of the device, material, and utility system comprise set I and the device's output pairs.Material and utility systems consist of sets XX represent the inflow and outflow of one material (utility system) in device e in the production part s , respectively.In the model, subscripts  and  represent the inflow and outflow properties of the material or utility system.

Objective function
The business requirement of refinery optimization implies optimizing the whole process from raw material procurement, production, and processing to product delivery to maximize the economic benefits of enterprises.Therefore, the following objective function is considered for the model to maximize profits: where m p is the selling price of material m , g p is the selling price of utility system g , m c is the purchase price of raw material m, and g c is the purchase price of utility system g .

Material balance limit
Any material must meet the material quantity balance constraint across the whole plant; that is, the sum of its purchased quantity and its quantity produced by all devices should equal the sum of the sales quantity of the material and the quantity of the material consumed by all devices: The atmospheric-vacuum distillation unit is the leading device in refinery enterprises.With crude oil as the main raw material, the output of each side of the atmospheric vacuum device is determined by the amount of processed crude oil and the yield of each fraction of crude oil after distillation.
The secondary unit is a general term for other processing units in refinery enterprises except for the atmospheric-vacuum distillation unit, including residual oil hydrogenation, catalytic cracking, catalytic reforming, hydrocracking, delayed coking, and other types.A set of secondary devices can have many production schemes, where the production process and the nature of raw materials determine each production scheme of raw material consumption and product yield.The delta-base structure of the secondary device should be established to describe the above process accurately [25,26].
where P is the set of alternative processing options, q is the quality attribute, d is the delta- base, ,, is the change rate of the product yield with the feed's quality attribute.The above constraints indicate that the product yield is based on the base and can change with the feed's quality attribute.
The feed of the atmospheric-vacuum distillation unit and the secondary device is usually mixed with multiple strands of materials with different properties.Limited by the production process, the properties of mixed raw materials should meet upper and lower ranges.Additionally, gasoline, diesel, and other final products produced by refinery enterprises are usually mixed with various components.The mixed product properties should meet specific upper and lower limits to meet national or industrial standards.In this paper, the above mixing processes are collectively referred to as blending and are characterized by a virtual device called blending pools.For each blending pool, the amount of material in and out should meet the following constraints:

Processing capacity of the device
The processing capacity of the device defines the maximum and minimum processing capacities of the production device.The device is not allowed to operate beyond the upper or lower capacities in the actual production process to ensure safe production.
In this paper, the processing capacity set of device is denoted by J .Generally, the consumption of primary raw materials characterizes the processing capacity of a device, then for each device, the processing capacity required for processing all major raw materials should not exceed the upper and lower limits; that is, it must meet the following constraints:

Physical properties limit
Since the quality attribute calculation between the lateral line and the blending discharge is similar, the blending equation is adopted to characterize the quality attribute calculation.Specifically, for the atmospheric-vacuum distillation unit, the quality attribute of one lateral line is equal to the weighted average of the quality attribute of the lateral line produced by each processed crude oil.For the blending pooling, the quality attributes of the discharge are equal to the weighted average of the quality attributes of each feed.
In this paper, Q is the set of quality attributes and ( ) , The nature of the secondary plant discharge, the primary feed source for the blending pool, usually varies with the processed raw material nature.It is necessary to establish the quality attribute transfer structure of the secondary device to describe the above process accurately [25,26]: ) where , , , properties in secondary devices are not subject to upper and lower limits.

Restrictions on the utility system
The utility system is a general term for auxiliary facilities to maintain the normal operation of refinery production and processing equipment, also known as supporting works, usually covering the water supply and drainage system, power supply system, compressed air system, steam power system, and wind power system.The production and consumption of the utility system occupy the main variable processing cost of the refinery, which should be modeled and optimized.To reflect the production and marketing process of the refinery utility system, the structures of procurement, utility sales and production, and consumption of the utility system are established in the optimization model of the refinery planning.Common utility systems include fresh water, recycled water, desalted water, deoxygenated water, power, steam, industrial air, instrument air, nitrogen, fuel gas, and fuel oil.In this scenario, the utility system is produced and consumed only in the atmospheric-vacuum distillation unit and secondary plant.The output or consumption of the utility system and the material consumption of the process are expressed with linear relations: ( ) where  represents the correspondence between the produced or consumed utility system and the number of processed materials.
The utility system of the whole plant should balance the output and consumption: 2.2.6.Operational limits In practice, raw material procurement and finished product sales should be within a given upper and lower range: The set of the utility system to be procured is gm G and the set of the utility system available for sale is xs G .The procurement and sales volume of the utility system should meet the upper and lower limits: Accordingly, the optimization model for this problem is established.The solution algorithm is presented in the next section.

Algorithm design
To facilitate the discussion of algorithm design, the optimization model is summarized in the following form: where a and b are parameters in Eq (1); the decision variable X corresponds to f is a linear equality constraint on X corresponding to Eqs (2)-( 4), ( 6) and ( 14)-( 16); 2 f is a linear inequality constraint on X corresponding to Eqs (7), ( 17) and ( 18); 1 h is a linear equality constraint on Y corresponding to Eq (13); 2 h is a linear inequality constraint on Y corresponding to Eq (12); and 1 g is a nonlinear equality constraint on X and Y corresponding to Eqs ( 5) and ( 8)- (11).
The optimization model is nonconvex and nonlinear because the chemical reactions can be generally described with a nonlinear model in actual industrial production, such as the delta-base structure in secondary devices and the quality attribute calculation of blending pools.The above problem coupled with the huge number of raw materials and final products makes it challenging to solve the model.Although many commercial or open-source solvers provide nonconvex optimization model solving, their solving efficiency is insufficient to deal with large industrial optimization problems.Therefore, to reduce the solving time of the model, a high-quality initial solution-finding method is proposed based on the problem structure and actual business logic, then this method is combined with the classical interior point method to solve the problem.The interior point method has a guarantee of convergence efficiency [27,28], especially if the initial solution is good.
The similarity between the initial and optimal solutions can significantly influence the efficiency of the interior point method.This paper generates the initial solution from two linear sub-models, whose solutions are combined to form the initial solution.The proposed algorithm consists of three parts: material (linear model) for solving the initial value corresponding to model ( 2P ), the solution of initial values of quality attributes (linear model) corresponding to questions ( 3 P ), and the entire model corresponding to model ( P ).Among them, the model ( 2P ) ignores the quality attribute based on the model P .Based on constraints such as the raw material purchasing prices, selling price of the product, upper and lower limits of purchase quantity, and upper and lower limits of sales, material balance, the raw material purchasing amount to maximize revenue, product sales, and the flow distribution between devices are determined.Model ( 3 P ) relaxes the nonlinear constraints based on model ( P ), and to avoid obtaining solutions that are grossly infeasible, we refer to the idea of penalty functions [29].Additionally, based on the material quantity calculated in the ( 2 P ) model, the quality attribute value that meets the quality attribute constraints is calculated.Finally, the model P takes the solutions of ( 2 P ) and ( 3 P ) as the initial points and employs the interior point method.The ( 2 P ) and g , respectively.The relaxation cost is taken as the objective function.Minimizing the relaxation cost ensures the satisfaction of the physical constraints as much as possible.The physical constraints can be satisfied by minimizing the slack cost under the determined raw material, product yield, and flow rate between the devices.
For the complete algorithm framework based on the interior point method, we present the solution approach in the form of the following pseudocode: Improved interior point method algorithm Initialization: model ( P ) is established according to actual business scenarios and parameters.

1)
Based on model P , the models ( 2 P ) and ( 3 P ) are generated 2) The optimal solution X is obtained by solving model ( 2 P ) 3) X is kept fixed, and the optimal solution Ŷ is obtained by solving model Taking 0 X and 0 Y as the initial solutions of the interior point method, the model ( P ) is solved.Output: the optimal solution * X and * Y

Case study
This section employs some actual data to establish the model and applies the above algorithm to solve the model.The simulation results indicate the superiority of the proposed algorithm to commercial solvers.

Experimental instances
The examples come from two refiners owned by an oil company.These refiners should consider more than 30 kinds of physical properties when optimizing the production plans, including density, sulfur content, acid content, residual carbon content, metal content, octane number, olefin content, aromatics content, freezing point, and viscosity.This paper selects the refining and petrochemical processes of two refineries A and B for the optimization process.Refinery A contains 36 kinds of physical properties and 79 sets of units, including 3 sets of atmospheric-vacuum distillation units, 50 sets of secondary units, and 26 sets of blending pools.Refinery B contains 54 kinds of physical properties and 58 sets of units, including 2 sets of atmospheric-vacuum distillation units, 40 sets of secondary units, and 16 sets of blending pools.After applying the proposed model, Case 1 (refinery A) has 7229 continuous variables and 3926 constraints, 1827 of which are bilinear.Case 2 (refinery A) has 7230 continuous variables and 3308 constraints, 1851 of which are bilinear.Case 3 (refinery B) includes 5365 continuous variables and 2339 constraints, 379 of which are bilinear.The above information is summarized in Table 3.

Result
As shown in Table 1, all three examples in the numerical experiment have more than 5000 continuous variables and more than 2000 constraints, many of which are bilinear.Since a nonlinear model with such a scale is extremely challenging for a general solver, this paper selects the Gurobi 9.5.1 for comparative evaluation.It improves the original spatial branch-and-bound algorithm and optimizes the solving efficiency of the pooling problem in the refinery.Gurobi's built-in spatial branch-and-bound algorithm can guarantee global optimality theoretically.The proposed algorithm is evaluated by combining the initial solution module with the interior point method solver IPOPT 3.14.5.The test environment is an Intel 2GHz i5 processor with 16 GB memory and MacOS operating system.The results are displayed in Table 2.For Cases 1 and 2, Gurobi cannot find feasible solutions within s, while the proposed improved interior point algorithm provides feasible solutions within 1 min.Meanwhile, the results of model 3 indicate that the proposed algorithm only needs 8 s to find a feasible solution, and the objective function value of the algorithm is nearly 10% higher than that obtained by Gurobi.Moreover, for Case 3, the theoretical upper bound of the spatial branch-and-bound is also obtained to evaluate the optimality of the proposed algorithm.The difference between the objective function value and the theoretical optimal value is less than or equal to 2.3%.Furthermore, in the experimental process, this paper also verifies the improvement degree of the initial solution module to the solving efficiency of the interior point method.As shown in Table 3, when the initial solution is not set in Cases 1-3 while directly employing the interior point method to solve the problem, a feasible solution cannot be obtained within 1000 s.It can be concluded that the designed initial solution module significantly improves the computational efficiency of the overall interior point method.

Conclusions
This paper employed the traditional P model to formulate the pooling structure in the complex refinery material flow network.Furthermore, an algorithm framework composed of an initial solution strategy and interior point method was designed to obtain and improve the efficiency of the solution.

XXXX
of the utility system flowing into the device G set of the utility system types H set of unit processing capacity s H set of unit processing capacity in link $s$ M set of material types s M set of material types that should be processed in link s Continued on next page of the device output utility system s m The amount of material m in link s (104 ton), { , } s gm xs  , ss mg XX The upper limit of material m (utility system g ) in link s (104 tons), { , } s gm xs  , ss mgXXThe lower limit of material m (utility system g ) in link s (104 tons), { , }The amount of material m flowing into device e in link s is the raw material consumption of device e in process s (104 tons), { , , }The amount of material m outflowing into device e in link s is the raw material consumption of device e in process s (104 tons), { , , } s cy ec th The amount of material m in link s flowing into device e by the processing scheme

qqq
value of the quality attribute flowing into material m in link s , { , , } value of the quality attribute outflowing material m in link s , { , , } value of the quality attribute of material m flowing into device e limits of the q value of the quality attribute of material m in link s .raw material procurement, atmospheric-vacuum distillation unit processing, secondary device processing, blending, and product sales.subscript g Consumption of the utility system m amount of material used


is the yield of the fraction out m corresponding to the crude oil cy in mM  .For each atmospheric-vacuum distillation unit, the amount of material in and out should meet the following requirements: , , ,

out e p m 
is the benchmark yield of material m processed by scheme p of refinery unit e , is the reference value of the quality attribute of feed, and , , , of all devices in link s .The upper and lower limits of processing capacity h of device e are denoted as s h J , s h

.
value of the quality attribute of material m in link s .The set of quality attributes satisfying the blending relation of the weight is denoted as weight Q and the set of quality attributes the volume blending relation is volume Q The quality attribute calculation in the blending process should satisfy the following equality constraints: coefficients of the transfer equation of the quality attribute q between the feed in m and the discharge out m of secondary device e .Physical

( 3 P 2 P 3 P ), 1 SL + and 1 SL
) models are described as follows: ) is the initial model P without quality attribute variables and constraints; 3 P is a model in which the initial model P fixes a certain set of material values and only considers quality attribute constraints.In the model ( − are the positive and negative relaxation variables of function 1 ) ) ) According to the requirements of the factory standards or the device's tolerance to the raw materials, the physical properties of the blended materials should meet specific upper and lower limits.aretheupper and lower limits of the q value of the quality attribute of material m in link s .The material's physical value should meet the following constraints: 

Table 1 .
The model scale.

Table 2 .
Comparison between the results obtained by the solver and the proposed algorithm.

Table 3 .
Comparison of the results of the interior point method algorithm with and without the initial solution.