Engineering 3 (2017) 188–201

Contents lists available at ScienceDirect

Engineering
j o u r n a l h o m e p a g e : w w w . e l s e v i e r. c o m / l o c a t e / e n g

Research
Smart Process Manufacturing—Article

Global Optimization of Nonlinear Blend-Scheduling Problems

Pedro A. Castillo Castillo a, Pedro M. Castro b, Vladimir Mahalec a,*
a
Department of Chemical Engineering, McMaster University, Hamilton, ON L8S 4L8, Canada
b
Center for Mathematics, Fundamental Applications and Operations Research, Faculty of Sciences, University of Lisbon, Lisbon 1749-016, Portugal

a r t i c l e i n f o a b s t r a c t

Article history: The scheduling of gasoline-blending operations is an important problem in the oil refining industry. This
Received 7 December 2016 problem not only exhibits the combinatorial nature that is intrinsic to scheduling problems, but also
Revised 16 February 2017 non-convex nonlinear behavior, due to the blending of various materials with different quality properties.
Accepted 20 February 2017
In this work, a global optimization algorithm is proposed to solve a previously published continuous-time
Available online 28 March 2017
mixed-integer nonlinear scheduling model for gasoline blending. The model includes blend recipe optimi-
zation, the distribution problem, and several important operational features and constraints. The algorithm
Keywords: employs piecewise McCormick relaxation (PMCR) and normalized multiparametric disaggregation tech-
Global optimization
nique (NMDT) to compute estimates of the global optimum. These techniques partition the domain of one
Nonlinear gasoline blending
of the variables in a bilinear term and generate convex relaxations for each partition. By increasing the num-
Continuous-time scheduling model
Piecewise linear relaxations ber of partitions and reducing the domain of the variables, the algorithm is able to refine the estimates of
the global solution. The algorithm is compared to two commercial global solvers and two heuristic methods
by solving four examples from the literature. Results show that the proposed global optimization algorithm
performs on par with commercial solvers but is not as fast as heuristic approaches.
© 2017 THE AUTHORS. Published by Elsevier LTD on behalf of the Chinese Academy of Engineering and
Higher Education Press Limited Company. This is an open access article under the CC BY-NC-ND
license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction imize or minimize a desired objective such as profit, total cost, lead
time, and so forth. Scheduling software and tools based on mathe-
Computing and implementing an optimal production schedule matical programming is becoming more usual in practice.
can reduce operational costs, increase profit margins, and avoid The scheduling of gasoline-blending operations is an important
deviations from environmental constraints [1]. However, complex and relevant industrial problem because gasoline accounts for 60%–
industrial plants can have multiple production, storage, and distri- 70% of the total profit of an oil refinery [2–4]. In a gasoline-blending
bution subsystems, several distinct raw materials and intermediate system, components from dedicated supply tanks are mixed in
and final products, and intricate connections between all these ele- blending tanks or in-line blenders and sent to product tanks. Blend-
ments that make scheduling a difficult decision-making process. ing tanks can either have flow in or out (Fig. 2), resembling batch
Scheduling problems typically deal with four main decisions [1]: operation. In contrast, in-line blenders operate in a continuous
① determining the required tasks to fulfill the corresponding ob- manner (Fig. 3). In addition to the four decisions mentioned earlier,
jectives, requirements, and/or demand targets; ② assigning each scheduling blend operations should also involve determining the
task to a processing unit or resource that is available in the network; blend recipes—that is, the amounts of components to mix such that
③ defining the sequence in which the tasks will be executed; and products’ quality properties meet given specifications.
④ timing the tasks—that is, determining when to start and stop The gasoline-blending system studied in this work is described
each one (Fig. 1). Optimal scheduling decisions are those that max- in Fig. 3. Gasoline blending is carried out by one or more continuous

* Corresponding author.
E-mail address: mahalec@mcmaster.ca

http://dx.doi.org/10.1016/J.ENG.2017.02.005
2095-8099/© 2017 THE AUTHORS. Published by Elsevier LTD on behalf of the Chinese Academy of Engineering and Higher Education Press Limited Company.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 189

Fig. 1. Main scheduling decisions.

Fig. 2. Batch-blending system. The variable t means time period.

quality giveaways [4,5].

Since the 1960s, there has been a significant effort to derive so-
called “blend indices.” These are nonlinear transformations of the
actual quality properties of the blend components, which then, as a
linear combination, can predict the quality of the blend. Even with
this approach, two issues remain:
(1) If a product is blended into a tank, there is always some ma-
terial in the tank (the so-called “tank heel”) left over from
previous blends. Any new blend must include the properties
of the material in the tank in the calculation of the new blend
recipe, which leads to a nonlinear blending model for multi-
period scheduling, even when using blend indices.
(2) Blend properties calculated from these indices are not 100%
accurate, and a cumulative annual quality giveaway of, for ex-
ample, octane can amount to very large cost increments. This
forces the use of nonlinear blending models for the calcula-
Fig. 3. General scheme of the continuous gasoline-blending system studied in this
tion of, for example, octane numbers.
work.
Mathematical models for scheduling problems are usually for-
mulated as mixed-integer linear programming (MILP). However, for
blenders. Each blender is connected to the sources of blend compo- gasoline blending, nonlinear behavior is intrinsic to the correspond-
nents. Blended material in some refinery configurations goes to a ing process and mixed-integer nonlinear programming (MINLP)
storage tank, while in other configurations, it can go directly to the needs to be employed for the sake of accuracy. Most nonlinear terms
pipeline. Since there are several grades of gasoline produced (e.g., are non-convex, making convex optimization techniques ineffective.
regular, medium, premium), the blender switches from blending A global optimization approach is thus required. Before describing
one grade to another. Each switch requires (partial) realignment of the proposed global optimization method, a brief review of previous
blend feeding lines, which leads to a loss of blending capacity. In work is presented in the following paragraphs.
addition, switching to a different range of quality properties often
requires resetting or recalibration of the analyzers used to measure 1.1. Literature review on refinery scheduling
them.
Some of the most important quality properties of gasoline are Scheduling models can be divided into two main categories based
research octane number (RON), motor octane number (MON), Reid on the treatment of the time domain: discrete- and continuous-
vapor pressure (RVP), density, sulfur, aromatics (Ar), and olefins (Ol) time formulations. In discrete-time models, the time horizon is di-
content. RON, MON, and RVP do not blend linearly; thus, considering vided into several time periods of known duration with fixed start
nonlinear blending rules for such quality properties in the sched- and end time. In continuous-time models, the time horizon is par-
uling model can increase the accuracy of the solution and reduce titioned into time slots whose duration will be determined by the
190 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201

