Online gasoline blending with EPA Complex Model for predicting emissions

Abstract

The empirical Complex Model developed by the US Environmental Protection Agency (EPA) is used by refiners to predict the toxic emissions of reformulated gasoline with respect to gasoline properties. The difficulty in implementing this model in the blending process stems from the implicit definition of Complex Model through a series of disjunctions assembled by the EPA in the form of spreadsheets. A major breakthrough in the refinery-based Complex Model implementation occurred in 2008 and 2010 through the use of generalized disjunctive and mixed-integer nonlinear programming (MINLP). Nevertheless, the execution time of these MINLP models remains prohibitively long to control emissions with our online gasoline blender. The first objective of this study is to present a new model that decreases the execution time of our online controller. The second objective is to consider toxic thresholds as hard constraints to be verified and search for blends that verify them. Our approach introduces a new way to write the Complex Model without any binary or integer variables. Sigmoid functions are used herein to approximate step functions until the measurement precision for each blend property is reached. By knowing this level of precision, we are able to propose an extremely good and differentiable approximation of the Complex Model. Next, a differentiable objective function is introduced to penalize emission values higher than the threshold emissions. Our optimization module has been implemented and tested with real data. The execution time never exceeded 1 s, which allows the online regulation of emissions the same way as other traditional properties of blended gasoline.

Keywords

emissions / reformulated gasoline / online control / global optimization

Introduction

The US Environmental Protection Agency ^{(EPA, 2015}) established reformulated gasoline (RFG) guidelines in reducing emissions from gasoline-powered vehicles to satisfy air quality requirements. In 1998, the EPA introduced the Complex Model, which is denoted herein as EPA-CM, to allow refineries calculate the nitrogen oxides (NO_x), volatile organic compounds (VOCs), and toxic air pollutants (TOX) based on 11 gasoline properties. The regulations embedded in the EPA-CM are available in print form ^{(EPA, 2015}). The EPA also publishes and updates an Excel spreadsheet ^{(EPA, 2016}) to calculate the emission reductions compared to a statutory baseline gasoline. This spreadsheet is selected to give the reference equations of the present study. In its original form, the EPA-CM is difficult to implement within mathematical optimization models. Since 1998, considerable effort has been expended (see Appendix for a bibliographical background) to implement the EPA-CM into refinery operations. The first breakthrough occurred in 2008 (^{Furman and Androulakis, 2008}) and 2010 ^{(Misener et al., 2010}), in which EPA-CM was represented as a mixed-integer nonlinear programming (MINLP) using disjunctive programming techniques. As discussed in the next section, the EPA-CM cannot be operated with the original properties of the target gasoline. More specifically, this model utilizes edge target and delta target values derived from the target property values by satisfying a number of “if–then” conditions, which introduce binary variables that must be managed along with continuous variables to mathematically model the emission functions. The MINLP by Furman and Androulakis (2008) is composed of 78 continuous variables and 21 binary variables. The authors’ signal computation times for the optimization solver extend approximately 600 s. The part of the EPA-CM model included in the more general pooling problem treated by ^{Misener et al. (2010)} requires 87 continuous variables and only 15 binary variables. The EPA-CM is part of the pooling problem, and its execution time has not been reported accurately. Nevertheless, an EPA-CM computation time of approximately 140 s is estimated in Appendix A.

The current state-of-the-art approaches cannot control the emission in online blending to the same extent as the typical properties of target gasoline. The present work addresses this problem. As mentioned previously, MINLP model computation times remain prohibitively long for online blending processes. Our online controller ANAMEL (^{Chèbre et al., 2010}) actually calls the optimization solver up to 16 times/min (because of the frequency of changes in component properties) to regulate some 40 properties of the target gasoline; thus, the optimal solution is returned by the solver in less than 4 s.

Our approach is completely different from all preceding contributions. We begin by introducing a differentiable approximation of the EPA-CM emission functions. To this end, sigmoid functions are used as an approximate differentiable variant of if–then logical conditions contained in the original model. The accuracy of this approximation is controlled by using the measured precision of gasoline properties. In this way, we avoid the need for binary variables and mixed constraints. Then, instead of imposing (non-convex) constraints on emission values, we minimize an objective function that penalizes emission values exceeding a given threshold. We thus obtain an optimization problem with a differentiable and non-convex objective function with respect to the classical linear constraints on recipe and gasoline properties.

This new control/optimization problem is implemented and tested with real data. The execution time of the global optimizer solver, which is written in MATLAB, is less than 1 s; this good result paves the way to controlling emissions online the same way as the other traditional properties of blended gasoline. This fast model is also suitable for refinery scheduling departments. The international patent (^{Janaqi et al., 2015}) for our process was awarded in July 2015.

This paper is organized as follows. Section 2 provides a bibliographical review. Section 3 describes the NO_x computation and lists useful definitions in the remainder of the article. Section 4 contains the theoretical development of our approach. Section 5 provides some implementation details and results. Finally, Section 6 presents the conclusion.

Bibliographical review

RFG fabrication is the subject of the first part of a negotiation between the US legislature and refiners in 1991. This negotiation aimed to decrease toxic emissions (i.e., NO_x, VOC, and TOX) of an RFG with respect to a baseline gasoline. This negotiation defined new standards that were subsequently introduced into the calculation models, namely, EPA Simple Model and EPA-CM (in effect since 1998). The following paragraphs chronologically review the papers discussing the economic, mathematical, and implementation aspects of the calculation models.

In 1994, ^{Nocca et al. (1994)} discussed the specifications to be imposed on RFG and the modifications to an upstream trial on a mixture to satisfy the emission constraints imposed by the EPA Simple and Complex Models.

At the same time, ^{Hirshfield and Kolb (1994)} demonstrated the benefits of linking EPA-CM with the mathematical models (i.e., linear and nonlinear) used in refineries. They stressed the need to formulate “simplified reduced forms” of EPA-CM that could easily be integrated with mathematical solvers. They did not succeed however to introduce simplified reduced (linear) forms because these forms, as signaled by ^{Treiber (1998)}, provide highly inaccurate approximations.

