A mathematical model of commercial production of petroleum products is considered, with an optimal control being the one minimizing deviations from the target formulation. The resulting optimization problem of mixed integer programming depends on two types of variables: 1) the flow rate of petroleum product; 2) a discrete value that takes the values 0 or 1, which characterizes the start of the certification in a given time interval in a given commercial tank before shipment to a given receiving object. To solve the problem of mixed integer programming, the following approach is proposed: at each iteration, the value of the certification vector is found using a genetic algorithm and registered. Then, by solving the optimization problem, the optimal control is determined, i.e. the flow rate of the petroleum product. A computational experiment is carried out to calculate the optimal schedule for the proposed process flow chart.
Commercial production of petroleum products (gasolines) is carried out by mixing various components in the
proportions specified by the formulation. The operation schedule of the oil plant equipment should streamline
the components mixing process and be made for a given calendar plan. At the same time, it should take into
account the current state of the equipment and minimally deviate from the estimated production schedule for a
given period. For the effective operation of oil refining enterprises, it is necessary to coordinate the dispatching
schedule with the actual state of production [Hus2021, Zyk2018, Zyk018, S
Petroleum products manufacturing is complicated by the following factors: – if more than one brand of gasoline is needed and all components are available, it is necessary to mix them so that there are no residues;
– it is necessary to take into account the change in the operating modes of some installations; – it is necessary to take into account the distribution of incoming and outgoing flows.
Besides: – the composition of crude oil may change; – needs and prices may change; – depending on uncertain factors (the time passing after the reactor shutdown, the activity of the catalyst, the air temperature, the temperature of the cooling water, etc.), the yields of petroleum products may change.
The schedule should be determined dynamically following the updated calendar task and changes in the actual state of the equipment. Changes can occur at any time, and the dispatching planning system must either find an effective solution or inform about the impossibility of performing an updated calendar task.
The dispatching schedule should be calculated since the components have to be included in the production more fully according to the total amount of the components arriving for commercial production for the entire period of the calendar schedule, provided that all specified restrictions are met. Such an optimal dispatching schedule will ensure a more coordinated operation of the refinery both for the period of the dispatching schedule and for the period of the calendar plan.
At the moment, owing to technical development, it is possible to create and implement advanced control systems for petroleum products manufacturing. The proposed algorithm for calculating the optimal dispatching schedule can be the basis of such a system that combines the work of specialists at different levels of business. 2
The first step in solving the problem of optimal scheduling is to formalize the blending processes in the form of a mathematical model. Consider the following description of the parameters of the mathematical model [Zyk2020, Sav2019, Zyk019].
T ypes = fpipe; pump; tank; pipec; tankg; stockg is a set of values that characterize the type of an object, the numbers of its elements corresponding to the numbers of the object: pipe is the type of object that represents the process piping in the diagram; pump is the type of object that represents a pumping station; tank is the type of an object that represents a tank; pipec is the type of object that represents a components’ blending unit; tankg is the type of object that represents a commercial tank; stock is the type of object that represents an outflow.
Range is the period for which the dispatching schedule is being made; Nproducts is the amount of the end products; Ninputs is the number of objects that are sources of components, p = 1; :::; Ninputs Nstock the number of receiving objects for specific end products, Nstock = Nproducts, s = 1; :::; Nstock; Ntankg is the number of commercial tanks, g = 1; :::; Ntankg ; Npipes is the number of pipes, b = 1; :::; Npipes; Ntanks is the number of component tanks, d = 1; :::; Ntanks; Ntimes is the number of time intervals, t = 1; :::; Ntimes; Nobjects is the number of objects involved in the model:
Nobjects = Ninputs + Npipes + Ntankg + Nstock; i = 1; :::; Nobjects; .
Nl is the number of connections between the model objects (equipment), l = 1; :::; Nobjects.
0f11 f N1l f 2 1 f 2 2 : : : fl2 : : : f N2l : : : f t
1 : : : f2t : : : : : : : : : flt : : : : : : : : : f Ntl : : : f1Ntimes 1 :: :: :: :f:2N:times CCC :: :: :: :f:lN:times CCCA : : : fNNltimes 0P ass1 BP ass2 BB: : : BBP asss B@: : :
P assNstock 1 C C C C C C A flt is the rate of the flow passing through the l-th connection for the t-th interval, l = 1; :::; Nl, t = 1; :::; Ntimes – The matrix P ass: P asss is a matrix in which the columns are the time intervalt, and the rows are the connections linking the commercial tanks g = 1; :::; Ntankg with the outflow s: 0P asss11
P asss1Ntanksg
P asss12 P asss22 : : : P asssl2 : : : : : : P assst1 : : : P assst2 : : : : : : : : : P assslt : : : : : : P asss2Ntanksg : : : P assstNtanksg : : : P asss1Ntimes 1 : : : P asss2Ntimes C : : : : : : CC : : : P assslNtimes CC : : : : : : CC : : : P asssNNttiamnkessg A
For the matrix elements, the condition is met: passsgt 2 f0; 1g. Next, for the genetic algorithm, the matrix is represented as a vector. To do this, first, each matrix is pulled into a vector. Then, the obtained vectors are concatenated. 3
In production, the scheduling process can take several days. Since there may be changes in the calendar task, as well as changes in the state of the equipment, the optimal schedule should be planned as quickly as possible. At the same time, it is necessary to calculate not only the optimal schedule, made with fixed terms of products’ certification but also the possible terms of accelerated preparation of petroleum products. A solution is required that can reduce the resources used by production through moving the period, lasting from the start of the production to the product certification, to an earlier date. The time of the petroleum product manufacturing should be also reduced, if possible. The optimal schedule should reflect the most productive process of manufacturing petroleum products.
Let us consider an approach to problem-solving, based on the use of a genetic algorithm.
Step 1: Initialization of the initial population. Each individual is a binary vector corresponding to a multidimensional matrix P ass that characterizes the beginning and end of the certification for commercial tanks.
Step 2: The initial values of matrix F k are generated, where k = 0.
Step 3: k = k + 1 is taken. The optimization problem is solved for an individual of the population P assk. In this case, solution F k is required that satisfies the constraints described earlier in the mathematical model [Zyk2020, Sav2019, Zyk019]. To assess the optimality of the solution, the values of the quality and mass matrices should be calculated.
Step 4: If there is a solution to the optimization problem, then the fitness of the individual P assk is evaluated. Otherwise, the transition to step 3 occurs.
Step 5: Selection, crossover, and mutation operations are carried out.
Step 6: The condition for stopping the genetic algorithm is checked. If the condition is met, then one individual or a set of individuals with the best fitness scores for all generations is taken as the answer. Otherwise, the transition to step 3 occurs.
To solve the optimization problem, the SciP y library with the optimize package, namely, the minimize method is used. The value of the parameter method = trust constr sets the algorithm selected for solving the optimization problem. The trust constr method is a trust region algorithm for solving optimization problems with constraints. It switches between the two implementations depending on the definition of the problem. This is the most universal minimization algorithm with constraints implemented in SciP y, and the most suitable for large-scale problems. For problems with equality constraints, it is the implementation of the Byrd OmojokunT rust Region SQP method. When inequality constraints are also set, the trust-constr method switches to the interior point method. The interior point algorithm, in its turn, solves the problem with inequality constraints by introducing stock variables and solving a sequence of problems by the method of barrier functions with equality constraints for gradually decreasing values of the barrier parameter. The sequential programming method (SQP ) with equality constraints is used to solve subtasks with an increase in the level of accuracy as the iteration approaches the solution [Byr99].
However, the direct application of the sequential quadratic programming method to the barrier method leads to inefficient primary steps that, as a rule, violate the positivity of weak variables, and, thus, are often interrupted by a restriction of the trust region. The formulation of the sequential quadratic programming iteration and the definition of the scaled trust region allow us to avoid too close an approximation to the boundary of the admissible domain. 5
The barrier method takes an admissible point as an initial one. In the optimal scheduling problem, selecting an initial point for the optimization algorithm in an admissible domain is a timeconsuming task. The following scheme for searching for an initial admissible point is proposed.
The algorithm for finding the initial admissible point (y0; x0) for the problem
y ! min vi(x)
y; i = 1; :::; m; hj (x) = 0; j = m + 1; :::; n: (1) is as follows.
Step 1: Initial point x0 is obtained that satisfies the equality constraints hj (x) = 0; j = (m + 1); :::; n. Step 2: y0 = maxfvi(x0) j i = 1; :::; mg is taken.
Step 3: The condition for stopping the algorithm is checked. If y0 < 0, then vi(x) 0; i = 1; :::; m, and the initial admissible point (y0; x0) for problem (1) is found.
Otherwise, starting from the initial point (y0; x0), the following minimization problem with equality constraints is solved.
y ! min vi(x)
y = 0; i = 1; :::; m; hj (x) = 0; j = m + 1; :::; n: until we get y < 0.
When the solution y < 0 is found, then vi(x ) < 0; i = 1; m, therefore, an admissible point for problem (1) is found. 6
To test the proposed solution of the optimization problem, a computational experiment was conducted. To do this, consider the process flow chart of the petroleum products manufacturing through blending components. This flow chart should reflect the basic logic of the sequence of actions taken during the commercial production of petroleum products (Figure 1).
The chart shows the main objects needed for modelling the commercial production of petroleum products. Each object has the values of its inherent characteristics. The description of the objects is given in Table 1.
The following is the description of the implemented solution. The first step is to solve the problem of finding an admissible point. The initial point is constructed to be passed to scipy.optimize.root function. The initial
Components’ blending unit Commercial tank Receiving object point is matrix F of flow rates along the connections of objects at each moment with the dimension of (56; 12) filled with zeros. The result of the scipy:optimize:root function is a point that satisfies a system of nonlinear equations consisting of a set of equality constraints.
The next step is the scipy:optimize:minimize function, which is used in the following form: fromscipy.optimize import minimize res = minimize(fun y, y f vec.ravel(), method=trust-constr, hess=lambda x: np.zeros((len(y f vec), len(y f vec))), constraints= [linear constraint, nonlinear constraint], bounds=bounds)
Vector f vec is obtained by pulling matrix F into a vector. Vector y f vec is obtained by adding variable y to vector f vec. The constraints of the problem of finding an admissible point are similar to those of the original optimization problem. Consider the imposed boundary constraints on the first value of vector y f vec, i.e. variable y.
To end the search for an admissible point, variable y must become a negative value during the optimization process. Even though the values y should lie within the boundaries ( 1; 0), scipy:optimize:minimize works better with specified boundaries ( 1; 0). In this case, the search for a solution stops when the value ofy becomes close to 1.
The difficulty of restrictions application lies in the way they are set. The scipy:optimize:minimize function with the method = trust constr parameter accepts LinearConstraint objects as linear constraints. In this object, the constraints are set by multiplication AF , where F is the velocity matrix of the form (n), and matrix Ahas the form (m; n). To create a LinearConstraint object, it is necessary to express F from the formulae contraints and then create matrix A.
As a result, the scipy:optimize:minimize function returns an admissible point for solving the main optimization problem. The found admissible point is taken to a function of the following type:
res op = minimize(fuction F, f vec.ravel(), method=trust-constr, constraints= [linear constraint, nonlinear constraint], bounds=bounds)
As a result of the function, we get a solution to the optimization problem. The solution satisfies the constraints specified in the mathematical model. The value of the target function at the optimal point is 1.525 E-12.
The operating time (Table 2) depends on the parameters of the scipy:optimize: minimize function that are taken into account. The function uses the values of the Hessian matrix to find an admissible point. For the considered constraints, the specified matrix is filled with zeros. During the experiment, the Hessian matrix is transferred to the scipy:optimize:minimize function of an additional optimization problem. When transferring the main optimization problem to the function, no improvement in the optimal point construction time was found in this experiment. In this case, the Hessian matrix is passed as a parameter when creating constraint objects for the main and additional optimization problems.
The solution needed for creating an optimal schedule can be obtained with the method used, combining the barrier methods, trust region strategies and sequential quadratic programming, with correctly selected parameters and a pre-calculated admissible initial point. The trust constr algorithm is effective because it avoids too close an approximation to the boundary of the admissible region in search of the optimal solution. The time to develop a solution is conditioned by the set of constraints for the specified process flow chart, as well as by the presence of non-linear constraints, i. e. restrictions on the sequence of filling and emptying of commercial tanks. 7
A computational experiment confirmed the effectiveness of using the scipy:opti mize:minimize function for the problem of optimal scheduling of commercial production of petroleum products. The following further work with this function is required: analysis of quality formulae for the problem of optimal scheduling of petroleum products; implementation of the constraints transformed into the required form due to the selected quality formula; and subsequent implementation of the algorithm for optimal scheduling proposed in the first section using the obtained solution. [Zyk018]
A. V. Zykina, O. N. Kaneva, V. N. Savkin, T. Yu. Fink. Calculation methods for quality indicators of commercial gasolines (petrochemicals). AIP Conference Proceedings, pp. 020030, 2019. Richard H Byrd [et al.]. An interior point algorithm for large-scale nonlinear programming. ISIAM Journal on Optimization, 9(4): 877–900, 1999.