optimization. While continuous-time models generate problems Lotero et al. [17] proposed another formulation of the multi-period
with fewer discrete variables than their discrete counterparts, they pooling problem. They denominated this discrete-time MINLP for-
are more complex to formulate and often feature many “big-M” mulation as a hybrid of a source-based model (similar to Castro’s
constraints that, due to their weak relaxations [6], compromise split-fraction model [18]) and a concentration-based model [19].
computational performance. More in-depth reviews of scheduling Redundant constraints were added to improve the linear relaxation,
formulations can be found in Refs. [1,7–9]. and the model was solved using a two-stage MILP-NLP approach.
Gasoline blending has been studied by many researchers due to The MILP was a relaxation of the original MINLP model. The NLP
its commercial importance and non-convex features, which makes model was obtained by fixing the integer variables of the original
it a suitable subject for testing different formulations and algorith- model to the values computed by the MILP. The algorithm adds
mic approaches. Operational constraints found in gasoline blending integer, optimality, or feasibility cuts to the MILP model at each it-
are related to the presence of multipurpose tanks and non-identical eration, and stops when the difference between the MILP and NLP
blenders, to different storage-tank policies (e.g., whether the simul- solutions is smaller than a pre-specified tolerance.
taneous receipt and delivery of material is allowed or not), and to Cerdá et al. [20] presented a continuous-time MILP formulation
practical aspects such as minimum blend sizes and minimum blend- that uses floating slots dynamically allocated to time periods while
er running and setup times. Not all of these constraints are consid- solving the problem. The model included most of the operation-
ered in published scheduling models. In some cases, blend recipes al constraints found in practice. Cerdá et al. [21] then extended
are assumed to be fixed (i.e., they cannot be optimized). The down- the model to handle nonlinear blending rules, thus formulating a
stream distribution or shipping problem (i.e., timing delivery tasks continuous-time MINLP model. An approximate MILP formulation
to fulfill the demand) is sometimes also part of the blend-scheduling was derived by replacing the nonlinear blending rules with linear
problem. blend indices. The values of the binary variables computed by this
Méndez et al. [2] presented both a discrete- and continuous-time MILP were fixed in the original MINLP, thus becoming an NLP that
MILP model to schedule gasoline-blending operations. An iterative was solved to find a near-optimal solution of the original problem.
method was employed to handle nonlinear blending rules while
preserving the linearity of the models. Several key operational con- 1.2. Literature review on global optimization
straints were omitted and the distribution problem was not consid-
ered. Global optimization of nonlinear non-convex problems has been
Jia and Ierapetritou [10] developed a continuous-time MILP mod- a subject of extensive research over the last three decades. Even
el to simultaneously schedule gasoline-blending tasks and distribu- though powerful commercial solvers have been developed [22,23],
tion operations. The linearity of the model was maintained by using there has been a continuous stream of advances in the field.
given preferred recipes. Their model was later extended to schedule Global optimization algorithms have in common the genera-
operations of the main processing units in an oil refinery [11]. tion of a convex relaxation of the problem, which provides a lower
Glismann and Gruhn [12] used a two-level approach based on bound to the value of the objective function, and a way to generate
discrete-time models. Blend recipes and production targets were feasible solutions (the upper bound). In addition, they have a meth-
computed first using a nonlinear programming (NLP) model. Then a od of bringing the lower and upper bounds together, so that the op-
MILP model was employed to solve the short-term scheduling prob- timality gap can be reduced to ε-tolerance.
lem using such recipes and targets. The scheduling model was based Computing a tight lower bound is absolutely critical. This may
on the resource-task-network (RTN) representation and did not con- involve replacing the original formulation with an equivalent that
sider multipurpose tanks or the delivery-scheduling problem. preserves the feasible space but has a stronger relaxation (the differ-
Li et al. [13] formulated a continuous-time MILP model featuring ent formulations for the pooling problem are a well-known example
a common time grid for all units (i.e., blenders and tanks). Li and [24]). Another option is to reorganize the model constraints and add
Karimi [3] extended and improved this MILP model by using unit- others that, although redundant in the original space, strengthen the
specific time grids and including most of the operational constraints relaxation. This procedure is known as the reformulation lineariza-
found in practice. Both models optimized blend recipes using blend tion technique [25]. The disadvantage is that there is no systematic
indices. Based on these two previous works, Li et al. [4] presented a way of knowing where to act in order to move toward a stronger
unit-specific continuous-time MINLP formulation where nonlinear relaxation.
terms arise from enforcing constant blending rates. A widely used method that guarantees convergence to the global
Castillo and Mahalec [14] developed a three-level decompo- optimal solution is known as spatial branch-and-bound [26,27]. It
sition algorithm through which recipe optimization can be done is an essential element of commercial solvers such as BARON [22],
using linear and/or nonlinear blending rules. They considered the ANTIGONE [23], Couenne [28], and SCIP [29]. Spatial branch-and-
distribution problem, blend-size threshold constraints, parallel bound works by iteratively reducing the domain of the variables,
non-identical blenders, swing tanks, and product-dependent setup one by one, which in turn improves the quality of the convex relaxa-
times. A discrete-time model was formulated for each level. The first tion. Note that if the initial relaxation is weak, due to the presence of
level optimized the blend recipes, the second level approximated many non-convex terms, convergence can be rather slow. It is thus
the production schedule, and the third level computed a detailed important to have good branching strategies and bound-tightening
blend-and-delivery schedule. Due to the large size of the schedul- techniques. Optimality-based bound tightening (OBBT) is an exam-
ing model at the third level for the entire horizon, it was solved in ple of the latter. Although typically applied only at the root node or
subintervals. Solutions computed by this approach were better and up to a limited depth [28], recent results have shown that applying
the execution times for large problems were two orders of mag- OBBT in every node may lead to considerably lower optimality gaps
nitude shorter than those from previous methods [3,13]. In their [30]. OBBT involves solving one minimization and one maximiza-
subsequent work, Castillo and Mahalec [15] introduced a signifi- tion problem for each variable (appearing in non-convex terms) in
cantly modified version of the continuous-time scheduling model order to generate tighter lower and upper bounds. It can be solved
from Li and Karimi [3] (with a smaller number of binary variables sequentially—which has the advantage of generating tighter variable
[16]) for dealing with the third level. Case studies with nonlinear bounds and the disadvantage of being computationally expensive
blend-scheduling problems were solved very close to global opti- when dealing with a large number of variables—or in parallel [31].
mality with short execution times. Bilinear terms are a common source of non-convexities in process
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 191

systems engineering. They can be relaxed using the McCormick shows the results obtained with the proposed algorithm and pro-
envelopes, considering either the full domain of the variables [32] vides a comparison with other methods. The paper ends with con-
or a reduced domain following partitioning [33,34]. Simultaneous clusions in Section 9.
domain partitioning involves adding a new set of binary variables
to the problem and guarantees global optimality in the limit of an 2. Problem statement
infinite number of partitions. Spatial branch-and-bound can thus
be avoided. One critical aspect concerns the scaling of problem size Given a blending system (i.e., storage tanks, blenders, and their
with the number of partitions. Earlier piecewise relaxation tech- interconnections; see Fig. 3), a scheduling horizon, a set of blend
niques [33] scale linearly, while recent formulations scale logarith- components and their corresponding supply and quality profiles
mically. Examples of the latter are described next. along the horizon, a set of products and their minimum and max-
Misener et al. [35] developed a global optimization algorithm imum quality property specifications, a set of delivery orders for
for the standard pooling problem and concluded that the logarith- each product, and the initial inventory levels, it is required to deter-
mic scheme is more advantageous for more than eight partitions. mine the blend recipes, the production and delivery sequences, and
Kolodziej et al. [19] proposed a MINLP formulation for the multi- the inventory profiles of all tanks, while minimizing the cost of the
period pooling problem, in which nonlinearities arise from using the blended materials plus the switching costs (i.e., number of blend
dynamic inventories of tanks as blend components. They employed runs, number of tanks delivering the same order, and product tran-
a radix-based discretization technique that partitions one variable sitions in the swing tanks) and the demurrage costs (i.e., late deliv-
in every bilinear term to obtain a MILP relaxation. This discretiza- eries).
tion technique is known as multiparametric disaggregation [36]. The following constraints are considered:
Castro [37] developed the normalized multiparametric disaggrega- (1) If a blender is to produce a product, it must blend at least a
tion technique (NMDT) [36], which works by discretizing the range minimum amount.
between a variable’s lower and upper bounds (belonging to [0, 1]). (2) A blender can produce at most one product at any time. Once
The advantage is that the number of partitions becomes the same it begins blending, it must operate for some minimum time
for every discretized variable, even if their domain is different (when before it can switch to another product.
using a global discretization level parameter). NMDT has been suc- (3) A blender requires a minimum setup time during a product
cessfully used to solve multi-period blending problems to global changeover.
optimality, both as a stand-alone approach [18,38,39] or integrated (4) A blender can feed at most one product tank at any time (in-
in a spatial branch-and-bound algorithm [30]. dustrial practice).
Overall, contributions from cited works have enabled optimal (5) Product tanks can only store one product at any time.
solutions of gasoline blend-scheduling problems up to a certain (6) Product tanks cannot receive and deliver material at the same
model size. However, when dealing with large-scale problems, the time.
computation and validation of global optimal solutions remain a The assumptions made in this work are:
difficult challenge. (1) The flowrate profile of each component from the upstream
process is piecewise constant.
1.3. Contributions of this work (2) The component quality profile is piecewise constant.
(3) Perfect mixing occurs in each blender.
This work presents a novel deterministic global optimization (4) There is only one tank for a given blend component.
algorithm to solve non-convex MINLP or NLP models with nonlin- (5) Only swing tanks can change their product service (i.e., change
earities that are strictly due to bilinear and/or quadratic terms. The from storing one product to storing another).
main features of this algorithm are: (6) Changeover times between products are negligible for swing
• The use of different linear and piecewise linear relaxation tech- tanks.
niques to derive convex relaxations of the original non-convex (7) For each blender, changeover times between product blend-
model; ing are product dependent but sequence independent.
• The collection of different feasible solutions from the convex (8) Each order involves only one product (any original order in-
relaxation, which provide starting points for a local nonlinear volving different products can be disaggregated into orders of
solver to find feasible solutions of the original model; a single product).
• A dynamic increase in the number of partitions for piecewise (9) All orders are fulfilled during the scheduling horizon.
linear relaxations; In summary, this problem considers the scheduling of blending
• The reduction of the domain of the variables involved in non- and delivery operations, recipe determination, and product alloca-
linear terms by means of an OBBT method; and tion of swing tanks along the scheduling horizon.
• The parallelization of the steps regarding computation of feasi-
ble solutions and the OBBT method. 3. Gasoline blend-scheduling model
The algorithm is tested on different gasoline blend-scheduling
examples. For this class of problems, the results show that the pro- The scheduling model employed in this work is the one pre-
posed algorithm is on par with or better than two leading commer- sented by Castillo and Mahalec [16]. It employs a continuous-time
cial global solvers. formulation, considers nonlinear blending equations, and does not
The rest of this paper is organized as follows: Section 2 describes allow simultaneous receipt and delivery by product tanks. This is
the problem statement and the assumptions made. Section 3 re- a non-convex MINLP model, and it will be denoted as model P (or
views the scheduling model employed in this work and presents the problem P). The main features of the scheduling model are high-
nonlinear equations used for octane blending. Section 4 briefly ex- lighted in this section.
plains the piecewise linear relaxations employed to compute the es- The scheduling model uses unit-specific time slots of varying
timates of the global solution. Section 5 describes the OBBT method length to determine when a specific task needs to be executed in
to reduce the domain of the variables involved in nonlinear terms. each unit (blenders and tanks in this case). We assign a sufficient-
Section 6 presents the steps of the global optimization algorithm. ly high number of time slots, which will likely be higher than the
Section 7 contains the data describing the test examples. Section 8 number of slots required for blending each grade. This ensures that
192 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201