In 1995, ^{Trierwiler (1995)} asserted that the emission functions cannot be (globally) linearized with an acceptable level of precision sufficient to be integrated into the refinery linear programming. He developed a tangent surface technique to locally linearize the EPA functions. We observed that this technique, although theoretically solid, can only be applied efficiently in the neighborhood of gasolines x whose emission functions are relatively linear. We found some actual cases with emission functions exhibiting strong variations around some gasolines x. The singular Hessian values at these points vary between 10^-1 and 10³. In such cases, the linear approximation is only valid over an extremely small zone around x; consequently, recipe modification by the solver is nearly zero.

In 1995, ^{Korotney (1995)} conducted a program to test whether some EPA-CM extrapolations might be useful. They showed in several cases that EPA-CM outputs emission values greater than the actual emissions. At first glance, this finding would seem to be extremely good given that if EPA-CM emissions were less than certain thresholds, then the actual emissions would surely comply with the legal requirements. The problem here is that, with such overestimations, the feasible zone narrows and blend cost rises.

During the same period, ^{Jackson and Vittachi (1995)} integrated the EPA Simple Emissions Model into the online control; their model only set bounds for TOX and exhausted benzene. Moreover, they estimated the impact of each component over the two pollutants and introduced this information into the online control. Their work did not integrate the EPA-CM.

^{Cason (1997)} performed a thorough analysis of the influence of emission constraints on the flexibility of gasoline blending. He demonstrated the impossibility of defining easy-to-respect bounds on the blend properties that would imply emissions in compliance with the EPA Standard because a large reduction in blend flexibility implies a needless cost increase. Mathematically speaking, introducing supplementary bounds on gasoline properties dramatically reduces the feasible zone. Cason wrote, “It is doubtful whether a refiner could consistently blend any gasoline with these properties… The industry is effectively stifled by mathematics.”

In the same year, ^{Hirshfield and Kolb (1997)} released a model for estimating the cost of emission reduction as required by the EPA Simple and Complex Models. They introduced a linear model of a virtual refinery that simulated the blending of five types of RFG under more restrictive emission constraints for Petroleum Administration for Defense Districts I, II, and III. Moreover, they demonstrated a monotonic increase in the cost of emission reduction for the five considered RFGs and proved that the major cost component is accounted for by the cost of refinery operation and the cost of credit required to update the technology in response to the new standards. Nevertheless, ancillary costs (due to more difficult management of basic stocks, tests and traceability, non-compliant blends, safety margins, and operational flexibility) are considered to lie between 12% and 15% of the total blend cost. Online blend control would reduce these ancillary costs to a considerable extent.

^{Treiber (1998)} indicated that the strategies organized to satisfy the EPA-CM constraints were initially aimed at reducing reid vapor pressure (RVP) and later at reducing sulfur content (SUL). The sale of RFG violating the statutory emission thresholds would involve civil or criminal penalties. They also indicated the gasoline properties most heavily influencing emissions (i.e., VOC(x) = f(RVP), NO_x(x) = f(SUL), and TOX(x) = f(BEN)). Two optimization strategies were proposed. The first strategy consisted of defining adequate thresholds on RVP, SUL, and benzene content (BEN) and then using a linear program. To counter the interactions between variables and the effect of the other variables on emissions, a wide safety margin was observed. This optimization step resulted in a 1.1% blend cost reduction. The second strategy used a nonlinear optimizer to obtain all the variable interactions and precise emission values. The safety margin was scaled back, and the finer optimization approach lowered cost by 1.8%.

^{Naman (1999)} used linear regression to predict emissions according to gasoline quality. Ideally, a “good” linear approximation of emissions would solve the problem of integrating emission constraints during the programming phase, as well as during the gasoline blending phase. He established two models. The first model was local and resulted in randomly generating 100 blends that followed a normal distribution with a mean equal to a typical gasoline and standard deviation equal to the laboratory measurement error. In other words, he replicated the same gasoline 100 times. The standard deviations of his results were as follows: TOX= 1.162, NO_x= 10.419, and VOC= 33.539. The second model generated gasolines uniformly over a wide min/max interval (closer to reality). In the latter case, the linear approximation errors were as follows: TOX= 8.991, NO_x = 101.314, and VOC= 112.474. These standard deviations are extremely large, thereby reconfirming our opinion regarding the strong nonlinearity of functions used to calculate EPA-CM emissions.

In 2008 (given the jump in time, one could surmise that the problems concerning EPA-CM had been definitively resolved), ^{Furman and Androulakis (2008)} tackled the critical problem of satisfying even more restrictive emission constraints while imposing a minimum number of changes in the blending process. They asserted that the current EPA-CM representation is extremely difficult to implement in refining operations or combine with gasoline design models. The main motivation behind their work was to derive a model capable of including the EPA-CM constraints within a mathematical programming formulation. The major difficulty herein lies in the EPA-CM formulation as a series of complicated disjunctive constraints that basically makes its incorporation into algebraic optimization models impossible. At present, a large array of typical “boutique” gasolines are available (up to 14 in number) that meet state and/or local air quality specifications. Such formulations use the technique of writing disjunctive constraints as a series of mixed binary/continuous variables, resulting in an optimization problem with 20 binary variables, 100 continuous variables, and 200 mixed constraints (for a blending problem with 10 components and 12 properties). Their objective function is a normalized distance between the target gasoline x and the baseline gasoline xb (see objective function F_FA(x), Eq. (5)). The indicated calculation times are on the order of 10 min in a 2.4 GHz AMD Athlon processor with 2 GB memory.

Misener et al. dedicated a suite of papers to the problem of pooling. ^{Misener et al. (2010)} addressed the problem of extended pooling with emission constraints. They employed and improved the technique of ^{Furman and Androulakis (2008)} to write the disjunctive constraints as a series of mixed constraints (i.e., binary/continuous), covering all EPA-CM scenarios with the IF List. Nevertheless, among the illustrative examples given, only the first one actually imposes useful reductions on the blended gasoline. This “utility” is deduced with regard to the reference baseline emissions (see Table 4, page 1449 of this reference), which has been copied in Table 1:

As shown in Table 1, only the C1 case is nearly correct (except for TOX), whereas the C2 and C3 cases violate the threshold emissions. Computation times vary from a few seconds to 5274 s. Calculating the time ascribed to the integration of EPA-CM is difficult. For all these examples, the branch-and-bound procedure visits 44 nodes, for a total time of 6249 s, thereby amounting to 142 s/node on average. In other words, the CPLEX+ MINOS solver required on average 142 s to find the min/max bounds for any given node.

