Blending Control for Fuel Production Yann Creff IFP Lyon - - PDF document

SLIDE 1

Blending Control for Fuel Production

Yann Creff IFP Lyon CEA-EDF-INRIA School Optimal Control: Algorithms and Applications May 30 - June 1st 2007 - INRIA Rocquencourt, France

Fuels are not directly extracted as parts of crude oils, but are produced by blending several

• components. These components come from transformations of petroleum cuts in refining units,

the cuts resulting from a first separation of the crude in atmospheric and vacuum distillation

• columns. Combining the available components in the right proportions gives the final blend the

properties required to cope with the standards, for gasolines or diesel fuels. The blend property control problem has received much attention in the last decades: it is a key point in refining

• perations. Under simplifying assumptions, that often hold, this problem can be treated as a

linear control problem. However, part of these simplifying assumptions do not hold for a large class of blending processes, when non linearities appear that are due to transport phenomena. We describe these nonlinearities and propose a way to deal with them for control purposes. The work presented in this paper is the result of a collaboration on blending control between Meriam Ch` ebre (Total), Nicolas Petit (ENSMP/CAS), Julien Barraud and the author (IFP).

1 Process description and control problem

In a blending process, n components are pumped, generally from storage tanks, to a blender (an example is given in figure 1). Downstream the blender, the product either fills a final product tank, or is directly sent off the refinery through a pipeline (a truck, a train...). The combination

• f the components properties makes the blend properties. The proportion of each component

is expressed on a volume basis, through the ratio of its volume flow rate to the blender total volume flow rate. These proportions, grouped in an n−dimensional vector u, form what is called the blend recipe. The blend properties form an m−dimensional vector y, as well as the components properties, denoted bi for component i. 1

SLIDE 2

From a stationary standpoint, the recipe ¯ u, the blend and component properties ¯ y and ¯ bi are linked by an algebraic relation. For most of the physical properties (flash point, cloud point,

• ctane number...), this relation is nonlinear. However, it is (almost) always assumed that for

each property, a transformation can be found that linearizes the blending rule. From now on, using such blending indices in place of the physical properties, we consider linear blending rules. Denoting ¯ B the (m × n) matrix obtained by the concatenation of the n column vectors ¯ bi, we then have ¯ y = ¯ B¯ u. Knowing the component properties, the available quantity of each component, and the required blend properties, the scheduling department of the refinery defines a recipe. The simplest control law (still used when no on-line continuous measurements of the blend properties are available), only consists in applying recipes, that is to say ensuring a proper regulation of each component ratio. This ratio control is a simple combination of flow controllers. The computed recipe includes margins: the corresponding product properties are set to better values than those required by the standards. So, the blend remains acceptable in case of inaccuracy on the blending rules and on the components properties. Some measurements are achieved during the blend, to modify the recipe whenever necessary, and compensate for drifts or sudden changes in the components properties. As the standards were hardened, to cope with new environmental regulations, it became clear that this kind of monitoring was not acceptable anymore. As the number of properties to control increased, providing sufficient margins became difficult: many blends had to be reprocessed. Of course, producing blends with properties beyond the prescribed standards (”give away”) also reduced the benefits. The progressive introduction of on-line analyzers totally modified the picture for blending control. Feedback is now possible, offering a much better operation of the blending process. In the next section, we present the blending control problem classically addressed today. We also briefly show how it is solved. This solution relies upon the fact that each component is directly linked to the blender through a dedicated pipe, which greatly simplifies the dynamics. However, all the blending flowsheets are not that simple. It often occurs that portions of pipes are shared between components. That induces dead volumes effects (called pre-blends). The blending dynamics then present specific delays. They are described in the third section. The fourth section proposes a solution for the control law to take them into account.

2 Solution without pre-blend

As all the components simultaneously reach the blender, for a current recipe u(t) and com- ponents properties ¯ B, the blend properties y(t) are given by y(t) = ¯ Bu(t). The available knowledge on the components properties being simply denoted B, Bu(t) is an estimate of the blend properties. The control objective is either to keep these properties at prescribed references yr, or to enforce bounds, ym ≤ y ≤ yM. 2

SLIDE 3

The recipe u(t) is the control vector. As it is made of proportions, we must have n

j=i ui(t) = 1.

Moreover, each element ui(t) is subject to minimum and maximum bounds um

j and uM j , that

summarize

• the minimum quantity that must be used for the current blend, to avoid the overflowing
• f the storage tank, filled by upstream refining units;
• the maximum quantity to be used for the current blend, to prevent component shortage

for the current and the future blends;

• the minimum flow rate through the pipe connecting the storage tank to the blender. This

non-zero flow rate corresponds to a started pump;

• the maximum flow rate for the pipe, given by the characteristics of the pump and by the

geometry of the pipe.

• the maximum increments and decrements between two executions of the algorithm.

The blend property vector y(t) is the output. We can group the blend properties associated to references in yr (r−dimensional, r ≤ m). They correspond to an (r × n) sub-matrix Br. Similarly, we can group the blend properties associated to minimum bounds in ym (resp. max- imum bound in yM). They correspond to a sub-matrix Bm (resp. BM). Note that it is a priori possible to associate a reference and bounds to the same blend property. Given an objective function J(u), the control problem writes min

u J(u)

           0 < um ≤ u ≤ uM < 1 n

i=1 ui = 1

Bru = yr Bmu ≥ ym BMu ≤ yM (1) The objective function is chosen to promote or avoid the introduction of some components. This can be achieved through linear or quadratic objective functions. For a restricted number

• f properties and components, the LP formulation for gasolines can be traced back to the early

60’s. Let us assume for the moment that this problem is feasible. If its solution is applied, except the case of perfect blending rules and perfectly known components properties, the measured properties differ from their estimations Bru, Bmu and BMu, used to compute the recipe u(t). So, feedback from the measurements y(t) is required. The most common way to proceed is to add biases to the estimations used in problem (1). These are computed by filtering, at first

• rder, the differences between the measurements and their unbiased estimations.

3

SLIDE 4

A priori, the value of m with respect to n is arbitrary. This of course can restrict the possible choices for yr, as feasibility implies r ≤ n − 1. However, even if the initial recipe provided by the scheduling department is compatible with the constraints, the feedback loop, modifying the estimations, can induce non-feasibilities. The way non-feasibilities are handled differentiates the algorithms (more or less easy tuning), as at least the solution of a relaxed problem must be provided each time the algorithm is run (its execution period is about 1 minute). The simplest way, also the hardest to tune, is to include all the constraints in the objective function and to set different weighting factors for each part of it. To summarize, the control problem (1) can be seen as an example of linear predictive control, where the predictions are particularly easy to compute, due to the nature of the system’s ”dynamics”. The situation is far different when pre-blends are considered, as detailed in the next sections.

3 Modelling the pre-blends

We call pre-blends portions of pipe shared by two or more components to reach the blender. A typical example is given by the flowsheet in figure 1. Pre-blends are present when storage tanks, located at various places throughout the refinery, are to be used in the same blending

• flowsheet. This is not common in new refineries, but frequently found in existing assets, either

for the production of final products, or for blenders located upstream transformation units. To illustrate how difficulties appear, let us consider the following single pre-blend in figure 2. There are 3 components, each one having properties bi. The first two components share a pipe volume V before they reach the blender. At time t, y(t) = (u1(t) + u2(t))yp(t) + u3(t)b3 = Qp(t) Q(t)

• yp(t) + u3(t)b3

where Q(t) is the total flow rate, Qp(t) is the flow rate from the pre-blend, and yp(t) represents the properties at the pre-blend outlet. The flow through a pipe can be considered as a pure transport phenomenon. If the flow rates change for components 1 or 2, Qp adapts instanta-

• neously. But it takes time, for the properties of the blend of components 1 and 2 to cross the

volume V . According to the notations in figure 2, we have yp(t) = u′

1(t−δ(t))b1+u′ 2(t−δ(t))b2 =

1 u1(t − δ(t)) + u2(t − δ(t))(u1(t−δ(t))b1+u2(t−δ(t))b2) The delay δ(t) is such that V = t

t−δ(t)

Qp(τ)dτ Finally, y(t) = (u1(t) + u2(t))yp(t, δ(t), b1, b2) + u3(t)b3 4

SLIDE 5

Tank 4 Tank 2 Tank 3

• Blender

Tank 5 Tank 1

Volume V1 Volume V3 Volume V2

Figure 1: A blending flowsheet with pre-blends in place of y(t) = u1(t)b1 + u2(t)b2 + u3(t)b3. The pre-blend dynamically limits the variations

• f the blend properties. This can be extended to more complicated flowsheets.

Table 1 represents the blending flowsheet in figure 3. It is a particular case that proposes a way to define a flowsheet with pre-blends. Each column corresponds to a particular product: it is either a component, an intermediate product (the outlet of a pre-blend), or the final product (the blend). Lines do not concern components. In the column corresponding to a product (except the final blend), there is a unique negative number, standing for the consumption of this product. The line where this negative number is located gives the nearest downstream product in the blending flowsheet: it is either an intermediate product or the final product. The fact that there is a unique negative number in each column means that no product is conveyed to the blender through two or more flows. This assumption always hold in practical cases. The entries in a column sum to

• -1 for components (negative to indicate consumption);
• 0 for intermediate products (produced and consumed);
• n for the final product (number of components in the final product).

5

SLIDE 6

✲ ✲

b3, u3 = 0.6 b2, u2 = 0.3, u′

2 = 0.75

b1, u1 = 0.1, u′

1 = 0.25

yp, Qp = 40 y, Q = 100 V

②

Blender Figure 2: Simple blender For the lines (intermediate or final product),

• negative numbers indicate which components and intermediate products compose the

product;

• the unique positive number is located in the column corresponding to the product. It

indicates the number of components entering its composition. So, the entries in a line sum to 0. To follow a component,

• look for the intermediate or final product that uses it;
• go to the column corresponding to this product, see how it is used . . . up to the final

product. For flowsheets without pre-blend, the array reduces to a single line and n+1 columns. The first n columns, corresponding to the components, are filled with 1, while the last one, corresponding to the product, contains n. Note that the sub-matrix obtained by eliminating the columns that correspond to components can always be chosen lower-triangular, providing a proper numbering

• f the pre-blends. This is helpful to develop blending simulators.

Let us consider a flowsheet with p pre-blends (1 to p) and denote:

• Qi(t) the volume flow rate of component i at time t, with i ∈ {1, . . . , n};

6

SLIDE 7

r r r r r r ✉

u1 u2 u3 u4 u5 u6

Pre-blend 1 Pre-blend 2 Pre-blend 3

Blender y Storage tank Figure 3: Blending flowsheet with 6 components

• Qn+i(t) the total volume flow rate through pre-blend i at time t, with i ∈ {1, . . . , p} (the

inlet and outlet flow rates are always equal);

• Q(t) the total volume flow rate of the blender at time t, Q(t) =

i=1,n Qi(t);

• Vi the dead volume corresponding to pre-blend i;
• bi the properties of component i, bE

j (t) (resp. bS j (t)) the inlet properties (resp. outlet

properties) of pre-blend j at time t. For a component, we consider bi = bE

i (t) = bS i (t).

To each component i, we associate a path Πi, defined by the series of pi pre-blends from this component to the blender. This is a series of pi different integers, that refer to the pre-blend numbering indices, Πi = {π1

i , π2 i , . . . πpi i }, with πj i ∈ {1, . . . , p} for all j ∈ {1, . . . , pi}. Πi = ∅,

i.e. pi = 0, means that component i is directly linked to the blender. For the case described in figure 3, we have:

• Π1 = {3}, p1 = 1;

7

SLIDE 8

comp 1 comp 2 comp 3 comp 4 comp 5 comp 6 PB1 PB2 PB3 P PB1 −1 −1 2 PB2 −1 −1 2 PB3 −1 −2 3 P −1 −2 −3 6 Table 1: Array representing the flowsheet in figure 3

• Π2 = {2, 3}, p2 = 2;
• Π3 = {2, 3}, p3 = 2;
• Π4 = ∅, p4 = 0;
• Π5 = {1}, p5 = 1;
• Π6 = {1}, p6 = 1.

For each pre-blend i, we define the set Γi of its qi inlet flow rates. This is a set of qi different integers, that refer to the numbering indices of the volume flow rates, Γi = {γ1

i , γ2 i , . . . γqi i },

with γj

i ∈ {1, . . . , n + p} for all j ∈ {1, . . . , qi}.

For figure 3, we have:

• Γ1 = {5, 6}, q1 = 2;
• Γ2 = {2, 3}, q2 = 2;
• Γ3 = {1, 6 + 3} = {1, 9}, q3 = 2.

For Πi = ∅, properties bi are simply weighted by Qi(t)/Q(t) in the relation giving the blend properties as a linear combination of the component properties. Now consider cases where Πi = ∅. For pre-blend πj

i , the total flowrate is Qπj

i (t) =

k∈Γπj

i

Qk(t). The inlet property bE

π1

i writes

bE

π1

i (t) =

• j∈Γπ1

i

Qj(t)bE

j

• j∈Γπ1

i

Qj(t) In this term, bi enters through Qi(t)

• j∈Γπ1

i

Qj(t)bi = Qi(t) Qπ1

i (t)bi

8

SLIDE 9

For the outlet, we have bS

π1

i (t) = bE

π1

i (t − δπ1 i (t)), the delay δπ1 i (t) being defined by

Vπ1

i =

t

t−δπ1

i (t)

Qπ1

i (τ)dτ.

(2) In bS

π1

i (t), bi enters through

Qi(t − δπ1

i (t))

Qπ1

i (t − δπ1 i (t))bi.

Similarly, for pre-blend π2

i , we have

bE

π2

i (t) =

• j∈Γπ2

i

Qj(t)bE

j

• j∈Γπ2

i

Qj(t) , i.e., for the bi term in bE

π2

i (t),

Qπ1

i (t)

Qπ2

i (t)

Qi(t − δπ1

i (t))

Qπ1

i (t − δπ1 i (t))bi.

For bS

π2

i (t) = bE

π2

i (t − δπ2 i (t)), bi appears as

Qπ1

i (t − δπ2 i (t))

Qπ2

i (t − δπ2 i (t))

Qi(t − δπ2

i (t) − δπ1 i (t − δπ2 i (t)))

Qπ1

i (t − δπ2 i (t) − δπ1 i (t − δπ2 i (t)))bi.

We see that compositions of pure delays appear on paths Πi. Define function ∆j

i(t) : t →

t − δπj

i (t), for all πj

i in Πi. The composition of these functions for a particular i is given by

∆k,j

i (t) ∆k i (∆j i(t)) : t → t − δπj

i (t) − δπk i (t − δπj i (t))

and ∆l,k,j

i

(t) ∆l

i(∆k,j i (t)). According to these notations, the term including bi in bS π2

i (t) writes

Qπ1

i (∆2

i (t))

Qπ2

i (∆2

i (t))

Qi(∆1,2

i (t))

Qπ1

i (∆1,2

i (t))bi.

At the outlet of the last pre-blend πpi

i , we have for bi in bS πpi

i (t)

Qπpi−1

i

(∆pi

i (t))

Qπpi

i (∆pi

i (t))

Qπpi−2

i

(∆pi−1,pi

i

(t)) Qπpi−1

i

(∆pi−1,pi

i

(t)) · · · Qπ1

i (∆2,...,pi

i

(t)) Qπ2

i (∆2,...,pi

i

(t)) Qi(∆1,2,...,pi

i

(t)) Qπ1

i (∆1,2,...,pi

i

(t)) Finally, for the blend, denoting Ui(t) the weight for bi Ui(t) = Qπpi

i (t)

Q(t) Qπpi−1

i

(∆pi

i (t))

Qπpi

i (∆pi

i (t))

Qπpi−2

i

(∆pi−1,pi

i

(t)) Qπpi−1

i

(∆pi−1,pi

i

(t)) · · · Qπ1

i (∆2,...,pi

i

(t)) Qπ2

i (∆2,...,pi

i

(t)) Qi(∆1,2,...,pi

i

(t)) Qπ1

i (∆1,2,...,pi

i

(t)). (3) 9

SLIDE 10

Then, as ui(t) Qi(t)/Q(t), with Ui(t) = ui(t) for Πi = ∅, the outlet writes y(t) =

n

• i=1

Ui(t)bi. (4) Dynamically, the blend properties still write as a linear combination of the components proper- ties, the weights being functions of the past and current component volume flow rates. As the past values of Qi(t) are known, it is possible to compute all the delays, through their implicit definition (equation (2)). So, the functions ∆j

i are defined. We can compute Ui(t) to provide

an estimation of y(t). This is used to design control laws.

4 Proposed solutions for pre-blends

4.1 ”Infinite horizon” control

As mentioned earlier, providing an estimation of the blend properties at current time is possible through (4). This estimation can be used to feedback a difference between the measurement and the estimation, to avoid control biases. To take into account the equality and inequality constraints on the blend properties, the situation is more complicated, because now the blend properties nonlinearly depend on the control. A first approach can be based upon the fact described below. Consider a case where there is a single pre-blend between the storage tank for component i and the blender. According to (4), we have Ui(t) = Qπ1

i (t)

Q(t) Qi(∆1

i (t))

Qπ1

i (∆1

i (t))

Denote t0 the current time. If the recipe u(t) and the total flow rate Q(t) are kept constant from t0, Qπ1

i (t) = Qπ1 i (∆1

i (t)) and Qi(t) = Qi(∆1 i (t)) for t ≥ Vπ1/Qπ1

i (t0). Then Ui(t) = Qi(t)/Q(t) =

ui(t). If we look beyond the time needed to flow through the pre-blend, the expression for Ui(t) is the expression corresponding to the case without pre-blend. This extends to any flowsheet. If u(t) and Q(t) are kept constant from t0, we can define a delay ˜ δ that corresponds to the longest (slowest) path linking a set of components (at least two for a pre-blend to exist) to the blender. For t ≥ t0 + ˜ δ, y(t) = n

i=1 ui(t)bi, as in the case without

pre-blend. Providing a convenient feedback is used, based upon (4) to properly synchronize the estimation and the measurements, we can still see (1) as a linear predictive control for a flowsheet including pre-blends. A single move is computed and the output constraints are enforced on a single point in the future, that corresponds to the stationary point. Indeed, in practice, linear predictive 10

SLIDE 11

controllers are often tuned in a similar way. This approach provides good performances in many practical cases. However, even if the problem keeps feasible, constraints can be transiently

• violated. This is critical when the product is directly sent off the refinery, so we try to find

implementable solutions to improve the situation.

4.2 ”Finite horizon” control

Our work on this topic is recent. What is described below corresponds to first results. As in the previous case, we solve a discretized version of the original problem: in practice, the blending control algorithm is run about every minute. Typical delays due to pre-blends range from 15 to 45 minutes. The duration of a batch blend (when the blend fills a final product tank) obviously depends on the tank volume, but is most often greater than 4 hours. Collecting the Ui(t) in an n−dimensional vector U(t), we can write U(t) = Γ(t)u(t) and y(t) = BΓ(t)u(t), where B and u(t) keep their definitions, introduced in section 1. The non-constant entries of Γ(t) are only functions of past values of u. For a proposed u(t0) at current time t0, we can compute all the delays (they depend on u(t0)) and fill Γ(t0) by picking past values of u at the right locations in historical tables. When predictions are done from a given trajectory u(t) for t ≥ t0, attention must be paid to the fact that the structure of Γ(t) varies as t increases. The involved mechanism has been mentioned in the previous section. For instance, in the simple case of figure 2, we have Γ(t0) =           u1(t0 − δ(t0)) u1(t0 − δ(t0)) + u2(t0 − δ(t0)) u1(t0 − δ(t0)) u1(t0 − δ(t0)) + u2(t0 − δ(t0)) u2(t0 − δ(t0)) u1(t0 − δ(t0)) + u2(t0 − δ(t0)) u2(t0 − δ(t0)) u1(t0 − δ(t0)) + u2(t0 − δ(t0)) 1           If u(t) = u(t0) for t ≥ t0, then, for t ≥ δ(t0), Γ(t) = I3, Ik denoting the k−dimensional identity

• matrix. The number of different structures to consider for Γ, from Γ(t0) to In, increases with

the number of pre-blends. These structures appear as time increases and components (or intermediate products) are released from dead-volumes. The times corresponding to switches for the structure of Γ can be computed for any given trajectory on u(t), as knowing such a trajectory is sufficient to compute all the delays in the system. So, it is possible to modify problem (1) into a problem where:

• the optimization variable is a set of Hc vectors u(tk), with k ∈ {0, 1, . . . , Hc − 1}. For an

execution period given by T, to each k corresponds an integer nk, increasing with k, with n0 = 0 and tk = t0 + nkT. 11

SLIDE 12

• the bounds on the variables are checked for each k. Bounds on control variations can be

included through linear constraints involving couples of successive ui(tk).

• the constraints are enforced at Hp times tl, with l ∈ {0, 1, . . . , Hp − 1}. To each l cor-

responds an integer nl such that tl = t0 + nlT, and a matrix Γ(tl), computed at the beginning of the function dedicated to the evaluation of constraints. Several Γ(tl) can have the same structure. They could also become equal for the same value of u(tk) and sufficiently large nl, but these cases are easily trapped and removed from the optimization problem. Of course, as Hc and Hp increase, the problem quickly becomes intractable in practice. So, we have decided to sequentially test different options. So far, only situations where Hc = 1 have been studied: we want each control applied on the process to correspond to a stationary point that respects the constraints. We first tried the simplest case, involving Γ(t0) and In. Compared to the solution of problem (1), it only adds a guarantee on the blend properties, just after the new control is applied. As forecast, we found that this was not sufficient. We then tried various combinations, more or less efficient, some of them including all the possible structures for Γ. Finally, we decided to uniformly distribute a given number of points on the horizon from current time t0 to t0 + ˜ δ. This is a common practice in predictive control. We have obtained good results, for example with 3 points for 1 pre-blend and 5 points for 3 pre-blends. An example is shown in the next section. However, we still lack a rule to robustly define the number of points and their locations. Our current work investigates this problem, with Hc = 1. We suspect that taking Hc = 1 could be the source of difficult problems: this point will be addressed later.

4.3 Some results

Simulation results are presented in figure (4). They show a comparison between the two control approaches described in the previous section, applied on a problem concerning the production

• f a diesel fuel. The simulation, that accurately represents the transport phenomena, is based
• n the description in table (2), and corresponds to figure 1.

comp 1 comp 2 comp 3 comp 4 comp 5 PB1 PB2 PB3 P PB1 −1 −1 2 PB2 −1 −1 2 PB3 −1 −2 3 P −1 −2 −3 5 Table 2: Array used for the simulations. It represents the flowsheet in figure 1 12

SLIDE 13

In both cases, the feedback uses the same synchronization between the estimations and the measurements (it is mandatory to get good performances). Initially, the distillation point and the cetane number are below their minima. The control laws, based on the same SQP algorithm, both find a stationary recipe that respects the given bounds. However, the finite horizon approach provides better transient results. It uses 5 ”coincidence points” uniformly distributed over the horizon [t0, t0 + ˜ δ] from current time t0. The delay ˜ δ corresponds to the longest path from components to the blender, for the recipe proposed by the optimizer. It is also interesting to notice that the stationary recipes are different in the two cases.

Properties

x−axis = time (min)

MAX Constraint MIN Constraint

Properties

x−axis = time (min)

WITH exact computation of the delays

Properties

x−axis = time (min)

WITHOUT pre−blend taken into account 10 20 30 Density 10 20 30 Cetane number 10 20 30 Sulfur 10 20 30 Distillation 360°C 10 20 30 Flash point

Figure 4: Comparison between ”infinite horizon” control and ”finite horizon” control 13

SLIDE 14

5 Conclusion

We have defined a mean to represent any flowsheet including pre-blends. It serves as a basis for the simulation of blending processes. We have derived an expression for the blend properties, that can be used to synchronize measurements and estimations, which is mandatory for the design of feedback loops. We have shown that the control problem in the case of pre-blends can be seen as an ”infinite horizon” linear predictive control problem. However, the solution of such a problem does not always give satisfactory results in practical cases. To go further, we have proposed a nonlinear solution, that can deal with constraints on the whole trajectory of the blend properties. Its current formulation provides better results, but requires computational times that are not really compatible yet with an online implementation. Work is undertaken to improve this situation. 14