there are sufficient degrees of freedom (enough available switches) variable delivery rates, blender- and product-specific setup times for
to meet the varying product-delivery schedule. the blenders (i.e., idle times for, e.g., cleaning or sensor recalibration
The start time of a unit slot is equal to the end time of the previ- purposes), maximum delivery rates from blend component tanks
ous one. The first unit slot starts at the beginning of the scheduling to the blenders, minimum blend size, and minimum running times
horizon, and the end of the last unit slot matches exactly the end for each blender and product. Other constraints include the material
of the horizon. Blending tasks begin at the start of a time slot, but balances, product composition specifications, product quality speci-
can finish before its end. Delivery tasks from product tanks can start fications, and linear and/or nonlinear blending equations.
and finish within the corresponding slot. It is assumed that com- The difficulty in solving this scheduling model to global optimal-
ponent tanks are continuously receiving material at some specified ity arises from the following factors:
rate (i.e., the supply profile). Time periods are used to delineate the • The significant number of discrete decisions that can be made,
points where changes occur in the supply rates and/or quality of which are directly related to the number of time slots, gasoline
blend components. Time slots are assigned to these time periods. grades, blenders, product tanks, and demand orders (the com-
A time slot must end within its assigned period. However, for com- binatorial nature of the problem);
ponent tanks, the last time slot of a period must end exactly at that • The inclusion of nonlinear blending equations (the non-convex
period’s boundary (in order to properly respect the changes in supply nature of the problem); and
rates and/or quality of blend components). Fig. 4 shows a graphical • All the considered operational constraints.
representation of these unit-specific time slots for a blending system Castillo and Mahalec [16] found that introducing constraints reg­
with two blend component tanks (CT1, CT2), one blender (B1), and ard­ing minimum blend cost and minimum switching cost can im-
two product tanks (PT1, PT2). Unit slots 1 and 2 are pre-assigned to prove the quality of the solution and reduce the execution times for
period 1, while slots 3 and 4 are pre-assigned to period 2. Note that small- to medium-size problems. The minimum blend cost is com-
the optimization has determined that slot 3 in the CT2 grid, and slot 4 puted using the approach delineated in Castillo et al. [40].
in the PT1 grid, have zero length. The nonlinear blending equations are presented next, since they
The objective of the scheduling model is to minimize the blend were rewritten in such a way that nonlinear terms are only bilinear
cost (i.e., materials cost), the switching cost associated with each or quadratic.
blend run, product changeovers in the swing tanks, the number of
“delivery runs” (i.e., the number of time slots used to deliver a spe- Nonlinear blending equations
cific order from a given tank), and the demurrage cost. Delivery runs Eq. (1) to Eq. (19) are the proposed reformulation of the ethyl
are penalized in order to avoid computing delivery schedules that RT-70 model for octane blending [5,41]. Bilinear terms appear in
deliver the same order from several tanks at the same time, and to Eqs. (1), (13), (14), and (18). Quadratic terms appear in Eqs. (15),
minimize intermittent deliveries of the same order from the same (16), (17), and (19). Main sets, subscripts, variables, and parameters
tank. are described next. Set I = {i} consists of the blend components,
Binary variables are employed in the model to determine, at each BL = {bl} is constituted by the blenders, N1 = {n} is the time slots,
time slot, the following discrete decisions: and set QN = {(θ, n)} represents the time slots associated with each
• Which product tank each blender is feeding (one variable for quality profile θ. Variable Vcomp(i, bl, n) indicates the volume of blend
each blender-tank connection); component i to blender bl during slot n. Variable Vblend(bl, n) is the
• What gasoline grade is stored in each product tank (one varia- volume being processed by blender bl during slot n. The volume
ble for each grade-tank pair); and fraction of component i going into blender bl during slot n is denot-
• What demand order each product tank is partially or complete- ed by variable r(i, bl, n). Parameter Qbc(i, e, θ) represents the value
ly fulfilling (one variable for each tank-order connection). of quality property e for blend component i and quality profile θ.
With these binary variables, other discrete decisions can be mod- sens(i, θ) is a parameter known as the octane number sensitivity; it
eled with 0–1 continuous variables, such as: is the difference between the octane numbers, that is, RON – MON,
• What gasoline grade each blender is producing; of blend component i and quality profile θ. The values for the ethyl
• The status of a blender (running or idle); RT-70 model coefficients are taken from Singh et al. [5] and are as
• The transition of a blender from running to being idle, or vice follows: a1 = 0.03224, a2 = 0.00101, a3 = 0, a4 = 0.0445, a5 = 0.00081,
versa; and a6 = −0.0645 × 10−4.
• When a new blend run starts;
• Product transitions in the blenders; and Vcomp ( i, bl , n ) = r ( i, bl , n )Vblend ( bl , n ) ∀i, bl , n ∈ N1 (1)
• Product changeovers in the swing tanks.
The scheduling model also considers variable blending rates,
RON
ravg ( bl , n ) = ∑ r ( i, bl , n ) Qbc ( i, e, θ )
i
 (2)
∀e = “RON”, bl , n ∈ N1, θ : (θ , n ) ∈ QN

MON
ravg ( bl , n ) = ∑ i r ( i, bl , n ) Qbc ( i, e, θ ) (3)
∀e = “MON”, bl , n ∈ N1, θ : (θ , n ) ∈ QN

Olavg ( bl , n ) = ∑ i r ( i, bl , n ) Qbc ( i, e, θ )
(4)
∀e = “Ol”, bl , n ∈ N1, θ : (θ , n ) ∈ QN

Aravg ( bl , n ) = ∑ i r ( i, bl , n ) Qbc ( i, e, θ )
(5)
∀e = “Ar”, bl , n ∈ N1, θ : (θ , n ) ∈ QN 

( bl , n ) = ∑ i r ( i,bl , n ) Qbc ( i, e, θ ) 
sq 2
Olavg (6)

∀e = “Ol”, bl , n ∈ N1, θ : (θ , n ) ∈ QN

( bl , n ) = ∑ i r ( i, bl , n ) Qbc ( i, e, θ ) 
sq 2
Aravg  (7)
Fig. 4. Representation of unit-specific time slots employed in the scheduling model. ∀e = “Ar”, bl , n ∈ N1, θ : (θ , n ) ∈ QN
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 193