To conclude this bibliographical review, these references can be divided into two major categories, namely, i) economic aspects and ii) mathematical/modeling/implementation aspects. From the first category, the optimization with constraints on emissions exerts significant economic impact on the cost of gasoline. The second category underscores the difficulty of integrating the EPA-CM constraints into the refining operations or combining them with gasoline model design. Even when such an integration step is performed, the calculation time of MINLP models remains impracticable for use on the online control of emissions and blending.

Notably, these contributions do not compare their results because the approaches are probably different and incomparable.

EPA-CM and useful notations and equations

A target gasoline x (more specifically, x∈R^p is a vector of p gasoline properties) is blended from a set of components characterized by their properties B₁,…,B_n∈R^p. The following subset of these properties is used to calculate emissions in EPA-CM:

OXY= Oxygen content of (any) target fuel, expressed as weight percent;

SUL= Sulfur content of (any) target fuel, expressed in parts per million by weight;

RVP= Reid vapor pressure of the target fuel, expressed in pounds per square inch;

E200= 200°F distillation fraction of the target fuel, expressed as volume percent;

E300= 300°F distillation fraction of the target fuel, expressed as volume percent;

ARO= Aromatics content of the target fuel, expressed as volume percent;

BEN= Benzene content of the target fuel, expressed as volume percent;

OLE= Olefin content of the target fuel, expressed as volume percent

MTB= Methyl tertiary butyl ether content of the target fuel, expressed as weight percent;

ETB= Ethyl tertiary butyl ether content of the target fuel, expressed as weight percent;

TAM= Tertiary amyl methyl ether content of the target fuel, expressed as weight percent;

ETH= Ethanol content of the target fuel, expressed as weight percent.

Other properties such as RON and MON are calculated and controlled; however, the above list is limited to those contributing to the calculation of NO_x, VOC(x) and TOX(x). A recipe is a vector u∈Rⁿ, such that u_i, i=1,…,n is the volume percentage of component Bi in blend x. Recipe u must verify the constraint . The properties of target blend x are a function of the component matrix B=[B1,…,B_n] and recipe u, with x=x(B,u). Over a short time horizon (e.g., a few minutes, as is the case herein), B can be considered as a constant; consequently, the only control variable is recipe u. The hydraulic and operational constraints actually yield the upper and lower bounds on the recipes:

(1)\underset{\_}{u}\le u\le \bar{u},{{\displaystyle \Sigma }}_i{{\displaystyle u}}_i=\mathrm{1.}

The regulatory constraints on the target gasoline are given as a set of upper and lower bounds on its properties: \underline{c}\le x=x\left(u\right)\le \overline{c}. Moreover, x(u) should satisfy the EPA regulatory constraints (see Appendix). The intersection of these constraints yields

(2)\max (\underset{\_}{c},\underset{\_}{x})\le x=x(u)\le \min (\bar{c},\bar{x}).

In our control process, when gasoline x(u) violates the constraints in Eq. (2) or is not optimal, recipe u is modified to optimize some general form of an objective function: F\left(u\right)=\alpha ·\Vert u-u^0\Vert +\beta ·\Vert x\left(u\right)-x^0\Vert Here, u⁰ is a baseline recipe calculated by the scheduling department, and x⁰ is an “ideal” gasoline that minimizes the giveaway.

The emission calculation entails an additional step. (This step is given in detail for NO_x in Appendix.) All constraints in Eqs. (1) and (2) are linear, except for the one on RVP. In practice, some transformation of RVP that blends linearly is introduced. The functions used to calculate emissions in EPA-CM (i.e., NO_x(u), VOC(u) and TOX(u) are continuous but non-differentiable and non-convex (see Appendix). As shown in Appendix, these emissions are feasible if they verify constraints:

(3)\mathrm{N}\mathrm{O}x(x(u))\le \overline{N}\mathrm{V}\mathrm{O}\mathrm{C}(x(u))\le \overline{V}\mathrm{T}\mathrm{O}\mathrm{X}(x(u))\le \overline{T}.

A recipe u now becomes feasible if it verifies Constraints (1), (2), and (3).

The non-convex and non-differentiable functions NO_x(x(u)), VOC(x(u)) and TOX(x(u)) serve to input the constraints in Eq. (3); hence, ^{Furman and Androulakis (2008)}, ^{Misener et al. (2010)}, and all preceding authors attested that EPA-CM is extremely difficult to implement for refinery operations or combine with existing models in predicting the quality of blended gasoline. ^{Furman and Androulakis (2008)} and ^{Misener et al. (2010)} introduced mixed (i.e., binary/continuous) models to integrate EPA-CM into the quality prediction and optimization sequence for blended gasoline. These authors applied the technique of disjunctive constraints (^{Raman and Grossmann, 1994}) and inserted binary variables to manage the “IF List.” More specifically, they introduced mixed (i.e., binary/continuous) constraints to obtain the correct values for y (edge target) and z (delta target) of Step 1 (see Appendix). To control the NO_x(u), VOC(u) and TOX(u) emissions, ^{Furman and Androulakis (2008)} optimized the following objective function with respect to the constraints in Eqs. (1) and (2):

(4){{\displaystyle F}}_{FA}(x)={{\displaystyle {{\displaystyle \Sigma }}_i\left({\displaystyle \frac{x(i)-x\mathrm{b}(i)}{\overline{x}(i)-\underline{x}(i)}}\right)}}^{\mathrm{2}}.

As discussed in Section 5, this optimization tends to find a blend as near as possible to baseline blend xb (see Appendix). However, the optimization may not verify the constraints in Eq. (3).

Differential approximation of EPA-CM

We continue with the example of calculating NO_x(u) for a given recipe u as a means of explaining our approach. The VOC and TOX calculations are relatively similar. The calculation of terms found in TOX(u) (i.e., exhausted benzene, non-exhausted benzene, acetaldehyde, formaldehyde, and butadiene) is simple because the delta target z is not utilized.

The disjunctive constraint technique used by Furman and Androulakis (2008) and Misener et al. (2010) is recalled in the next section. This technique is illustrated by providing the equivalent form of the IF List (see Appendix). Afterward, we develop a new way of writing the IF List. Finally, a differential objective function is introduced to manage the emission levels.

Disjunctive formalism

Disjunctive programming was first developed in 1974 by ^{Balas (2010)}. Its MINLP model was built by ^{Raman and Grossmann (1994)}. ^{Furman and Androulakis (2008)} and ^{Misener and Floudas (2009)} extensively used the MINLP formalism of disjunctive programing in reformulating EPA-CM as an optimization problem. We examine the following conditions on the IF List:

If x(ARO)≥18 and x(ARO)≤36.8, then y(ARO)=x(ARO), z(ARO)=0.

If x(ARO)<18, then y(ARO)=18, z(ARO)= x(ARO)–18.

If x(ARO)>36.8, then y(ARO)=36.8, z(ARO)= x(ARO)–36.8.

The disjunctive constraints for these conditions are as follows (in this case, the x, y, z values correspond to the ARO property):

\left[\begin{array}{l}x\langle \mathrm{18}\\ y=\mathrm{18}\\ z=x-\mathrm{18}\end{array}\right]\vee \left[\begin{array}{l}x\rangle \mathrm{36.8}\\ y=\mathrm{36.8}\\ z=x-\mathrm{36.8}\end{array}\right]\vee \left[\begin{array}{l}\mathrm{18}\le x\le \mathrm{36.8}\\ y=x\\ z=\mathrm{0}\end{array}\right].

These conditions can now be written as a series of mixed linear constraints. The domain of valid properties for the EPA-CM functions is determined by the min/max bounds of x given by the constraints in Eq. (2).

^{Furman and Androulakis (2008)} introduced four binary variables b₀, b₁, b₂, b₃∈{0,1} to translate the disjunctive constraints on ARO. This example will be analyzed herein. The values x, y, z pertain to the ARO property.

General constraints on ARO: Variable b₀

This variable is introduced to yield the edge target y of x(ARO) used in all emission calculations.

(b0.1)\mathrm{y-18}\le \left(\bar{x}-\mathrm{18}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{0}}),

(b0.2)\mathrm{18-y}\le \left(\mathrm{18}-\underset{\_}{x}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{0}}).

We obtain b0=1 if and only if x(ARO)<18. When b0=1, the result of constraints in Equations (b0.1) and (b0.2) is y=18. If b0=0, then these two constraints simply state that the edge target y lies in the box [\underline{x},\overline{x}]. However, this finding is implied by Eq. (2) because the edge target interval is included in the interval [\underline{x},\overline{x}] for all the gasoline properties of EPA-CM.

Edge target and delta target constraints on ARO: Variable b₁; Case x(ARO)<18

These constraints are unique to the calculation of NO_x(u) . Here, b₁=1 if and only if x≤18. This equivalence is maintained by Constraint (b1.1). For b₁=1, y=18 from Constraints (b0.1) and (b0.2). Constraints (b1.2) and (b1.3) result in z=x–18. To simplify the notation, we introduce the min/max values for delta target\underset{\_}{z},\bar{z}.

(b1.1)x-\mathrm{18}\le \left(\bar{x}-\mathrm{18}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{1}}),

(b1.2)\mathrm{z-}x+\mathrm{18}\le \left(\bar{z}-\underset{\_}{x}+\mathrm{18}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{1}}),

(b1.3)x-\mathrm{18}-z\le \left(\bar{x}-\mathrm{18}-\underset{\_}{z}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{1}}).

If b1=0, Constraint (b1.1) simply states that x\le \bar{x}. This equation is a redundant statement and follows from Eq. (2). Constraints (b1.2) and (b1.3) are simplified to z-x\le \left(\bar{z}-\underset{\_}{x}\right) and \left(\bar{x}-\underset{\_}{z}\right), respectively. These constraints are a consequence of \underset{\_}{x}\le x\le \bar{x} and \underset{\_}{z}\le z\le \bar{z} and are thus redundant.

Edge target and delta target constraints on ARO: Variable ; Case 18≤x≤36.8

Here, b₂=1 if and only if 18≤x≤36.8.

(b2.1)\mathrm{18-}x\le \left(\mathrm{18}-\underset{\_}{x}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{2}}),

(b2.2)x-\mathrm{36.8}\le \left(\bar{x}-\mathrm{46}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{2}}),

(b2.3)y-x\le \left(\bar{x}-\underset{\_}{x}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{2}}),

(b2.4)x-y\le \left(\bar{x}-\underset{\_}{x}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{2}}),

(b2.5)z\le \left(\bar{z}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{2}}),

(b2.6)-z\le \left(-\bar{z}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{2}}).

If b₂=1 or equivalently if 18≤x≤36.8, then Constraints (b2.1) and (b2.2) imply that 18≤x≤36.8. Constraints (b2.3) and (b2.4) yield y=x. Moreover, Constraint (b2.5) with (b2.6) yields z=0.

Edge target and delta target constraints on ARO: Variable ; Case 36.8<x

Here, b₃=1 if and only if 36.8≤x. This equivalence is a consequence of Constraints (b3.1) and (b3.2).

(b3.1)\mathrm{36}\mathrm{.8-}x\le \left(\mathrm{36.8}-\underset{\_}{x}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{3}}),

(b3.2)y-\mathrm{36.8}\le \left(\bar{x}-\mathrm{36.8}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{3}}),

(b3.3)\mathrm{3}\mathrm{6}.\mathrm{8}-y\le \left(\mathrm{36.8}-\underset{\_}{x}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{3}}),

(b3.4)z-x+\mathrm{3}\mathrm{6}.\mathrm{8}\le \left(\bar{z}-\underset{\_}{x}+\mathrm{36.8}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{3}}),

(b3.5)x-\mathrm{3}\mathrm{6}.\mathrm{8}-z\le \left(\bar{x}-\mathrm{36.8}-\underset{\_}{z}\right)(\mathrm{1}-{{\displaystyle b}}_{\mathrm{3}}).

If b3=1, then Constraint (b3.3) yields y=36.8, and Constraint (b3.4) with (b3.5) gives z=x–36.8. When b3=0, these constraints are automatically verified (to be redundant) as a consequence of Eq. (2) and \underset{\_}{z}\le z\le \bar{z}. The following constraint is then added to limit the number of possible cases to just one:

(b4){{\displaystyle b}}_{\mathrm{1}}+{{\displaystyle b}}_{\mathrm{2}}+{{\displaystyle b}}_{\mathrm{3}}\ge \mathrm{1.}

Consequently, to ensure the correct case for the edge target and delta target in the NO_x(u) calculation, ^{Furman and Androulakis (2008)} introduced the binary variables b₀, b₁, b₂, and b₃, along with 17 mixed constraints: (b0.*), (b1.*), (b2.*), (b3.*), and (b4).