sensavg ( bl , n ) = ∑ i r ( i, bl , n ) sens ( i, θ ) ∀i, bl , n ∈ N1  (8) bilinear term (denoted as the discretized variable). The number of
partitions is denoted as NP, and it is assumed that all discretized

RON
sensavg ( bl , n ) = ∑ i r ( i, bl , n ) Qbc ( i, e, θ ) sens ( i, θ )  (9) variables have the same number of partitions. PMCR has a linear
∀e = “RON”, bl , n ∈ N1, θ : (θ , n ) ∈ QN relation between NP and the number of extra binary variables re-
quired per discretized variable, while NMDT exhibits a logarithmic

MON
sensavg ( bl , n ) = ∑ i r ( i, bl , n ) Qbc ( i, e, θ ) sens ( i, θ ) (10) relation. For a more detailed explanation of these methods, the
∀e = “MON”, bl , n ∈ N1, θ : (θ , n ) ∈ QN reader is encouraged to review Refs. [16,37].
The resulting MILP model is denoted as model PR and is a relax-
Qpr ( bl , e, n ) = ravg
RON
( bl , n ) + a1  sensavg
RON
( bl , n ) − rsavg
RON
( bl , n ) ation of problem P. This means that the optimal solution of model
 +a2 Olavg
sq
( bl , n ) − Ol 2avg ( bl , n )  PR is a valid estimate of the global solution of P (in the minimiza-
(11) tion case, this will be a lower bound, LB). Moreover, an estimate
avg ( bl , n ) − 2 ⋅ Ar 3avg ( bl , n ) + Ar 4avg ( bl , n ) 
+ a3  Ar 2sq  of the best possible solution of model PR is a valid estimate of the
∀e = “RON”, bl , n ∈ N1 global optimum of P. Therefore, even if model PR is not solved to
optimality by a MILP solver within a given allocated time, a new
Qpr ( bl , e, n ) = ravg
MON
( bl , n ) + a4  sensavg
MON
( bl , n ) − rsavg
MON
( bl , n ) estimate of the global solution can still be found. The larger the
number of partitions, the closer model PR is to model P; see Fig. 5
+a5 Olavg
(12)
sq
( bl , n ) − Ol 2avg ( bl , n ) for an illustration with an example involving a single discretized
variable.
avg ( bl , n ) − 2 ⋅ Ar 3avg ( bl , n ) + Ar 4avg ( bl , n ) 
+ a6  Ar 2sq 
If the relaxation is tight, then its optimal solution will be very
∀e = “MON”, bl , n ∈ N1 close to the original optimum. Hence, a strategy to find a feasible
solution to the original problem P (in the minimization case, this
RON
rsavg ( bl , n ) = ravgRON ( bl , n ) sensavg ( bl , n ) ∀bl , n ∈ N1  (13)
will be an upper bound, UB) is to initialize P with the optimal solu-
tion of model PR. Since some MILP solvers, such as CPLEX, can store
MON
rsavg ( bl , n ) = ravgMON ( bl , n ) sensavg ( bl , n ) ∀bl , n ∈ N1  (14)
multiple feasible solutions to the MILP problem, potentially leading
Ar 2avg ( bl , n ) = Aravg ( bl , n )
2
∀bl , n ∈ N1  to different solutions of P due to the different starting points, we use
(15)
a multi-start strategy in parallel fashion. Note that, for practical rea-
Ol 2avg ( bl , n ) = Olavg ( bl , n )
2
∀bl , n ∈ N1  (16) sons related to the speed and robustness of commercial solvers, it is
more convenient to solve NLP models instead of MINLPs. This is the
avg ( bl , n ) = Aravg ( bl , n )
Ar 2sq sq 2
∀bl , n ∈ N1  (17) reason why the values of the binary variables are fixed, converting
problem P (MINLP) into PF (NLP). The compact notations of models P,
Ar 3avg ( bl , n ) = Aravg
sq
( bl , n ) Ar 2avg ( bl , n ) ∀bl , n ∈ N1 (18) PR, and PF are as follows:
Model P:
Ar 4avg ( bl , n ) = Ar 2avg ( bl , n )
2
∀bl , n ∈ N1 (19) min f 0 ( x, y )
s.t. f m ( x, y ) ≤ 0 ∀m ∈ M {0}
4. Piecewise linear relaxations
f m ( x, y ) = ∑ ( i , j )∈BLT aijm xi x j + Bm x + Cm y + d m ∀m ∈ M
As mentioned in Section 1, the use of piecewise linear relaxations 0 ≤ xL ≤ x ≤ xU
is becoming more widespread due to the maturity of MILP solvers. x ∈  lx , y ∈ {0, 1}
ly

Piecewise McCormick relaxation (PMCR) and the NMDT will be

employed in this work. These techniques replace each bilinear term Model PR:
in model P with a single variable, thus linearizing the correspond-
ing equations. This single variable is then subject to various linear min f 0R ( x, y )
constraints, which add extra continuous and binary variables to the s.t. f mR ( x, y ) ≤ 0 ∀m ∈ M {0}
model. If equal to 1, these extra binary variables activate a specific f mR ( x, y ) = ∑ ( i , j )∈BLT aijm wij + Bm x + Cm y + d m ∀m ∈ M
interval of the domain (i.e., partition) of one of the variables in the
g nR ( x, w, v, z ) ≤ 0 ∀n ∈ N
g nR ( x, w, v, z ) = H n′ x + An′ w + Bn′ v + Cn′ z + d n′ ∀n ∈ N
0 ≤ xL ≤ x ≤ xU
x ∈ lx , y ∈ {0, 1} , w ∈ lw , v ∈ lv , z ∈ {0, 1}
ly lz

Fig. 5. Accuracy of the relaxation (f0R) with respect to exact representation (f0) of the boundaries of a feasible region increases with the number of partitions (maximization prob-
lem). (a) 10 partitions; (b) 100 partitions.
 min f 0R ( x, y )
s.t. f mR ( x, y ) ≤ 0 ∀m ∈ M {0}
194 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201
f mR ( x, y ) = ∑ ( i , j )∈BLT aijm wij + Bm x + Cm y + d m ∀m ∈ M
g nR ( x, w, v, z ) ≤ 0 ∀n ∈ N objective function is to maximize this variable. In order to compute
g nR ( x, w, v, z ) = H n′ x + An′ w + Bn′ v + Cn′ z + d n′ ∀n ∈ N new bounds, the extra constraint added imposes the condition that
0 ≤ xL ≤ x ≤ xU the value of the relaxed version of the original objective function,
R
that is, f0 (x, y), must be at least as good as the current best feasible
x ∈ lx , y ∈ {0, 1} , w ∈ lw , v ∈ lv , z ∈ {0, 1}
ly lz

Model PF: solution.

min f 0F ( x ) Note that models PR and PRB can use different relaxations. In
this work, model PRB employs standard McCormick envelopes [32]
s.t. f mF ( x ) ≤ 0 ∀m ∈ M {0}
and integrality requirements on variables y are dropped, thus reduc-
f mF ( x ) = ∑ ( i , j )∈BLT aijm xi x j + Bm x + Cm yˆ + d m ∀m ∈ M ing PRB to linear programming (LP). The lower and upper bounds
0 ≤ xL ≤ x ≤ xU of variable xh are updated with the optimal solutions of the corre-
x ∈  lx sponding LP model. Compact notation of model PRB is shown below
for a minimization problem.
Note that, in this section, set M = {m} represents all the origi- Model PRB:
nal constraints, set N = {n} represents all the constraints required
by the piecewise linear relaxation technique, and set BLT = {(i, j)}
xhL = min xh (x U
h = max xh )

represents all the bilinear terms. Variables x and y are the original s.t. f 0
R
( x, y ) ≤ UB
continuous and binary variables, respectively, and v and z are the f m ( x, y ) ≤ 0 ∀m ∈ M {0}
R