All the disjunctive constraints introduced by ^{Furman and Androulakis (2008)} require 20 binary variables, 100 continuous variables, and 200 mixed constraints. ^{Misener et al. (2010)} achieved this result using 15 binary variables, 56 continuous variables, and 100 mixed constraints. A common characteristic of these approaches is the large number of redundant constraints. As shown in the ARO example, as their values vary, only a small subset of constraints will be of any utility, whereas the others will be useless for the solver. Our approach does not make use of binary variables and has only a few redundant constraints or none at all, as discussed in the next section.

Functional form of the IF List

Our approach is based on the simple fact that the binary condition t≤r is equivalent to S(t, r)=0, where S(t, r) is the function

S(t,r)=\{\begin{array}{l}\mathrm{0}\hspace{0.17em}t\le r\\ \mathrm{1}\hspace{0.17em}t\rangle r\end{array}.

One implementation of this function might use the sign function:

sign(t-r)=\{\begin{array}{l}-\mathrm{1}\hspace{0.17em}t\le r\\ \mathrm{1}\hspace{0.17em}t\rangle r\end{array},

(5)S(t,r)=\mathrm{0.5}(\mathrm{1}+\mathrm{sign}(t-r)).

S(t, r) has been used to calculate the correct edge target and delta target values for x(ARO) given in the previous section. We simply set m=18 and M=36.8. Then,

(6)y=(1-S(x,m))\cdot m+S(x,m)\cdot (\mathrm{1}-S(x,M))\cdot x+S(x,M)\cdot M,

(7)\mathrm{z=}(\mathrm{1}-S(x,m))\cdot (x-m)+S(x,M)\cdot (x-M).

The equivalence of these expressions with the result of the IF List is straightforward.

Case x<m. S(x,m)=0 and S(x,M)=0, which gives y=m and z=x–m.

Case x≥m and x≤M. S(x,m)=1 and S(x,M)=0, which gives y=x and z=0.

Case x>M. S(x,m)=1 and S(x,M)=1, which gives y=M and z=x–M.

Consequently, Eqs. (2), (6), and (7) are all equivalent to the set of mixed constraints (b*.*) discussed in Section 4.1. The correct values of y and z can thus be calculated online, which simplifies Step 1 of the emission calculation (see Appendix). Steps 2 and 3 are merely composite functions with Eqs. (6) and (7). Intrinsically, no additional constraints are needed to find the correct edge target and delta target values.

We now introduce a differentiable approximation of S(t, r) via the sigmoid.

(8)SC(t,r,a)=\mathrm{0.5}(\mathrm{1}+\tanh (a(t-r))).

The coefficient a={\displaystyle \frac{\mathrm{d}(SC(t,r,a))}{\mathrm{d}t}} for t=r can be selected, such that \left|SC(t,r,a)-S(t,r)\right|\langle \epsilon for any given ε everywhere except within a small interval t\in [r-\delta ,r+\delta ]. We consider the level of precision when measuring any property t of gasoline x in our calculations. Selecting a, such that 2δ is less than the measurement error for any property of the blended target gasoline required for the emission calculation, is in fact straightforward. The deviation in emissions calculated using Eq. (8) compared to those calculated with Eq. (5) is controlled by the level of measurement accuracy of the properties of x.

Gradient of emission functions

We now calculate the gradient of NO_x(u) versus recipe u and gasoline properties x=x(u). This process is a simple calculus exercise by following the chain rule for derivative of composite functions. Nevertheless, we indicate it here for the sake of having a self-contained paper. The gradient calculation for VOC(u) and TOX(u) follows a similar manner.

{\displaystyle \frac{\partial NOx(u)}{{{\displaystyle \partial u}}_k}}={{\displaystyle \Sigma }}_{j=\mathrm{1}}^p{\displaystyle \frac{\partial NOx(x)}{{{\displaystyle \partial x}}_j}}\cdot {\displaystyle \frac{{{\displaystyle \partial x}}_j}{{{\displaystyle \partial u}}_k}},k=\mathrm{1},\mathrm{...},n.

The partial derivative {\displaystyle \frac{{{\displaystyle \partial x}}_j}{{{\displaystyle \partial u}}_k}} for each property xj as a function of recipe u is required. For the linear blending law properties, this partial derivative is simply the entry B(j,k) of the component matrix. Edge target y and delta target z are both defined as functions of x; consequently,

{\displaystyle \frac{\partial NOx(x)}{{{\displaystyle \partial x}}_i}}={{\displaystyle \Sigma }}_{j=\mathrm{1}}^p\left({\displaystyle \frac{\partial NOx(y,z)}{{{\displaystyle \partial y}}_j}}\cdot {\displaystyle \frac{{{\displaystyle \partial y}}_j}{{{\displaystyle \partial x}}_i}}+{\displaystyle \frac{\partial NOx(y,z)}{{{\displaystyle \partial z}}_j}}\cdot {\displaystyle \frac{{{\displaystyle \partial z}}_j}{{{\displaystyle \partial x}}_i}}\right).

Steps 2 and 3 (Appendix) of the calculation of each term of NOx(y,z)={{\displaystyle w}}_N\cdot N(y)\cdot FN(y,z)+{{\displaystyle w}}_H\cdot H(y)\cdot FN(y,z)) can now be recovered.

{\displaystyle \frac{\partial NOx(y,z)}{{{\displaystyle \partial y}}_i}}={{\displaystyle w}}_N\cdot \left({\displaystyle \frac{\partial NO(y)}{{{\displaystyle \partial y}}_j}}\cdot FN(y,z)+{\displaystyle \frac{\partial FN(y,z)}{{{\displaystyle \partial y}}_j}}\cdot N(y)\right)+{{\displaystyle w}}_H\left({\displaystyle \frac{\partial H(y)}{{{\displaystyle \partial y}}_j}}\cdot FH(y,z)+{\displaystyle \frac{\partial FH(y,z)}{{{\displaystyle \partial y}}_j}}\cdot H(y)\right).

Given the exponential form of N(y) and H(y), we obtain

{\displaystyle \frac{\partial NOx(y,z)}{{{\displaystyle \partial y}}_i}}={{\displaystyle w}}_N\cdot N(y)\cdot \left({\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{{{\displaystyle \partial y}}_j}}\cdot FN(y,z)+{\displaystyle \frac{\partial FN(y,z)}{{{\displaystyle \partial y}}_j}}\right)+{{\displaystyle w}}_H\cdot H(y)\left({\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{{{\displaystyle \partial y}}_j}}\cdot FH(y,z)+{\displaystyle \frac{\partial FH(y,z)}{{{\displaystyle \partial y}}_j}}\right),

{\displaystyle \frac{\partial NOx(y,z)}{{{\displaystyle \partial z}}_i}}={{\displaystyle w}}_N\cdot {\displaystyle \frac{\partial FN(y,z)}{{{\displaystyle \partial z}}_j}}\cdot N(y)+{{\displaystyle w}}_H\cdot {\displaystyle \frac{\partial FH(y,z)}{{{\displaystyle \partial z}}_j}}\cdot H(y).

For the next step, the partial derivatives on y and z are calculated. To avoid cumbersome notations, the following sections denote the coordinates of vectors y and z by y and z indexed by gasoline properties. For example, y:=y(OXY) indicates that y is defined as the y(OXY) value of edge target y.

y:=y(OXY),Z:=Z(OXY).

The following numbers are found in (EPA, 2016):

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=\mathrm{0.0018571},

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=-\mathrm{0.00913},

{\displaystyle \frac{\partial FN(y,z)}{\partial y}}={\displaystyle \frac{\partial FH(y,z)}{\partial y}}={\displaystyle \frac{\partial FN(y,z)}{\partial z}}={\displaystyle \frac{\partial FH(y,z)}{\partial z}}=\mathrm{0},

y:=y(SUL),Z:=Z(SUL).

We find

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=\mathrm{0.00069205}-\mathrm{2}\cdot \mathrm{6.6263}\cdot {\mathrm{10}}^{-\mathrm{7}}\cdot y,

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=\mathrm{0.000252},

{\displaystyle \frac{\partial FN(y,z)}{\partial y}}=z\cdot (-\mathrm{0.00000133}),

{\displaystyle \frac{\partial FN(y,z)}{\partial y}}=\mathrm{0},

{\displaystyle \frac{\partial FN(y,z)}{\partial z}}=(-\mathrm{0.00000133}\cdot y+\mathrm{0.000692}),

{\displaystyle \frac{\partial FN(y,z)}{\partial z}}=\mathrm{0.000252},

y:=y(RVP),z:=z(RVP).

We also find

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=\mathrm{0.0090744},

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=-\mathrm{0.013973},

{\displaystyle \frac{\partial FN(y,z)}{\partial y}}={\displaystyle \frac{\partial FH(y,z)}{\partial y}}={\displaystyle \frac{\partial FN(y,z)}{\partial z}}={\displaystyle \frac{\partial FH(y,z)}{\partial y}}=\mathrm{0},

y:=y(E\mathrm{200}),z:=z(E\mathrm{200}),

which yields

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=\mathrm{0.00093065},

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=-\mathrm{0.000931},

{\displaystyle \frac{\partial Fn(y,z)}{\partial y}}={\displaystyle \frac{\partial Fh(y,z)}{\partial y}}={\displaystyle \frac{\partial Fn(y,z)}{\partial z}}={\displaystyle \frac{\partial Fh(y,z)}{\partial y}}=\mathrm{0},

y:=y(E\mathrm{300}),z:=z(E\mathrm{300}),

which yields

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=\mathrm{0.00084596},

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=-\mathrm{0.004009},

{\displaystyle \frac{\partial Fn(y,z)}{\partial y}}={\displaystyle \frac{\partial Fh(y,z)}{\partial y}}={\displaystyle \frac{\partial Fn(y,z)}{\partial z}}={\displaystyle \frac{\partial Fh(y,z)}{\partial y}}=\mathrm{0},

y:=y(ARO),z:=z(ARO).

We now have

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=\mathrm{0.0083632}-\mathrm{2}\cdot \mathrm{0.00011905}\cdot y,

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=\mathrm{0.007097005}-\mathrm{2}\cdot \mathrm{0.000079951}\cdot y,

{\displaystyle \frac{\partial FN(y,z)}{\partial y}}=z\cdot (-\mathrm{0.000238}),

{\displaystyle \frac{\partial FH(y,z)}{\partial y}}=z\cdot (-\mathrm{0.0001599}),

{\displaystyle \frac{\partial FH(y,z)}{\partial z}}=(-\mathrm{0.000238}\cdot y+\mathrm{0.0083632}),

{\displaystyle \frac{\partial FH(y,z)}{\partial z}}=-\mathrm{0.0001599}\cdot y+\mathrm{0.007097},

y:=y(OLE),z:=z(OLE),

which leaves us with

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{1}}(y)}{\partial y}}=-\mathrm{0.0027735}+\mathrm{2}\cdot \mathrm{0.00036652}\cdot y,