extra continuous and binary variables, respectively, required by f mR ( x, y ) = ∑ ( i , j )∈BLT aijm wij + Bm x + Cm y + d m ∀m ∈ M
the relaxation strategy. Variable wij is the continuous variable that
replaces the bilinear term xi xj. Scalars lx, ly, lw, lv, and lz represent g kRB ( x, w ) ≤ 0 ∀k ∈ K
the size of vectors x, y, w, v, and z, respectively. Parameters xL and xU g kRB ( x, w ) = H k′ x + Ak′ w + d k′ ∀k ∈ K
are respectively the lower and upper bounds of the x variables. Note 0 ≤ xL ≤ x ≤ xU
that quadratic terms can be treated as bilinear terms. x ∈  lx , y ∈ [ 0, 1] , w ∈  lw
ly

5. Tightening bounds on the variables

OBBT consists of solving these LP models for all the variables
Model PR becomes tighter (i.e., closer to model P) as the number involved in nonlinear terms, in a parallel framework, to reduce exe-
of partitions of the discretized variables is increased. However, an cution times. Thus, the bounds will generally be weaker than when
increase in the number of partitions produces an increment in the solving the problems sequentially. Since the number of instances
size of model PR and, after a certain number of partitions, model PR to solve can be very large, instances are solved in different blocks.
can become computationally intractable. Therefore, another tech- These blocks are defined by a maximum number of problems to
nique is required in order to avoid the necessity of a large number be solved in parallel. After one block is solved, the corresponding
of partitions. In this work, an OBBT method is employed [34,42]. The bounds are updated and then the next block is solved. Fig. 7 shows
idea is to reduce the domain of the variables involved in nonlinear the flowchart of the OBBT method. Note that OBBT is applied only
terms by computing new bounds of these variables by solving two once per variable.
optimization problems (a maximization problem and a minimiza-
tion problem per variable). This is done after a new and better feasi- 6. Global optimization algorithm
ble solution to P is computed. After reducing the domain of the vari-
ables, model PR becomes closer to P without increasing the number The steps of the proposed global optimization algorithm are
of partitions, as shown in Fig. 6. presented next for a minimization problem. Fig. 8 shows the corre-
The mathematical model used in OBBT is denoted as model PRB, sponding flowchart. Note that the algorithm can be applied to any
which is constructed as a relaxation of P, but with a different ob- MINLP problem with nonlinearities exclusive to bilinear or quadratic
jective function and an extra constraint. To compute a lower bound terms.
of variable xh, that is, xhL, the objective function is to minimize this (1) Initialize algorithm parameters. Define the number of parti-
U
variable. To compute an upper bound of variable xh, that is, xh , the tions to be used {NP1, NP2, …, NPlast} and set NP = NP1. Set the

Fig. 6. Accuracy of the relaxation increases when the domain of the variables involved in nonlinear terms is reduced. (a) 10 partitions with x ∈ [0, 4.5]; (b) 10 partitions with
x ∈ [2.25, 2.7].
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 195

lower bound LB = −∞, upper bound UB = +∞, total number of

iterations counter ITtotal = 1, iterations with the same number
of partitions ITsameNP = 1, maximum number of total iterations
max
ITtotal , maximum number of iterations with the same number
max
of partitions ITsameNP max
, maximum total time TIMEtotal , and mini-
mum relative tolerance ε.
(2) Lower bound computation. Solve MILP model PR using the
CPLEX solver with parallel and solution pool options active.
Update LB with the best possible solution from CPLEX, if this
value is greater than the previous LB.
(3) Upper bound computation. Use the solutions stored in the
CPLEX solution pool as starting points for NLP model PF.
Solve NLP model PF instances in parallel using a local non-
linear solver. Update UB if any of the computed solutions is
feasible and has a smaller objective function value than the
previous UB.
(4) Update optimality gap. The following formula is used in this
step: OptGap = (UB – LB)/LB × 100.
max
(5) Check termination criteria. Stop if OptGap ≤ ε, if ITtotal = ITtotal , if
max
the total execution time is equal to or greater than TIMEtotal , or
if the number of partitions has already reached NPlast. Other-
wise, continue to Step 6.
(6) If the upper bound UB did not improve in Step 3, or if ITsameNP
max
= ITsameNP , continue to Step 7. Otherwise, reduce the domain
of the variables involved in nonlinear terms using the OBBT
method described in Section 5; set ITtotal = ITtotal + 1 and ITsameNP
= ITsameNP + 1, and then go back to Step 2.
(7) Increase the number of partitions to the next specified value.
Set ITtotal = ITtotal + 1 and go back to Step 2.
Although the main elements of the algorithm have already been
proposed (e.g., PMCR, NMDT, OBBT), the novelty is related to the way in
Fig. 7. Flowchart of the OBBT method. which they are implemented. More specifically: ① the CPLEX solution

Fig. 8. Flowchart of the global optimization algorithm.

196 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201

pool is used to store starting points for model PF, ② instances of convert the blend components’ RBN values to RON values before
model PF are solved in parallel, ③ OBBT is applied to blocks of var- the optimization runs (i.e., they do not appear in the optimization
iables and in a parallel framework, and ④ no branching strategy problems).
is employed.
RBN = RON + 11.5 0 ≤ RON ≤ 85 (20)
7. Case studies RBN = exp ( 0.0135RON + 3.422042 ) RON > 85(21)

The tests in this paper consist of Examples 4, 8, 12, and 14 from Table 1 describes the size of the blending system examples. In-
Li and Karimi [3]. The difference in this work is that the ethyl RT- formation about the periods, their duration, their corresponding
70 models are considered for blending RON and MON properties (as time slots, and the orders that can be delivered within such periods
described in Section 3.1) instead of blend indices. RON index corre- is presented in Table 2. RON and MON values and their respective
lations from Li et al. [13], shown in Eq. (20) and Eq. (21), were used specifications are shown in Table 3. Table 4 presents the statistics
to compute the actual RON values from the blend indices given by regarding the size of model P when not using the constraints on the
Li and Karimi [3]. RBN denotes the research octane number blend minimum blend and switching costs. When using such constraints,
index. MON values were assumed in this work and the correspond- four equations are added to the model (minimum blend cost, mini-
ing minimum product specifications were set equal to zero; there- mum number of delivery runs, minimum number of blend runs, and
fore, MON constraints will not be active at the optimal solution. minimum number of product changeovers in the swing tanks). Note
Quality properties of blend components are assumed to be known that the size of the blending system and its corresponding schedul-
(recall Section 2); therefore, Eq. (20) and Eq. (21) are only used to ing model increase from Example 4 to Example 14.

Table 1
Summary of the blending system examples.
Example ID Number of blenders Number of orders Number of products Number of product tanks Number of quality properties
4 1 15 4 11 9
8 2 20 4 11 9
12 2 35 5 11 9
14 3 45 5 11 9

Table 2
Periods, duration, time slots, and orders that can be delivered in each period.
Example ID Period Duration (h) Slots Orders that can be delivered
4 1 100 1, 2 O1–O7, O12–O15
2 92 3, 4 O8–O11
8 1 80 1, 2 O1–O7, O12–O19
2 70 3, 4 O8
3 42 5, 6 O8–O11, O20
12 1 50 1–3 O1–O7, O12, O13, O15, O19, O33
2 50 4–6 O14–O18, O27, O28, O33
3 50 7–9 O8, O21, O24, O29–O32, O34, O35
4 42 10–12 O8–O11, O20, O22, O23, O25, O26
14 1 50 1–3 O1–O7, O12, O13, O15, O19, O26
2 50 4–7 O14–O18, O26
3 50 8–10 O8, O21, O24, O27–O31, O45
4 42 11–13 O8–O11, O20, O22, O23, O25, O32–O44

Table 3
RON and MON values and specifications.
Property Blend components Product specifications [min, max]
C1 C2 C3 C4 C5 C6 C7 C8 C9 P1 P2 P3 P4 P5
RON 75 90.3 95.6 97.3 83 100 115 118 81 [95, 200] [96, 200] [94, 200] [90, 200] [98, 200]
MON 66 80.8 80.5 91.7 74 100 109 100 72 [0, 200] [0, 200] [0,200] [0, 200] [0, 200]

Table 4
Statistics of model P.
Example ID Number of equations Number of variables Number of binary variables Number of bilinear terms
4 6 207 2 503 433 168
8 9 297 3 323 553 336
12 23 087 8 170 1 317 672
14 32 574 10 693 1 628 1 092
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 197

8. Results BARON does not find a feasible solution for Examples 12 and 14
within 1 h. Based on the log files generated by commercial solvers,
All the examples were solved on a computing machine Intel ® it seems that BARON relies more on the branching procedure, while
Core™ i7-4710HQ central processing unit (CPU), 2.50 GHz, 12 GB ANTIGONE focuses more on the MILP relaxation and bound-tight-
random-access memory (RAM), Windows 10 (8-core). The global ening steps (as our proposed algorithm does). Feasible integer solu-
optimization algorithm was implemented in Python 2.7. The Py- tions for scheduling problems may be found only at deep nodes in
thon script generates general algebraic modeling system (GAMS) the branch-and-bound tree [43], which can be of significant size
files with the corresponding mathematical models, which are then when the number of binary variables is large.
solved by employing the GAMS-Python application program inter- In Examples 4 and 8, the algorithm and the commercial solvers
face (API). The selected solvers were CPLEX 12.6 for model PR and compute similar optimality gaps. For Examples 12 and 14, BARON
model PRB, and CONOPT 3 for model PF. cannot compute an optimality gap (no feasible solution was found),
The global optimization algorithm termination criteria were as and GO-PMCR obtains a lower optimality gap than GO-NMDT and
follows: 0.01% optimality gap or 3600 s (1 h). There was no limit on ANTIGONE.
the total number of iterations, nor on the iterations with the same Both commercial solvers and the proposed algorithm did not
number of partitions. The numbers of partitions in model PR when solve all four examples to the desired tolerance within 1 h. The
using PMCR were {2, 4, 8, 16, 32}, and when using NMDT were {10, times reported in Table 5 are the times in which the best solution
100, 1000}. was found. GO-PMCR and GO-NMDT require shorter times than both
The termination criteria for the MILP problems (instances of commercial solvers. GO-NMDT is significantly faster than GO-PMCR
model PR) were: an optimality gap of 0.01% or 600 s. The CPLEX only in the small-sized Example 4. Note that, for the number of
parallel option was active (in deterministic mode) with a maximum partitions selected, the size of the MILP relaxation grows faster with
number of threads equal to 8. In addition, the CPLEX solution pool GO-NMDT than with GO-PMCR. Therefore, MILP relaxations are of-
option was active with a maximum pool capacity of 30 and the ten faster to solve to optimality with GO-PMCR. However, GO-PMCR
replacement option that generates diverse solutions. Thus, a maxi- will require more iterations. Based on the three factors considered
mum of 30 instances of model PF was solved per iteration using the (i.e., best solution found, optimality gap, and time to best solution),
GAMS parallel computing grid facility. CONOPT 3 default termina- GO-PMCR shows the best performance. Fig. 9 shows the total num-
tion criteria were used. ber of binary variables in the relaxation of the scheduling model
The termination criteria for the LP problems (instances of model (i.e., model PR) when using PMCR and NMDT, at each iteration of
PRB) were: optimality or 60 s. A maximum number of 100 LP prob- the algorithm and for each example. It shows that, for the selected
lems were solved in parallel using the GAMS parallel computing grid partition values, PMCR requires fewer binary variables in the first
facility. 4–5 iterations than NMDT at any iteration. This is expected since the
For comparison purposes, the global commercial solvers BARON partitions when using PMCR are fewer than 8 in those iterations,
15.9 [33] and ANTIGONE 1.1 [34] were employed to solve the orig- and NMDT starts with 10. Note that flat sections in the curves in-
inal model P using the same termination criteria as the proposed dicate that the OBBT method was applied instead of increasing the
global optimization algorithm. number of partitions.
Section 8.1 presents the results obtained when not including the The major differences between the proposed algorithm and BAR-
constraints on the minimum blend cost and minimum switching ON are:
cost, while Section 8.2 shows the results when such constraints are • The use of piecewise relaxation methods for bilinear terms in-
added to the model. A comparison with previously published heu- stead of standard McCormick envelopes; and
ristic methods is included in Section 8.3. • Dynamically increasing the number of partitions in order to
tighten the MILP relaxation, instead of implementing a branch-
8.1. Not using constraints on the minimum blend and switching costs ing strategy.
These features seem to be adequate for the scheduling problems
A comparison of the results obtained by the proposed algorithm presented here. We do not claim that our proposed algorithm will al-
with those obtained by commercial solvers is presented in Table 5. ways be better than BARON when solving a different type of problem.
For simplicity, we refer to our proposed algorithm as GO-PMCR Our proposed algorithm and ANTIGONE perform similarly, but
when it uses piecewise McCormick relaxation to construct model differ mainly in the following areas:
PR and as GO-NMDT when it employs the NMDT. • The use of NMDT as a piecewise relaxation technique;
ANTIGONE, BARON, and GO-PMCR compute the same solution for • How the number of partitions is increased in each iteration;
Examples 4 and 8. GO-NMDT computes the same solution for Exam- • When and how to apply OBBT; and
ple 4, but the final solution for Example 8 is slightly higher. GO-PMCR • The use of the CPLEX solution pool to store feasible solutions
computes better solutions than GO-NMDT and ANTIGONE in all ex- from the MILP relaxations and use them as starting points to
amples. In this work, this is because GO-PMCR can use more distinct solve the NLP problem.
numbers of partitions (2, 4, 8, 16, 32) than GO-NMDT (10, 100, 1000); Finally, ANTIGONE and BARON can handle more than just biline-
thus, it generates more feasible solutions from the MILP relaxation. ar and quadratic terms; in addition, they apply other mathematical

Table 5
Summary of results (not using constraints on the minimum blend and switching costs).
Example ID Best solution found (1000 USD) Optimality gap (%) Time to best solution (s)
ANTIGONE BARON GO-PMCR GO-NMDT ANTIGONE BARON GO-PMCR GO-NMDT ANTIGONE BARON GO-PMCR GO-NMDT
4 4 633 4 633 4 633 4 633 6.90 6.37 6.31 6.37 2 350 930 616 120
8 8 203 8 203 8 203 8 223 9.53 9.53 9.19 9.49 1 708 3 273 714 1 391
12 16 650 NF 15 408 15 440 20.51 NA 14.00 14.22 1 631 3 600 1 411 1 438
14 21 360 NF 21 316 31 639 12.50 NA 12.31 40.92 3 600 3 600 2 911 2 904
NF = not found; NA = not available.
198 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201

techniques capable of improving performance (e.g., reformulation- Overall, GO-PMCR shows the best performance once again. For il-
linearization technique, cutting planes, feasibility-based bound lustration purposes, the blend and delivery schedules computed for
tightening, different branching strategies, etc.). Example 14 by the algorithm using PMCR are shown in Fig. 10 and
Fig. 11, respectively.
8.2. Using constraints on the minimum blend and switching costs
8.3. Comparison with heuristic methods
For this case, the results computed by the algorithm and com-
mercial solvers are presented in Table 6. The most notable differenc- In this section, the proposed algorithm is compared with previ-
es with respect to Table 5 are the smaller optimality gaps and the ously published heuristic methods [15,21]. Table 7 [15,21] contains
shorter times for Examples 4 and 8. the best solution found by those methods and the time required to
Both commercial solvers and the algorithm find similar solutions compute such solutions. Note that heuristic methods do not com-
for Examples 4 and 8. BARON still does not find a feasible solution for pute an optimality gap since they aim to find close-to-optimal solu-
Examples 12 and 14 within the allocated time. GO-PMCR computes tions very rapidly, and they do not spend time estimating and refin-
better solutions than GO-NMDT for Examples 12 and 14, which in ing the value of the global optimal solution. These heuristic methods
turn has better solutions than ANTIGONE. Note that the addition of are tailored to the examples used in this work. These methods con-
bounds caused an increment in the number of solutions for Example 8. struct the final solution by decomposing the original problem into
Since the solutions from Section 8.1 are still feasible even with the different levels, each one with different accuracy and complexity.
inclusion of the constraints regarding minimum blend and switching Short execution times are achieved by solving the least complex
costs, this suggests that such constraints are affecting the solvers. This level first and then, in each subsequent level, fixing the values of the
effect is also observed in ANTIGONE in Examples 12 and 14. most important variables to those from the previous level’s solution.
Regarding optimality gaps, most of the same observations as in The objective function of the scheduling model employed in this
Section 8.1 can be made. Similar optimality gaps are computed by all work is the same as the one used by Castillo and Mahalec [15]. This
methods for Examples 4 and 8. For Examples 12 and 14, GO-PMCR objective function penalizes each individual blend run, even when
calculates lower optimality gaps than GO-NMDT and ANTIGONE. the same product is being blended in contiguous blend runs. On the
Both commercial solvers and the proposed algorithm solve Ex- other hand, Cerdá et al. [21] did not penalize the number of indi-
ample 4 to the desired tolerance; BARON is the slowest. For Example vidual blend runs, but only penalized the product transitions in the
12, the time to the best solution required by GO-PMCR is larger than blenders. We show the adjusted values of the solutions reported by
that required by GO-NMDT; however, it must be considered that Cerdá et al. [21]; that is, individual blend runs are penalized.
GO-PMCR computes a better solution. All methods find the same solution for Example 4. In general, all
the methods compute very similar solutions for the remaining ex-
amples. Solutions from Cerdá et al. [21] have higher costs for Exam-
ples 8, 12, and 14 because they did not originally penalize individual
blend runs. The method from Cerdá et al. [21] might compute sim-
ilar solutions to those from the other methods if it used the same
objective function.
Heuristic methods still require smaller execution times than the
proposed global optimization algorithm. This is expected, because
those methods do not involve as many steps as global optimization
techniques. The proposed global optimization algorithm does not
find solutions of the same quality as quickly as the two selected heu-
ristic methods. To compute feasible solutions in each iteration, our
proposed algorithm needs to first solve a MILP (i.e., model PR). The
solution of the MILP is the most time-consuming step, thus reducing
the speed required to compute new feasible solutions. Moreover, the
small number of partitions at the beginning of the algorithm may
result in weak MILP relaxations, which generate starting points for
the NLP models that are far from the global optimum.
These results indicate the need to improve the corresponding step
to compute feasible solutions, or to simply integrate heuristic meth-
ods into the proposed algorithm.