{\displaystyle \frac{{{\displaystyle \partial n}}_{\mathrm{2}}(y)}{\partial y}}=-\mathrm{0.0027603}+\mathrm{2}\cdot \mathrm{0.000366}\cdot y,

{\displaystyle \frac{\partial FN(y,z)}{\partial y}}=z\cdot (\mathrm{0.000733}),

{\displaystyle \frac{\partial FH(y,z)}{\partial y}}=z\cdot (\mathrm{0.000732}),

{\displaystyle \frac{\partial FN(y,z)}{\partial z}}=\mathrm{0.000733}\cdot y-\mathrm{0.002774},

{\displaystyle \frac{\partial FH(y,z)}{\partial z}}=\mathrm{0.000730}\cdot y-\mathrm{0.00276.}

For any property p different from SUL, OLE, and ARO, we obtain y(p)=x(p) and z(p)=0. For these properties, {\displaystyle \frac{\partial y(p)}{\partial x(p)}}=\mathrm{1} and {\displaystyle \frac{\partial z(p)}{\partial x(p)}}=\mathrm{0}. For the SUL, OLE, and ARO properties, the calculation of these derivatives requires computing the derivative SC'(t, r, a) of SC(t, r, a). For each of these properties, we assume

p=SUL, m=10, M=450, a=prec(SUL) (with the SUL precision given by the user);

p=OLE, m=3.77, M=19, a=prec(OLE) (with the OLE precision given by the user);

p=ARO, m=18, M=36.8, a=prec(ARO) (with the ARO precision given by the user);and then calculate