9. Conclusions

Fig. 9. Number of binary variables in model PR at each iteration of the algorithm. In this work, we presented a global optimization algorithm that

Table 6
Summary of results (using constraints on the minimum blend and switching costs).
Example ID Best solution found (1000 USD) Optimality gap (%) Time to best solution (s)
ANTIGONE BARON GO-PMCR GO-NMDT ANTIGONE BARON GO-PMCR GO-NMDT ANTIGONE BARON GO-PMCR GO-NMDT
4 4 633 4 633 4 633 4 633 0.01 0.01 0.01 0.01 26 296 30 14
8 8 207 8 204 8 206 8 204 0.05 0.02 0.04 0.02 557 1 218 103 140
12 23 590 NF 15 384 15 403 34.80 NA 0.01 0.13 3 333 3 600 2 674 742
14 23  520 NF 21 270 21 360 9.68 NA 0.13 0.55 1 636 3 600 2 574 2 845
NF = not found; NA = not available.
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 199

Fig. 10. Blend schedule computed for Example 14 by the proposed algorithm using PMCR. Kbbl is short for kilobarrel, 1 kbbl = 158.9873 m3.

Fig. 11. Delivery schedule computed for Example 14 by the proposed algorithm using PMCR.

Table 7
Comparison with heuristic methods.
Example ID Best solution found (1000 USD) Time to best solution (s)
Castillo and Cerdá et al. [21] Cerdá et al. [21] GO-PMCR GO-NMDT Castillo and Cerdá et al. [21] GO-PMCR GO-NMDT
Mahalec [15] adjusted values Mahalec [15]
4 4 633 4 613 4 633 4 633 4 633 3 0.4 30 14
8 8 203 8 163 8 223 8 206 8 204 6 7.5 103 140
12 15 403 15 342 15 442 15 384 15 403 17 31.0 2 674 742
14 21 263 21 181 21 301 21 270 21 360 24 21.0 2 574 2 845

can solve MINLP problems with bilinear and quadratic terms. The number of partitions is increased during the algorithm.
algorithm computes estimates of the global solution by constructing To avoid a rapid increase in the model size due to a large number
and solving MILP problems that are relaxations of the original prob- of partitions, an OBBT method is used. The MILP relaxation will be
lem obtained by using either PMCR or NMDT. These methods discre- closer to the original problem if the number of partitions stays the
tize the domain of one of the variables of a bilinear term into several same but the domain of the variables is reduced. The OBBT method
partitions, and introduce extra binary and continuous variables into solves several LPs in a parallel setting.
the model. To improve the estimates of the global optimum, the The CPLEX solution pool is active and stores different feasible
200 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201

solutions found during the branch-and-bound procedure to solve Continuous variables

the MILP relaxation. These solutions are then used as starting points Aravg (bl, n) Volumetric average of the aromatics content of the
for a nonlinear solver to find feasible solutions to the original prob- processed material by blender bl during slot n
lem. This step is also parallelized. Arsqavg (bl, n) Volumetric average of the squared value of the aro-
We showed that the proposed algorithm can be used to schedule matics content of the processed material by blender bl
gasoline-blending operations, taking into consideration the distri- during slot n
bution problem and the most important operational constraints. We Ar2avg (bl, n) Squared value of Aravg (bl, n)
employ a continuous-time MINLP scheduling model [16] where the Ar2sqavg (bl, n) Squared value of Arsqavg (bl, n)
ethyl RT-70 models are used for octane blending. Ar3avg (bl, n) Product of Arsqavg (bl, n) and Ar2avg (bl, n)
The proposed algorithm was compared with two commercial Ar4avg (bl, n) Squared value of Ar2avg (bl, n)
solvers and two heuristic methods. The elements under evaluation Olavg (bl, n) Volumetric average of the olefins content of the pro-
were: the best solution found, the corresponding optimality gap, cessed material by blender bl during slot n
and the time to best solution. The proposed algorithm with PMCR Ol sqavg (bl, n) Volumetric average of the squared value of the olefins
showed a better performance than with NMDT. In our large-sized content of the processed material by blender bl during
examples, the proposed algorithm with either PMCR or NMDT per- slot n
formed better than the commercial global solvers. This result shows Ol2avg (bl, n) Squared value of Ol sqavg (bl, n)
that further research on this algorithm may be very promising. Both Qpr (bl, e, n) Value of quality property e of the processed material
selected heuristic methods provided good solutions in shorter ex- by blender bl during slot n
ecution times than the global algorithms. This result indicates that r (i, bl, n) Volume fraction of blend component i going into
the step to compute feasible solutions can still be improved. blender bl during slot n
We tested the performance of the algorithm by solving the sched- r (bl, n) Volumetric average of the motor octane number of the
uling model for two scenarios: ① not including lower bounds on the processed material by blender bl during slot n
blend cost and switching costs, and ② including such bounds. The r (bl, n) Volumetric average of the research octane number of
first problem is harder to solve and is representative of a kind of mod- the processed material by blender bl during slot n
el one may write without diligently trying to reduce the search space rs (bl, n) Product of r (bl, n) and sensavg (bl, n)
as much as possible. Adding a tight lower bound to the blending cost rs (bl, n) Product of r (bl, n) and sensavg (bl, n)
as a constraint, as well as adding the lower bound to the switching sensavg (bl, n) Volumetric average of the octane number sensitivity
costs, enables algorithms to search smaller spaces and improves their of the processed material by blender bl during slot n
performance. This result also indicates that the relaxations are still sens (bl, n) Volumetric average of the octane number sensitivity
not tight enough. Future work will include the derivation and addi- times the motor octane number
tion of redundant and symmetry-breaking constraints; the testing sens (bl, n) Volumetric average of the octane number sensitivity
of a different relaxation scheme for quadratic, cubic, and higher times the research octane number
order terms (e.g., outer approximation); and the modification of the Vblend (bl, n) Volume being processed by blender bl during slot n
bound-tightening method in order to speed up the algorithm. Vcomp (i, bl, n) Volume of blend component i transferred to blender bl
during slot n
Acknowledgements
References
Support by Ontario Research Foundation, McMaster Advanced
Control Consortium, and Fundação para a Ciência e Tecnologia (In- [1] Harjunkoski I, Maravelias CT, Bongers P, Castro PM, Engell S, Grossmann IE, et al.
Scope for industrial applications of production scheduling models and solution
vestigador FCT 2013 program and project UID/MAT/04561/2013) is
methods. Comput Chem Eng 2014;62:161–93.
gratefully acknowledged. [2] Méndez CA, Grossmann IE, Harjunkoski I, Kaboré P. A simultaneous optimization
approach for off-line blending and scheduling of oil-refinery operations. Comput
Compliance with ethics guidelines Chem Eng 2006;30(4):614–34.
[3] Li J, Karimi I. Scheduling gasoline blending operations from recipe determination
to shipping using unit slots. Ind Eng Chem Res 2011;50(15):9156–74.
Pedro A. Castillo Castillo, Pedro M. Castro, and Vladimir Mahalec [4] Li J, Xiao X, Floudas CA. Integrated gasoline blending and order delivery opera-
declare that they have no conflict of interest or financial conflicts to tions: Part I. Short-term scheduling and global optimization for single and multi-
period operations. AIChE J 2016;62(6):2043–70.
disclose. [5] Singh A, Forbes JF, Vermeer PJ, Woo SS. Model-based real-time optimization of
automotive gasoline blending operations. J Process Contr 2000;10(1):43–58.
Nomenclature [6] Joly M, Pinto JM. Mixed-integer programming techniques for the scheduling of
fuel oil and asphalt production. Chem Eng Res Des 2003;81(4):427–47.
[7] Floudas CA, Lin X. Continuous-time versus discrete-time approaches for sched-
Sets and subscripts uling of chemical processes: A review. Comput Chem Eng 2004;28(11):2109–29.
BL = {bl} Blenders [8] Sundaramoorthy A, Maravelias CT. Computational study of network-based
mixed-integer programming approaches for chemical production scheduling.
E = {e} Quality properties Ind Eng Chem Res 2011;50(9):5023–40.
I = {i} Blend components and corresponding storage tanks [9] Maravelias CT. General framework and modeling approach classification for
N1 = {n} Time slots chemical production scheduling. AIChE J 2012;58(6):1812–28.
[10] Jia Z, Ierapetritou M. Mixed-integer linear programming model for gasoline
QN = {(θ, n)} Time slot n is associated with the period with quality
blending and distribution scheduling. Ind Eng Chem Res 2003;42(4):825–35.
profile θ [11] Jia Z, Ierapetritou M. Efficient short-term scheduling of refinery operations based
on a continuous time formulation. Comput Chem Eng 2004;28(6–7):1001–19.
[12] Glismann K, Gruhn G. Short-term scheduling and recipe optimization of blend-
Parameters
ing processes. Comput Chem Eng 2001;25(4–6):627–34.
a1, a2, ..., a6 Coefficients for the ethyl RT-70 model [13] Li J, Karimi I, Srinivasan R. Recipe determination and scheduling of gasoline
Qbc (i, e, θ) Value of quality property e of blend component i dur- blending operations. AIChE J 2010;56(2):441–65.
ing quality profile θ [14] Castillo PAC, Mahalec V. Inventory pinch based, multiscale models for inte-
grated planning and scheduling—Part II: Gasoline blend scheduling. AIChE J
sens (i, θ) Octane number sensitivity, i.e., octane difference RON 2014;60(7):2475–97.
– MON for blend component i during quality profile θ [15] Castillo PAC, Mahalec V. Inventory pinch gasoline blend scheduling algo-
 P.A. Castillo Castillo et al. / Engineering 3 (2017) 188–201 201

rithm combining discrete- and continuous-time models. Comput Chem Eng ametric disaggregation. Optim Methods Softw. Epub 2016 Dec 13.
2016;84:611–26. [31] Castillo PC, Castro PM, Mahalec V. Global optimization algorithm for large-
[16] Castillo PAC, Mahalec V. Improved continuous-time model for gasoline blend scale refinery planning models with bilinear terms. Ind Eng Chem Res
scheduling. Comput Chem Eng 2016;84:627–46. 2017;56(2):530–48.
[17] Lotero I, Trespalacios F, Grossmann IE, Papageorgiou DJ, Cheon MS. An MILP- [32] McCormick GP. Computability of global solutions to factorable noncon-
MINLP decomposition method for the global optimization of a source based vex programs: Part I—Convex underestimating problems. Math Program
model of the multiperiod blending problem. Comput Chem Eng 2016;87:13–35. 1976;10(1):147–75.
[18] Castro PM. New MINLP formulation for the multiperiod pooling problem. AIChE [33] Karuppiah R, Grossmann IE. Global optimization for the synthesis of integrated
J 2015;61(11):3728–38. water systems in chemical processes. Comput Chem Eng 2006;30(4):650–73.
[19] Kolodziej SP, Grossmann IE, Furman KC, Sawaya NW. A discretization-based ap- [34] Castro PM. Tightening piecewise McCormick relaxations for bilinear problems.
proach for the optimization of the multiperiod blend scheduling problem. Com- Comput Chem Eng 2015;72:300–11.
put Chem Eng 2013;53:122–42. [35] Misener R, Thompson JP, Floudas CA. APOGEE: Global optimization of standard,
[20] Cerdá J, Pautasso PC, Cafaro DC. A cost-effective model for the gasoline blend generalized, and extended pooling problems via linear and logarithmic parti-
optimization problem. AIChE J 2016;62(9):3002–19. tioning schemes. Comput Chem Eng 2011;35(5):876–92.
[21] Cerdá J, Pautasso PC, Cafaro DC. Optimizing gasoline recipes and blending opera- [36] Kolodziej S, Castro PM, Grossmann IE. Global optimization of bilinear programs
tions using nonlinear blend models. Ind Eng Chem Res 2016;55(28):7782–800. with a multiparametric disaggregation technique. J Glob Optim 2013;57(4):
[22] Tawarmalani M, Sahinidis NV. A polyhedral branch-and-cut approach to global 1039–63.
optimization. Math Program 2005;103(2):225–49. [37] Castro PM. Normalized multiparametric disaggregation: An efficient relaxation
[23] Misener R, Floudas CA. ANTIGONE: Algorithms for continuous/integer global for mixed-integer bilinear problems. J Glob Optim 2016;64(4):765–84.
optimization of nonlinear equations. J Glob Optim 2014;59(2):503–26. [38] Castro PM, Grossmann IE. Global optimal scheduling of crude oil blending oper-
[24] Boland N, Kalinowski T, Rigtering F. New multi-commodity flow formulations for ations with RTN continuous-time and multiparametric disaggregation. Ind Eng
the pooling problem. J Glob Optim 2016;66(4):669–710. Chem Res 2014;53(39):15127–45.
[25] Sherali HD, Alameddine A. A new reformulation-linearization technique for bi- [39] Castro PM. Source-based discrete and continuous-time formulations for the
linear programming problems. J Glob Optim 1992;2(4):379–410. crude oil pooling problem. Comput Chem Eng 2016;93:382–401.
[26] Ryoo HS, Sahinidis NV. A branch-and-reduce approach for global optimization. J [40] Castillo PAC, Mahalec V, Kelly JD. Inventory pinch algorithm for gasoline blend
Glob Optim 1996;8(2):107–38. planning. AIChE J 2013;59(10):3748–66.
[27] Smith EMB, Pantelides CC. Global optimization of nonconvex MINLPs. Comput [41] Healy WC, Maassen CW, Peterson RT. A new approach to blending octanes. In:
Chem Eng 1997;21(Suppl):S791–6. Proceedings of the 24th Midyear Meeting of American Petroleum Institute’s Di-
[28] Belotti P, Lee J, Liberti L, Margot F, Wächter A. Branching and bounds tightening vision of Refining; 1959 May 27; New York, US; 1959. p. 132–136.
techniques for non-convex MINLP. Optim Methods Softw 2009;24(4–5):597– [42] Castro PM, Grossmann IE. Optimality-based bound contraction with multipar-
634. ametric disaggregation for the global optimization of mixed-integer bilinear
[29] Achterberg T. SCIP: Solving constraint integer programs. Math Program Comput problems. J Glob Optim 2014;59(2):277–306.
2009;1(1):1–41. [43] Kallrath J. Planning and scheduling in the process industry. OR Spectrum
[30] Castro PM. Spatial branch-and-bound algorithm for MIQCPs featuring multipar- 2002;24(1):219–250.