{\displaystyle \frac{{{\displaystyle \partial y}}_p}{{{\displaystyle \partial x}}_p}}=-m\cdot {SC}^{\text{'}}({{\displaystyle x}}_p,m,a)+M\cdot {SC}^{\text{'}}({{\displaystyle x}}_p,M,a)+{SC}^{\text{'}}({{\displaystyle x}}_p,m,a)\cdot (\mathrm{1}-SC({{\displaystyle x}}_p,M))\cdot {{\displaystyle x}}_p-SC({{\displaystyle x}}_p,m)\cdot {SC}^{\text{'}}({{\displaystyle x}}_p,M,a)\cdot {{\displaystyle x}}_p+S({{\displaystyle x}}_p,m)\cdot (\mathrm{1}-SC({{\displaystyle x}}_p,M)),

{\displaystyle \frac{{{\displaystyle \partial z}}_p}{{{\displaystyle \partial x}}_p}}=-({{\displaystyle x}}_p-m)\cdot {SC}^{\text{'}}({{\displaystyle x}}_p,m,a)+(\mathrm{1}-SC({{\displaystyle x}}_p,m))+SC({{\displaystyle x}}_{\mathrm{2}},M)+({{\displaystyle x}}_p-M)\cdot {SC}^{\text{'}}({{\displaystyle x}}_p,M,a).

This expression represents the last link in the chain rule for derivatives. These analytical gradients are extremely easy to calculate as functions of recipe u and properties x. This crucial point serves to control emission variations online the same way as the other traditional gasoline properties.

Objective function to control emissions

We have now defined constraints on recipe u and properties x. We have also defined constraints on emissions in Eq. (3). Given a vector of positive penalties \pi =[{{\displaystyle \pi }}_N,{{\displaystyle \pi }}_V,{{\displaystyle \pi }}_T], the previous constraints in Eq. (3) are satisfied if and only if a recipe u exists, such that

F(u,\pi )={{\displaystyle \Sigma }}_{R\in \{N,V,T\}}{{\displaystyle \pi }}_R\cdot SG\left(R,\bar{R}\right)\cdot \left(R(u)-\bar{R}\right)=\mathrm{0.}

The idea behind this condition is to minimize the function F(u, π) with respect to the linear constraints on recipes in Eq. (1) and properties in Eq. (2). In this manner, we are not changing the online control framework for traditional gasoline properties but merely add a new objective function.

Blend instance

For a given component matrix B, a blend instance is defined by the constraints in Eqs. (1) and (2) and the penalty vector π.

Feasible blend instance

A blend instance is feasible if and only if the optimization problem

(9)\{\begin{array}{l}\min F(u,\pi )={{\displaystyle \Sigma }}_{R\in \{N,V,T\}}{{\displaystyle \pi }}_R\cdot SG\left(R,\bar{R}\right)\cdot \left(R(u)-\bar{R}\right)\\ \underset{\_}{u\le u\le \bar{u}}\\ \max \left(\underset{\_}{x},\underset{\_}{c}\right)\le x(u)\le \min \left(\bar{x},\bar{c}\right)\end{array}

has a feasible (optimal) solution u₀, such that F(u₀, π)=0.

In the following, we denote the cost function as F₁(u)=c^Tu and the giveaway function as {{\displaystyle F}}_{\mathrm{2}}(u)=\Vert x(u)-x^{\circ }\Vert. The optimization problem in Eq. (9) is central to our online emission control system. The resolution of this problem has been implemented in MATLAB (see Section 5). The general framework for optimization and online emission control is given below. In implementing the instances of this framework, Steps 2 and 3 can be combined.

Online optimization and control process:

Step 1: Define a blend instance

Choose a component matrix B. Define the constraints in Eqs. (1) and (2). Choose the penalty vector π>0.

Step 2: Feasibility

Solve the optimization problem in Eq. (9). Let u₀ be its optimal solution.

If F({{\displaystyle u}}_{\mathrm{0}},\pi )\rangle \mathrm{0} // Emission threshold violated

Change the blend instance; GO TO Step 1;

Otherwise

GO TO Step 3 to identify an optimal recipe.

Step 3: Optimization

3.1. Compute the values of F₁(u₀) and F₂(u₀).

3.2. Find the optimal solution (S*, T*, u*) to the following problem:

(10)\{\begin{array}{l}\begin{array}{l}\min G(S,T,u)=a\cdot S+b\cdot T+w\cdot F(u,\pi )\\ {{\displaystyle F}}_{\mathrm{1}}(u)\le S\cdot {{\displaystyle F}}_{\mathrm{1}}({{\displaystyle u}}_{\mathrm{0}})\\ {{\displaystyle F}}_{\mathrm{2}}(u)\le T\cdot {{\displaystyle F}}_{\mathrm{2}}({{\displaystyle u}}_{\mathrm{0}})\end{array}\\ \mathrm{0}\le S,T\le \mathrm{1}\\ \underset{\_}{u}\le u\le \bar{u}\\ \max \left(\underset{\_}{x},\underset{\_}{c}\right)\le x\le \min \left(\bar{x},\bar{c}\right)\end{array}.

3.3 IfF(u^{*},\pi )=\mathrm{0}

STOP: The optimal recipe has been found;

Otherwise

Find a recipe {{\displaystyle u}}_{\mathrm{1}}\in [{{\displaystyle u}}_{\mathrm{0}},u^{*}], such that F(u₁, π)=0;

{{\displaystyle u}}_{\mathrm{0}}\leftarrow {{\displaystyle u}}_{\mathrm{1}}; GO TO Step 3.1.

In the following section, we present the results from an actual refinery.

Results

The quick calculation of the gradient of F(u, π) is extremely important for all nonlinear optimization solvers. The optimization problem in Eq. (9) can actually feature up to 80 linear constraints. A standard branch-and-bound global optimizer is used to solve Eq. (9) given that the objective function is non-convex. This branching occurs within the domain of recipe u in verifying 0≤u≤1 and Σu_i=1. In practice, the dimension of u (i.e., the number of components) does not exceed 20 in typical applications. These recipe domain properties make allows the inclusion of a branch-and-bound procedure stop criterion that depends on the diameter of cells created by the branching procedure. Hence, reducing the diameter of the box containing the optimal recipe to a value lying below a given threshold (10^-8 for our application) is rather easy. All computations are performed on a machine with a Duo CPU, 2 GHz Intel processor, and 3 GB memory. All our optimization instances are run in less than 1 s. These results serve to implement the online emission control strategy into the TOTAL ANAMEL software (ANAMEL is the software used in several TOTAL Group refineries to control blended product quality by means of real-time blend recipe optimization). Preliminary tests with these new emission control features embedded into ANAMEL and approximate values for distillation figures are conducted at the Leuna refinery (Germany) in 2015.

For comparison of results, the optimal values and the emission performance of the following three objective functions are calculated:

(i) Cost function F₁(u)=c^Tu;

(ii) Giveaway function as the distance {{\displaystyle F}}_{\mathrm{2}}(u)=\Vert x(u)-x^{\circ }\Vert to an ideal gasoline x⁰;

(iii) Normalized distance to the baseline gasoline F_FA of Furman and Androulakis (2008) in Eq. (4).

We now provide a typical example with data from an American refinery. The EPA-CM parameters are as follows in Table 2.

We then assign arbitrary penalties to illustrate our approach. Table 3 recalls the baseline RFG properties for the given season and region:

The target blend is 57,000 m³ of RFG. The properties of six components (B₁, ..., B₆) are listed in Table 4.

The two bottom rows list the lower and upper bounds on recipes. The two columns on the right describe the target zone. Table 5 presents the following optimization results for all objective functions implemented: i) F₁: optimal value of the economic function (i.e., blend cost); ii) F₂: optimal value of the objective (ii); iii) F₃: optimal value of the F_FA function in Eq. (4) optimized by ^{Furman and Androulakis (2008)}; and iv) F₄: optimal value of the objective function in Eq. (9) (emission control). The optimal recipes of each optimization are given in columns u₁, u₂, u₃ and u₄. The corresponding RFG gasolines appear in columns x₁, x₂, x₃ and x₄. The emissions of these gasolines are shown in columns E₁, E₂, E₃, and E₄. The cost or giveaway optimizations lead to blends with high emission levels. The optimization of F₃ improves emission levels but still violates the VOC emission threshold. For this blend instance, only the optimization of F₄ yields the gasoline x₄ with emissions E₄, which lie far below the emission thresholds. The computation time for this instance is 0.6 s.

The values in bold indicate the results that satisfy the emissions thresholds.

The comparison of execution times is performed with respect to the approach of ^{Furman and Androulakis (2008)}. We do not have data to implement the pooling problem of ^{Misener and Floudas (2009)}.

Conclusions

A differentiable approximation of EPA-CM is developed in this study. This approximation is based on a computation chain that offers ease of implementation on online emission control. Only a new objective function should be added to our general control framework for traditional gasoline properties. Our model has been implemented as a global optimization problem of a differentiable objective function with respect to linear constraints. This approach is relatively straightforward in refinery operations. The execution time of our model is several orders of magnitude faster than the MINLP models developed by ^{Misener and Floudas (2009)} and ^{Furman and Androulakis (2008)}. Our typical actual example has demonstrated the simplicity of our model for emission computations. This model is now patented (^{Janaqi et al., 2015}) and implemented using our ANAMEL online controller.

References

[1] Balas E (2010). Disjunctive programming. In: 50 Years of Integer Programming 1958–2008: From the Early Years to the State-of-the-Art, 283–340

[2] Cason W W (1997). The illusory flexibility of the complex model: a graphical analysis of the implications of the complex model for refinery blending flexibility. In: NPRA Annual Meeting. San Antonio, Texas: NPRA Annual Meeting, 16–18

[3] Chèbre M, Creff Y, Petit N (2010). Feedback control and optimization for the production of commercial fuels by blending. Journal of Process Control, 20(4): 441–451

[4] EPA (2015). Complex emission model.

[5] EPA (2016). Template spreadsheet for calculating emissions from combustion plants.

[6] Furman K C, Androulakis I P (2008). A novel MINLP-based representation of the original complex model for predicting gasoline emissions. Computers & Chemical Engineering, 32(12): 2857–2876

[7] Hirshfield D S, Kolb J A (1994). Minimize the cost of producing reformulated gasoline by integrating EPA’s Complex Model into refinery linear programming (LP) models. NPRA Annual Meeting, 64–67

[8] Hirshfield D S, Kolb J A (1997). The economics of gasoline reformulation: refining economics, emissions standards, and the Complex Model. In: NPRA Annual Meeting, San Antonio, USA

[9] Jackson J, Vittachi K (1995). CITGO Corpus Christi refinery gasoline blender upgrade. In: NPRA Annual Meeting, San Francisco, USA

[10] Janaqi S, Chebre M, Pitollat G (2015). Method and device for monitoring induced properties of a mixture of components, in particular emission properties, 69.

[11] Korotney D J (1995). Reformulated gasoline effects on exhaust emissions: Phase III; investigation on the effects of sulfur, olefins, volatility, and aromatics and the interactions between olefins and volatility or sulfur. SAE Special Publications, 179–186

[12] Misener R, Floudas C A (2009). Advances for the pooling problem: modeling, global optimization, and computational studies Survey. Applied and Computational Mathematics, 8(1): 3–22

[13] Misener R, Gounaris C E, Floudas C A (2010). Mathematical modeling and global optimization of large-scale extended pooling problems with the (EPA) complex emissions constraints. Computers & Chemical Engineering, 34(9): 1432–1456

[14] Naman B T (1999). Linear models help refiners develop RFG [(reformulated gasoline)] recipes. Oil & Gas Journal, 97(8): 64–66

[15] Nocca J L, Forestiere A, Cosyns J (1994). Diversify process strategies for reformulated gasoline. Fuel Reformulation, 4(5): 18–22

[16] Raman R, Grossmann I E (1994). Modelling and computational techniques for logic based integer programming. Computers & Chemical Engineering, 18(7): 563–578

[17] Treiber S (1998). RFG [(reformulated gasoline)]: the challenge to conventional blending technology. In: Gulf Publishing, eds. Hydrocarbon Processing 2nd International Conference on Process Optimization. Houston: Hydrocarbon Processing 2nd International Conference on Process Optimization, 101–103

[18] Trierwiler L D (1995). Representing the simple and complex models in linear programming with respect to reformulated gasoline. In: American Energy Week ’95 “Pipelines, Terminals Storage, and Reformulated Fuels” Conference’ Houston, USA

RIGHTS & PERMISSIONS

The Author(s) 2017. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)