Computers and Chemical Engineering 35 (2011) 2750–2766

Contents lists available at ScienceDirect

Computers and Chemical Engineering

journal homepage: www.elsevier.com/locate/compchemeng

A new Lagrangian decomposition approach applied to the integration of refinery

planning and crude-oil scheduling
Sylvain Mouret a , Ignacio E. Grossmann a,∗ , Pierre Pestiaux b
a
Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, PA, 15213, USA
b
Total Refining & Marketing, Research Division, 76700 Harfleur, France

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

Article history: The aim of this paper is to introduce a methodology to solve a large-scale mixed-integer nonlinear
Received 20 September 2010 program (MINLP) integrating the two main optimization problems appearing in the oil refining indus-
Received in revised form 27 February 2011 try: refinery planning and crude-oil operations scheduling. The proposed approach consists of using
Accepted 25 March 2011
Lagrangian decomposition to efficiently integrate both problems. The main advantage of this technique
Available online 2 April 2011
is to solve each problem separately. A new hybrid dual problem is introduced to update the Lagrange mul-
tipliers. It uses the classical concepts of cutting planes, subgradient, and boxstep. The proposed approach
Keywords:
is compared to a basic sequential approach and to standard MINLP solvers. The results obtained on a
Refinery planning
Crude-oil scheduling
case study and a larger refinery problem show that the new Lagrangian decomposition algorithm is more
Mixed-integer nonlinear programming robust than the other approaches and produces better solutions in reasonable times.
Lagrangian decomposition © 2011 Elsevier Ltd. All rights reserved.

1. Introduction solution. To the best of our knowledge, Xpress-SLP (FICO) is the only
standalone solver implementing an SLP algorithm.
The oil refining industry is a prolific field for the application of A major issue with refinery planning is that most models are
mathematical programming techniques (see Bodington & Baker, single-period models where the refinery is assumed to operate in
1990). Refinery operators have to make decisions on the logis- the same state over the whole planning period (typically 1 month).
tics operations and on crude blend qualities taking into account Therefore, the planning solution is used as a tactical goal for refin-
a large number of crude-oils, finished products such as liquified ery operators rather than as an operational tool. In particular, CDU
petroleum gas, gasoline, diesel fuel, and a wide variety of high (Crude Distillation Unit) feedstock decisions returned by the refin-
flexibility production units involving many different chemical pro- ery planning problem are usually not applicable in the field due
cesses. Furthermore, the economic impact of optimizing operations to crude logistics constraints. These are described in the crude-
can be very significant (see Kelly & Mann, 2003). oil operations scheduling problem, which includes unloading from
The refinery planning problem often involves the pooling prob- crude-oil tankers, preparation of crude-blends, and CDU feed charg-
lem (see Audet et al., 2004; Floudas & Visweswaran, 1993; Foulds ing. This problem has been addressed since the late 90s by many
et al., 1992; Misener & Floudas, 2009; Quesada & Grossmann, 1995), different research groups (see Jia et al., 2003; Lee, Pinto, Grossmann,
which has been addressed since the early 80s and usually con- & Park, 1996; Moro and Pinto, 2004; Mouret et al., 2009).
sists of optimizing feedstocks, unit settings, as well as final product Although, integration of planning and scheduling has recently
blending and shipping. Some examples of nonlinear refinery plan- been addressed in the context of multiproduct continuous and
ning problems including pooling constraints and nonlinear process batch production plants (see Erdirik-Dogan & Grossmann, 2008;
models can be found in Pinto and Moro (2000), Li, Hui, and Li (2005), Maravelias & Sung, 2009), very little work has been done towards
and Alhajri, Elkamel, Albahri, and Douglas (2008). Although com- the integration of planning and crude-oil scheduling problems in
mercial solvers such as GRTMPS (Haverly Systems), PIMS (Aspen the context of refineries. To the best of our knowledge, Produc-
Tech), and RPMS (Honeywell Hi-Spec Solutions) implement in- tion Scheduler (Honeywell Hi-Spec Solutions) is the only existing
house successive linear programming (SLP) algorithms to solve this commercial tool that integrates simultaneously both problems.
problem (see Zhang et al., 1985), any standard NLP solvers can also The main difficulty arises from the fact that, in this particu-
be used although they may not guarantee global optimality of the lar case, the planning model is not an aggregate scheduling
model. Thus, the decomposition methods developed for batch
and continuous plants are not directly applicable to refineries.
∗ Corresponding author. Tel.: +1 412 268 2230; fax: +1 412 268 7139. In particular, refinery planning and crude scheduling correspond
E-mail address: grossmann@cmu.edu (I.E. Grossmann). to two different problems solely linked through CDU feedstocks.

0098-1354/$ – see front matter © 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/j.compchemeng.2011.03.026
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2751

Fig. 1. Basic refinery planning system.

Therefore, instead of using a hierarchical decomposition, a spa- 2.1. Refinery planning problem
tial Lagrangian decomposition is preferred. The reader may refer
to Fisher (1985) and Guignard (2003) for extensive reviews The refinery planning problem can be regarded as a flowsheet
on Lagrangian relaxation and decomposition techniques. These optimization problem with multiple periods during which the
approches have been applied to many industrial problems such refinery system is assumed to operate in steady-state. Due to exten-
as production planning and scheduling integration (see Li & sive stream mixing, the model for each period is based on a pooling
Ierapetritou, 2010) or multiperiod refinery planning (see Neiro & problem that is extended in order to include process models for
Pinto, 2006). each refining unit. The different periods in the model are con-
The content of this paper is organized as follows. The planning nected through material inventories, which are usually constrained
and scheduling problems are stated in Section 2 as well as the by storage limitations. These inventories allow product transfers
full-space integrated problem. In Section 3, a Lagrangian decom- between periods and final stock optimization. In this work, we con-
position scheme based on the dualization of CDU feedstock linking sider a single-period planning model based on a pooling problem
constraints is presented. A new hybrid method is introduced to inspired from the literature (see for instance Adhya, Tawarmalani,
solve the Lagrangian relaxation in Section 4. A heuristic algorithm & Sahinidis, 1999). A basic refinery planning system is represented
is developed to obtain good feasible solutions for the integrated in Fig. 1. A set of crudes i ∈ I are to be mixed in different types of
full-space problem in Section 5. In Section 6, we provide technical crude-oil blends j ∈ J (e.g. low-sulfur and high-sulfur blends), each
remarks in order to further develop the practical details of the pro- associated to a specific CDU operating mode. For each mode and
posed approach. Numerical illustrations of the proposed approach each crude, several distillation cuts are obtained with different
are presented in Section 7 for the case study and Section 8 for a fixed yields. These crude cuts are then blended into intermediate
larger refinery problem. Section 9 concludes the paper. pools which are used to prepare several final products. Therefore,
the refinery planning system is composed of the following ele-
ments:

2. Problem statement • One input stream for each selected crude i ∈ I and each type of
crude blend j ∈ J
In this section, we present the main concepts involved in the • One CDU with fixed yields for each distillation cut
refinery planning and crude-oil scheduling problems. The full- • Set of distillation cuts k ∈ K
space model is detailed as well as from case study examples. • One pool for each type of crude blend j ∈ J and each cut k ∈ K

Fig. 2. Refinery planning case study.

2752 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

• One intermediate stream between each pool (j, k) ∈ J × K and each Dl is the maximum demand in product l
final product l ∈ L Zlp is the maximum specification for quality p of product l
• Set of final products l ∈ L
• One sales stream for each final product l ∈ L Fig. 2 displays the pooling structure of a case study with corre-
• Set of qualities p ∈ P sponding data for crudes (A, B), blends (X, Y), distillation cuts (M, N),
and final products (P, Q, R, S). Two different qualities are considered
The yield of crude i ∈ I in distillation cut k ∈ K, when processed in this case.
in crude blend j ∈ J is assumed to be fixed and is denoted by ˛ijk . In In the remainder of the paper, we consider the following NLP,
terms of stream qualities, it is assumed that distillation cuts have which is a simplified version of the planning model (PP ).
fixed qualities while pool qualities are calculated by bilinear quality ⎧ T
balance constraints. A pure flow-based model is used to formulate ⎪
⎨ max VP xS
the pooling problem as shown below. CDU flowrate limitations are s.t. fP (xF , xI , xS ) ≤ 0
(PP )
considered independent of the operating mode and are enforced ⎪
⎩ gP (xI ) ≤ 0
globally for all crudes processed during the period. xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S|

max pl xSl (sales revenue maximization)
l∈L

 ij
s.t. 0≤ xF ≤ C i i∈I (crude availability)
j∈J

 ij
FR · H ≤ xF ≤ FR · H (CDU flowrate limitations)
i∈I j∈J

ijk ij
x1 = ˛ijk · xF (i, j, k) ∈ I × J × K (CDU yield calculation)
 ijk
 jkl
x1 = x2 (j, k) ∈ J × K (pool mass balance)
i∈I l∈L
 ijkp ijk jkp
 jkl
q1 x1 = q2 x2 (j, k, p) ∈ J × K × P (pool quality balance, nonlinear)
i∈I l∈L
 jkl
x2 = xSl l∈L (product mass balance)
j∈J k∈K

xSl ≤ Dl l∈L (maximum product demand)

 jkp jkl
q2 x2 ≤ Z lp xSl (l, p) ∈ L × P (product quality requirement, nonlinear)
j∈J k∈K

ij ijk jkl jkp

The nomenclature used is as follows:
xF , x1 , x2 , xSl ≥ 0, q2 ∈ R
• VP is the market value of final products
The nomenclature used is as follows:
• xF is a set of continuous variables representing CDU feedstock
• Variables: quantities over the single planning period
• xS is a set of continuous variables representing final products sales
ij • xI is a set of intermediate continuous variables (e.g. pool quantity
xF is a variable representing the amount of crude i selected for CDU
distillation in blend j and quality variables)
ijk • fP (xF , xI , xS ) ≤ 0 is the set of linear constraints (e.g. material balance
x1 is a variable representing the amount of cut k extracted from
constraints)
crude i in blend j
jkl • gP (xI ) ≤ 0 is the set of nonlinear constraints (e.g. quality balance
x2 is a variable representing the flow of material between pool (j,
constraints)
k) and product l
jkp
q2 is a variable representing the quality p of pool (j, k) 2.2. Crude-oil scheduling problem
xSl is a variable representing the amount of final product l sold
The crude-oil scheduling problem deals with the unloading,
• Parameters: transfer and blending operations executed on crude-oil tankers and
crude-oil inventories. The goal is to sequentially prepare multiple
ijkp
q1 is the (fixed) quality p of cut k extracted from crude i in blend crude blends, which are defined by specific property requirements.
j Each type of crude blend corresponds to a specific CDU operating
p is the market value of final products mode. Different objectives have been studied, namely minimiza-
Ci is the amount of available crude i tion of logistics costs (see Lee et al., 1966) or maximization of profit
˛ijk is the yield of cut k extracted from crude i in blend j (see Mouret, Grossmann, & Pestiaux, 2009; Mouret, Grossmann, &
H is the planning horizon Pestiaux, 2010). In this work, the objective is to minimize the total
[FR, FR] is the bounds on CDU flowrate replacement cost of the crudes that are selected for distillation. The
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2753

Table 1
Crude-oil scheduling data for case study.

Scheduling horizon 8 days

Vessels Arrival time Composition Amount of crude

Vessel 1 0 100% A 100
Vessel 2 4 100% B 100

Storage tanks Capacity Initial composition Initial amount of crude

Tank 1 [0, 100] 100% A 25
Tank 2 [0, 100] 100% B 75

Charging tanks Capacity Initial composition Initial amount of crude

Tank 1 (mix X) [0, 100] 80% A, 20% B 50
Tank 2 (mix Y) [0, 100] 20% A, 80% B 50

Crudes Property 1 (sulfur concentration) Crude unit cost

Crude A 0.01 7
Crude B 0.06 6

Crude mixtures Property 1 (sulfur concentration) Maximum number of batches

Crude blend X [0.015, 0.025] 4
Crude blend Y [0.045, 0.055] 4

Unloading flowrate [0, 50] Transfer flowrate [0, 50]

Distillation flowrate [5, 50] Minimum duration of distillations 1 day

In the remainder of the paper, we consider the following MINLP,

which is a simplified version of the scheduling model PS .
⎧
⎪ max −VCT yF
⎨
s.t. fS (yB , yC , yF ) ≤ 0
(PS )
⎪ gS (yC ) ≤ 0
⎩
yB ∈ {0, 1}|B| , yC ∈ R|C| , yF ∈ R|F|
Fig. 3. Crude-oil scheduling system 1 (Lee et al., 1996). The nomenclature used is as follows:

replacement cost is the cost of replacing the crude once it has been • VC is the replacement cost of crude-oils (usually based on market
processed. The crude-oil schedule must satisfy inventory capac- value)
ity limitations, crude tankers arrival dates as well as the following • yF is a set of continuous variables representing total CDU feed-
logistics constraints: stock quantities over the scheduling horizon
• yC is a set of continuous variables representing other continuous
(i) only one berth is available at the docking station for crude decisions (e.g. timing decisions)
tanker unloadings, • yB is a set of binary variables representing sequencing decisions
(ii) inlet and outlet transfers on tanks must not overlap, (see Mouret et al., 2011)
(iii) a tank may charge only one CDU at a time, • fS (yB , yC , yF ) ≤ 0 is the set of linear constraints (e.g. scheduling
(iv) a CDU can be charged by only one tank at a time, constraints)
(v) CDUs must be operated continuously throughout the schedul- • gS (yC ) ≤ 0 is the set of nonlinear stream composition constraints
ing horizon. (see Mouret et al., 2009)

Fig. 3 shows the refinery system corresponding to problem 1 2.3. Full-space problem
introduced in Lee et al. (1996). Table 1 displays the dimensionless
data for this example. Besides a different objective function, the Given the importance of crude selection for refinery optimiza-
example is modified by introducing a minimum duration of 1 day tion, the refinery planning problem and the crude-oil scheduling
for distillation operations. Therefore, due to crude blend alternative problem should ideally be optimized simultaneously. This can only
sequencing, at most 4 batches of each crude mix can be processed be done by solving an integrated full-space MINLP problem, denoted
in 8 days. The scheduling problem is formulated using the Multi- (P), which aims at optimizing all refinery decisions subject to plan-
Operation Sequencing (MOS) time representation as introduced in ning, scheduling, and linking constraints.
Mouret et al. (2010) (see Appendix B). ⎧
⎪ max VPT xS − VCT yF
⎪
⎪
⎪
⎪ s.t. fP (xF , xI , xS ) ≤ 0
⎪
⎪ gP (xI ) ≤ 0
⎪
⎨
fS (yB , yC , yF ) ≤ 0
(P)
⎪ gS (yC ) ≤ 0
⎪
⎪
⎪
⎪ yF − xF = 0
⎪
⎪ xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S|
⎪
⎩
yB ∈ {0, 1}|B| , yC ∈ R|C| , yF ∈ R|F|

The integrated objective is to maximize profit defined by final

product sales revenues minus crude-oil replacement costs. The
linking constraint yF − xF = 0 ensures consistency of planning and
scheduling decisions in terms of CDU feedstock quantities. More
Fig. 4. Economic interpretation of the Lagrangian decomposition. precisely, it ensures that the amounts of crudes selected for distil-
2754 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

Fig. 5. General iterative primal-dual algorithm.

Fig. 6. (a and b) Plots of the feasible space of (P̂DK+1 ).

lation are identical in the planning and scheduling solutions. Also, Robertson, Palazoglu, and Romagnoli (2010) proposed a multi-
to be consistent in time, it is considered that the planning and level approach consisting of approximating the refinery planning
scheduling horizons have identical lengths. model by multiple linear regressions that are then used in the
In the remainder of the paper, we use the notation v(P) to denote crude-oil scheduling problem for the minimization of the total
the optimal objective value for problem (P). logistics and production costs. The method is applied to a case
study comprising two different crudes. Although computationally
effective, the use of linear regressions may not be sufficient for
3. Lagrangian decomposition scheme
the the global optimization of highly nonlinear refinery planning
models.
The full-space problem (P) is a large-scale MINLP, which con-
In this work, we present an integration approach based on
tains many binary variables from the crude-oil scheduling problem
Lagrangian decomposition, which is a special case of Lagrangian
and many non-convex constraints from the refinery planning
relaxation (Guignard, 2003). The idea is to build a relaxed version
model. Due to convergence issues and the presence of many poten-
of the full-space problem, which is decomposable, and therefore,
tial local optima, standard MINLP solvers for convex optimization,
much easier to solve. In particular, the decomposition procedure
such as AlphaECP, Bonmin, DICOPT, KNITRO, or SBB, may fail solv-
is based on the dualization of the linking constraint yF − xF = 0. The
ing the model or return poor solutions. Global MINLP solvers, such
relaxed problem (PR ()), composed of NLP and MINLP models, is
as BARON, Couenne, or LINDOGLOBAL, are in principle able to solve
defined by removing this constraint and penalizing its violations by
the problem but they may require prohibitive computational times.
adding the term T (yF − xF ) to the objective function. The parame-
Therefore, a specific solution strategy needs to be developed to
address this problem.

Fig. 7. Iterative primal-dual algorithm with heuristic step.

 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2755

the scheduling system. Therefore, the spatial Lagrange decomposi-

tion procedure applied to this problem can be seen as introducing
a crude market between the planning and scheduling systems (see
Fig. 4). The planning system acts as a consumer who buys crude
from the market, while the scheduling system acts as a producer
who sells it to the market. It is clear that for fixed prices both actors
are independent, which explains why the two corresponding sub-
problems can be solved in parallel.
Although computationally convenient, this decomposition pro-
cedure does not solve the original full-space problem. In particular,
it is well-known that (PR ()) is a relaxation of the full-space prob-
lem, thus v(PR ()) > v(P). However, one can search for a Lagrange
Fig. 8. Crude-oil scheduling and multi-period refinery planning integration. multiplier  that minimizes v(PR ()) in order to get as close as possi-
ble to v(P). This problem is called the dual problem and the function
ter  is a Lagrange multiplier whose value is fixed prior to solving  → v(PR ()) is often called the Lagrangian function.
the model and adjusted iteratively.
(PD ) min v(PR ())
⎧ T 
⎪ max VPT xS − T xF + ( − VC ) yF
⎪
⎪
⎪
⎪ s.t. fP (xF , xI , xS ) ≤ 0 Duality theory establishes that v(PD ) − v(P) ≥ 0 (see Guignard,
⎪
⎨ gP (xI ) ≤ 0 2003). In some cases (e.g. nonconvex models), we may have
(PR ()) fS (yB , yC , yF ) ≤ 0 v(PD ) − v(P) > 0 and this difference is called dual gap. Our hope
⎪
⎪ gS (yC ) ≤ 0
⎪
⎪ is that this dual gap is small enough so that valuable information
⎪
⎪ xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S| can be inferred to generate near-optimal heuristic solutions (see
⎩
yB ∈ {0, 1}|B| , yC ∈ R|C| , yF ∈ R|F| Section 5).

As already mentioned, problem (PR ()) is easier to solve as it can

4. Solution of the dual problem
be decomposed into two subproblems (PP ()) and (PS ()).

v(PR ()) = v(PP ()) + v(PS ()) Several approaches have been proposed in the literature in order
to solve the dual problem associated with the Lagrangian relax-
The subproblem (PP ()), an NLP, is a modification of the original
ation. A classical approach is the subgradient method proposed
refinery planning problem (PP ) as it consists of assigning crude costs
by Held and Karp (1971, 1974). Many researchers have used and
 to the CDU feedstock variables xF . For a given crude i, increasing i
improved this technique over the years (see Bazaraa and Sherali,
will decrease the incentive to select this crude for distillation pro-
1981; Camerini et al., 1975; Fisher, 1981). This approach is pre-
cessing. On the other hand, decreasing i will increase the incentive
ferred as it usually predicts very good Lagrange multiplier updates.
to select it.
⎧ However, special care must be taken in order to insure conver-
T T
⎪
⎨ max VP xS −  xF gence and it requires a good strategy for defining and updating the
s.t. fP (xF , xI , xS ) ≤ 0 subgradient step size. Another approach that theoretically displays
(PP ())
⎪
⎩ gP (xI ) ≤ 0 better convergence properties has been introduced by Cheney and
xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S| Goldstein (1959) and Kelley (1960). It is often denoted as the cutting
plane method. In practice, this approach usually takes a long time
The subproblem (PS ()), an MINLP, is a modification of the orig- to converge as many iterations are necessary in order to obtain
inal crude-oil scheduling problem (PS ) as it consists of assigning good Lagrange multiplier updates. A refinement of this approach is
crude values  to the CDU feedstock variables yF . For a given crude the boxstep method presented in Marsten, Hogan, and Blankenship
i, increasing i will increase the incentive to select this crude for (1975). It allows obtaining better updates for the Lagrange multi-
blending and distillation processing. On the other hand, decreasing plier during early iterations while keeping the same convergence
i will decrease the incentive to select it. properties. Other refinements of previous approaches include the
⎧ T
⎪
⎨ max ( − VC ) yF
bundle method (Lemaréchal, 1974), the volume algorithm (Barahona
s.t. fS (yB , yC , yF ) ≤ 0 & Anbil, 2000), and the analytic center cutting plane method (Goffin,
(PS ())
⎪
⎩ gS (yC ) ≤ 0 Haurie, & Vial, 1992).
yB ∈ {0, 1}|B| , yC ∈ R|C| , yF ∈ R|F| All the above approaches are based on an iterative solution
procedure between the primal and the dual world. Fig. 5 gives a
On the whole, the Lagrange multiplier  can be seen as a crude schematic description of this algorithm. The first step consists of
purchase cost for the planning system, and as a crude sales value for initializing the Lagrange multipliers. Problem-specific strategies,
often based on the economic interpretation of , exist in order
to provide good initial values. Then, at each iteration the relaxed
problem is solved and a primal solution is obtained. If a stopping cri-
terion is satisfied, the algorithm converges. Otherwise, the Lagrange
multipliers are updated for the next iteration. The definition of the
stopping criterion depends on the approach used and can, in certain
cases, ensure convergence to an optimal dual solution ∗ .
In this work, we introduce a new hybrid method to update the
Lagrange multipliers. It is based on the three concepts of cutting
planes, subgradient and boxstep. Cutting planes are valid con-
straints for the dual problem that are generated at each iteration.
They are used to record valuable dual information to be used dur-
Fig. 9. Disaggregated CDU feedstocks synchronization. ing later iterations. A subgradient provides a descent direction for
2756 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

Fig. 10. Complete algorithm implementation.

Table 2
Lagrangian iterations statistics (6 priority-slots, NLP = SNOPT).

Iteration Pure dual Hybrid dual Step size Modified relaxation Modified heuristic Original heuristic CPU time
v(PDK ) v(P̂DK ) ˛K v(P̃R (K )) v(P̃H (yBK )) v(PH (yBK ))
1 –a –b –b 645.000 400.942 –c 3s
2 –a 389.942 1 689.049 564.000 560.500 9s
3 –a 547.929 1 637.063 592.368 592.368 13 s
4 –a 582.790 1 603.377 592.368 592.368 18 s
5 –a 591.172 1 609.526 586.191 –c 24 s
6 –a 590.629 0.851 598.380 562.881 –c 29 s
7 580.717 591.324 −0.617 600.171 586.191 –a 33 s
8 592.368 592.785 0.390 594.774 586.191 –a 46 s
9 592.369 592.369 1 592.372 592.368 592.368 50 s
a
Unbounded LP.
b
Not available.
c
Locally infeasible NLP.

Table 3
Lagrangian iterations statistics (7 priority-slots, NLP = SNOPT).

Iteration Pure dual Hybrid dual Step size Modified relaxation Modified heuristic Original heuristic CPU time
v(PDK ) v(P̂DK ) ˛K v(P̃R (K )) v(P̃H (yBK )) v(PH (yBK ))
1 – – – 645.000 393.470 – 3s
2 – 381.562 1 751.494 568.077 568.077 10 s
3 – 552.076 1 622.785 592.368 – 16 s
4 – 574.735 1 614.178 592.368 592.368 22 s
5 – 583.166 1 602.218 592.368 592.079 50 s
6 – 588.279 1 617.268 592.368 592.079 56 s
7 – 592.856 0.911 600.752 592.368 592.368 66 s
8 – 592.835 0.517 595.264 592.368 – 81 s
9 – 592.288 0.357 595.292 592.368 592.368 95 s
10 592.368 592.369 1 592.369 592.368 592.079 101 s
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2757

Fig. 11. Lagrangian iteration objective values (6 priority-slots, NLP = SNOPT). Fig. 13. Lagrange multiplier updates (6 priority-slots, NLP = SNOPT).

the dual problem, while a boxstep is defined to allow deviations by the parameter ˛ ¯ > 0. Note that ˛ is allowed to take negative
from this direction within a specified domain. The combination values. Variable ı defines a deviation from the subgradient step.
of these techniques ensures good convergence properties, while The parameter ı̄ > 0 is the maximum deviation in each multiplier
providing efficient Lagrange multiplier updates. At iteration K + 1, direction and defines a full-dimensional boxstep. Both subgradient
the Lagrange multiplier is updated to the solution of the following and boxstep concepts are simultaneously embedded in constraint
restricted LP dual problem with subgradient-based boxstep. (SG+TR). Note that the parameters ˛ ¯ and ı̄ can be heuristically
⎧ updated at each iteration. In practice, our computational experi-
⎪
⎪
min 
ments have shown that using fixed values for ˛ ¯ and ı̄ is a reasonable
⎪
⎨ s.t.  ≥ VPT xSk − VCT yFk + T (yFk − xFk ) ∀k = 1 . . . K (CP k )
v(PD ) − v(PR (K )) K choice.
(P̂DK+1 )  = K + ˛ (yF − xFK ) + ı (SG + TR)
⎪ ||yFK − xFK ||2 Fig. 6 displays the feasible space (grey area) of (P̂DK+1 ). The pro-
⎪
⎪  ∈ R,  ∈ R |F|
⎩ |F|
jection on the space of Lagrange multipliers (1 , 2 ) depicts the
˛ ∈ ] − ∞, ˛]
¯ , ı ∈ [−ı̄, ı̄ shape of the subgradient-based boxstep. In the space of (1 , ), CPk
The variables  and  are classically used in the pure cutting represents the projection of the cutting plane generated at itera-
plane approach. The pure cutting plane restricted LP dual prob- tion k. Note that the feasible space of (P̂DK+1 ) contains (1 = K1 , 2 =
lem (PDK+1 ) consists of minimizing  subject to constraints (CPk ), K2 ,  = v(PR (K1 , K2 ))). In both plots (a) and (b), 1 corresponds to a
k = 1 . . . K only: lower bound on multiplier 1 induced by the boxstep constraints.
 The stopping criterion for this hybrid strategy is identical to
min  the pure cutting plane method and is based on the Lagrangian
(PDK+1 ) s.t.  ≥ VPT xSk − VCT yFk + T (yFk − xFk ) ∀k = 1 . . . K (CP k ) gap between the relaxed primal problem and the restricted dual
 ∈ R,  ∈ R|F| problem:
This problem is always unbounded during early iterations (i.e. v(PR (K )) − v(PDK ) ≤ ε (1)
v(PDK ) = −∞), so it cannot be used directly to update the Lagrange
multipliers. This issue can be solved by defining a bounded feasi- In the pure cutting plane approach (without constraint (SG+TR)),
ble set for the multipliers based on their interpretation (e.g. lower the optimal value of the restricted dual problem v(PDK ) iteratively
and upper bounds). However, in this work, an iteratively updated approximates v(PD ) from below since it involves the minimization
boxstep based on the subgradient step is used so that the restricted of a relaxation of v(PD ).
dual problem (P̂DK+1 ) is bounded. v(PDK ) ≤ v(PD ) ∀K (2)
The subgradient step is classically defined as  = K + ˛(v(PD ) −
v(PR (K )))/||yFK − xFK ||2 (yFK − xFK ) where v(PD ) can be estimated using Therefore, the stopping criterion (1) ensures convergence to an
a heuristic solution for (P). However, instead of heuristically updat- ε-optimal solution of the dual problem (PD ) as:
ing the step size, it is optimized using variable ˛, which is bounded
v(PR (K )) − v(PDK ) ≤ ε ⇒ v(PR (K )) − v(PD ) ≤ ε (3)

In the proposed approach, when the pure restricted dual prob-

lem (PDK ) becomes bounded, the stopping criterion (1) is used to
check convergence while the hybrid restricted dual problem (P̂DK ) is
used to update the Lagrange multipliers. Finite convergence prop-
erties for the pure cutting plane and boxstep methods in the context
of mixed-integer linear programming are reported in Frangioni
(2005) and Marsten et al. (1975), respectively. For the rest of the
paper, it is assumed that practical convergence of the proposed
hybrid method can be achieved in the context of integrated refinery
planning and scheduling.

5. Heuristic solutions

In this section, a classical adaptation of the primal-dual algo-

rithm is presented in order to obtain solutions that satisfy all
Fig. 12. Lagrangian iteration objective values (7 priority-slots, NLP = SNOPT). constraints of the full-space problem, including linking constraints.
2758 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

Table 4
Comparative performance of several MINLP algorithms.

6 priority-slots 7 priority-slots

MINLP solver Objective value CPU time Optimality gap Objective value CPU time Optimality gap

Proposed (CONOPT) 592.368 37 s 0% 592.368 94 s 0%

Proposed (SNOPT) 592.368 50 s 0% 592.368 101 s 0%
Proposed (IPOPT) 592.368 244 s 0% 592.368 833 s 0%

Sequential (BARON) 545.000 9s –b 545.000 10 s –b

DICOPT (CONOPT) 545.000 5s –b 592.368 7s –b

DICOPT (SNOPT) 592.368 429 s –b 592.368 6s –b
DICOPT (IPOPT) 568.077 54 s –b 592.368 44 s –b

AlphaECP (CONOPT) 512.324 67 s –b 545.000 120 s –b

AlphaECP (SNOPT) 512.324 67 s –b 545.000 395 s –b
AlphaECP (IPOPT) 512.324 69 s –b 545.000 175 s –b

SBB (CONOPT) 592.368 267 s –b 592.368 +1000 s –b

SBB (SNOPT) –a +1000 s –b –a +1000 s –b
SBB (IPOPT) –a +1000 s –b 493.536 +1000 s –b

LINDOGLOBAL 568.077 +1000 s 11.9% 532.857 +1000 s 17.1%

BARON 592.170 +1000 s 7.3% 400.000 +1000 s 37.5%
a
No solution found.
b
Not available.

The heuristic algorithm consists of fixing binary variables yB

from the crude-oil scheduling formulation to their values yBK in the
solution of the relaxed problem (PR (K )). As a consequence, the full-
space problem (P) reduces to a continuous NLP, denoted (PH (yBK )),
and can then be efficiently solved. If it is feasible and its (local) opti-
mal solution is better than the previous incumbent, PLB is updated.
Otherwise, PLB is left unchanged.

⎧
⎪ max VPT xS − VCT yF
⎪
⎪
⎪
⎪ s.t. fP (xF , xI , xS ) ≤ 0
⎪
⎪
⎪
⎨ gP (xI ) ≤ 0
K
fS (yB , yC , yF ) ≤ 0
(PH (yBK ))
⎪
⎪ gS (yC ) ≤ 0
⎪
⎪ yF − xF = 0
⎪
⎪
⎪
⎪ xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S|
Fig. 14. Lagrange multiplier updates (7 priority-slots, NLP = SNOPT).
⎩
yC ∈ R|C| , yF ∈ R|F|

As explained by Frangioni (2005), the solution of the Lagrangian

dual problem yields primal information that can be used to gener- In this heuristic solution, fixing binary variables yB to yBK cor-
ate good heuristic solutions for (P). In this paper, a heuristic step is responds to fixing the selection and sequencing of operations for
introduced in the iterative algorithm to produce valid lower bounds the crude-oil scheduling system. Therefore, when solving problem
PLB (see Fig. 7). This induces a second stopping criterion based on the (PH (yBK )), the nonlinear solver has the opportunity to re-optimize
dual gap: v(PR (K )) − P LB ≤ ε. If either of the two stopping criteria all other continuous decisions such as quantities, blend recipes,
is satisfied, the algorithm converges and returns PLB . and timing of operations (start time, duration, and end time). It
is crucial to note that the timing of operations can only be re-
optimized if these decisions are handled by continuous variables.
For instance, in Mouret et al. (2010), the discrete-time represen-
tation, denoted by MOS-FST, uses binary variables to determine
the timing decisions whereas all other representations, denoted
by MOS, MOS-SST, and SOS, use continuous variables instead. This
shows a clear advantage, although not intuitive, of continuous-
time scheduling formulations over discrete-time representations.
In particular, the latter might be inefficient in the context of this
work as it would decrease the flexibility of the heuristic algorithm
to find good, or at least feasible solutions.
Overall, it is interesting to note that using this approach, the
iterative primal-dual algorithm that solves the Lagrangian dual
problem acts as a discrete solution generator that suggests poten-
tially good discrete solutions for the full-space problem. In other
words, it searches the optimal selection and sequencing of opera-
tions for the crude-oil scheduling system. Additionally, it provides
Fig. 15. Blend compositions in solutions obtained with 6 priority-slots. an upper bound for the global optimal solution.
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2759

Fig. 16. Planning model for larger refinery problem.

6. Remarks either exchange crudes between the two systems or sell and buy to
and from the crude market (see Fig. 4). From this observation, it is
6.1. CDU feedstocks and Lagrange multipliers natural to use the crude costs defined in the crude-oil scheduling
problem as initialization for the Lagrange multipliers. For the case
The number of Lagrange multipliers highly depends on the study, we use the following initial values:
CDU feedstock possibilities. Exactly one multiplier is needed
for each feasible combination of crude and type of crude 1(X,A) = 1(Y,A) =7
blend (or corresponding CDU operating mode). In the case 1(X,B) = 1(Y,B) =6
study, there are 2 different crudes which can both be blended
in any of the 2 different types of blends, so 4 Lagrange
multipliers are needed to solve the dual problem. Typical 6.2. Multi-period refinery planning
large-scale refineries may need 50 and up to 100 Lagrange
multipliers. In this work, the refinery planning problem is expressed
The optimal value of the Lagrange multipliers correspond to the over a single period for which CDU feedstocks are synchro-
optimal marginal costs of the linking constraints for the convexified nized with the crude-oil scheduling problem. Even though it is
problem (for more details in the case of MILP models, see Frangioni, often computationally critical to increase the time horizon for
2005). Therefore, the Lagrange multipliers can be seen as the opti- scheduling problem, this can easily be done in refinery plan-
mal pricing strategy between the crude-oil scheduling system and ning models by introducing additional time periods. Therefore,
the refinery planning system. In other words, the optimal Lagrange one can define 2 or more time periods for the refinery plan-
multipliers are crude prices for which it is marginally equivalent to ning model and synchronize CDU feedstocks for the first period
only, as shown in Fig. 8. In this particular case, the refin-
ery planning decisions for the second period, including CDU
Table 5
Crude cut prices and specification for larger refinery problem. feedstock decisions xF , are made without taking into account
crude-oil scheduling constraints. This methodology allows making
Crude cut Price Specifications

LPG 8.5 None

Naphta 8.0 Specific gravity ∈[0.72, 0.775]

Motor octane number ≥45
Research octane number ≥45
Sulfur weight content ≤120 ppm

Kerosene 7.0 Specific gravity ∈[0.775, 0.84]

Freeze point ≤−40 ◦ C

Diesel 8.0 Specific gravity ∈[0.82, 0.86]

Cetane number ≥48
Cloud point ≤4 ◦ C
Sulfur weight content ≤2800 ppm

Residue 6.5 None Fig. 17. Crude-oil scheduling system 3 (Lee et al., 1996).
2760 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

Table 6
Crude-oil scheduling data for larger refinery problem.

Scheduling horizon 12 days

Vessels Arrival time Composition Amount of crude

Vessel 1 0 100% A 500
Vessel 2 4 100% B 500
Vessel 3 8 100% C 500

Storage tanks Capacity Initial composition Initial amount of crude

Tank 1 [0, 1000] 100% D 200
Tank 2 [0, 1000] 100% E 200
Tank 3 [0, 1000] 100% F 200

Charging tanks Capacity Initial composition Initial amount of crude

Tank 1 (mix 1) [0, 1000] 100% G 300
Tank 2 (mix 2) [0, 1000] 100% E 500
Tank 3 (mix 3) [0, 1000] 100% F 300

Crudes Property 1 (sulfur concentration) Crude unit cost

Crude A 0.01 6
Crude B 0.085 6.5
Crude C 0.06 5.5
Crude D 0.02 7.2
Crude E 0.05 6.7
Crude F 0.08 6.2
Crude G 0.03 7.5

Crude mixtures Property 1 (sulfur concentration) Maximum number of batches

Crude mix 1 [0.025, 0.035] 6
Crude mix 2 [0.045, 0.065] 6
Crude mix 3 [0.075, 0.085] 6

Unloading flowrate [0, 50] Transfer flowrate [0, 50]

Distillation flowrate [5, 50] Minimum duration of distillations 1 day

short-term scheduling decision while considering the long-term feedstocks. Fig. 9 depicts the synchronization of disaggregated CDU
economic impacts of these decisions, which cannot be done with batches.
detailed long-term scheduling models due to their computational
complexity. 6.4. Handling nonlinearities in crude-oil scheduling model

6.3. CDU feedstocks aggregation Although, the relaxed problem (PR ()) is decomposable, it is
not easy to solve it to global optimality. In particular, two major
An important issue with the proposed approach comes from the issues arise. First, the crude-oil scheduling problem (PS ()) corre-
fact that the CDU feedstocks for the linking period are aggregated. sponds to an MINLP due to the presence of nonlinear composition
In the optimal solution, the crude-oil operations schedule pre- constraints. In Mouret et al. (2009) an MILP relaxation is derived
pares several batches for each type of crude blends. Then, for each by simply dropping these nonlinear constraints. The solution can
of these blend types, all the corresponding batches are accumu- then be refined by fixing the binary variables and solving the
lated and the refinery planning solution determines the processing reduced NLP, similarly to the heuristic approach presented in Sec-
decisions for the aggregated batch. This approximation may lead tion 5. Results show that the solution obtained is close to the global
to sub-optimality or even technical infeasibility of the solutions optimum as it tends to satisfy the relaxed nonlinear constraints.
obtained. This problem can be solved by postulating exactly one Therefore, a similar methodology is used in this work. Instead of
period for each batch and synchronizing all the corresponding CDU simply dualizing the linking constraints, the nonlinear scheduling

Table 7
Lagrangian iterations statistics for larger refinery problem (6 priority-slots, NLP = CONOPT).

Iteration Pure dual Hybrid dual Step size Modified relaxation Modified heuristic Original heuristic CPU time
v(PDK ) v(P̂DK ) ˛K v(P̃R (K )) v(P̃H (yBK )) v(PH (yBK ))
1 – – – 266.226 – – 43 s
2 – −2.166 1 417.311 257.943 222.339 102 s
3 – 253.546 1 411.657 258.488 245.306 173 s
4 – 281.076 0.857 325.906 258.106 – 235 s
5 – 260.746 1 288.587 258.104 – 298 s
6 – 256.301 1 273.197 258.516 245.467 372 s
7 – 256.385 1 271.415 258.027 250.989 459 s
8 – 258.091 1 264.974 257.961 233.713 539 s
9 – 258.847 0.861 265.149 257.961 233.713 610 s
10 – 257.623 1 267.541 257.740 228.230 682 s
11 – 259.500 0.868 261.237 257.750 247.971 769 s
12 – 259.170 1 261.512 257.750 247.971 850 s
13 – 259.126 0.759 264.109 258.027 250.989 913 s
14 – 259.121 1 263.708 258.306 – 979 s
15 – 259.592 0.774 260.857 258.516 238.763 1045 s
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2761

Table 8 6.5. Handling nonlinearities in the refinery planning model

Optimal Lagrange multipliers for each crude and each CDU mode.

Crude CDU mode Initial value (=Price) Optimal value In order to obtain valid upper bounds when solving (PR ()) or
1 6.0 6.169 (P̃R ()), the refinery planning problem (PP ()) has to be solved to
A 2 6.0 6.207 global optimality. Although global optimization of industrial large-
3 6.0 7.695 scale refinery planning models is still not achievable, the refinery
1 6.5 6.971 planning case study presented in Section 2 is solvable by the global
B 2 6.5 6.987 NLP solver BARON in reasonable time. However, it should be noted
3 6.5 7.026 that it is critical to provide tight bounds for the quality variables
jkp
1 5.5 6.062 q2 . In particular, based on the structure of the pooling system, we
C 2 5.5 6.009 use the following bounds:
3 5.5 6.054
ijkp jkp ijkp
1 7.2 7.186 minq1 ≤ q2 ≤ maxi ∈ I q1
i∈I
D 2 7.2 7.262
3 7.2 7.962

1 6.7 6.859 6.6. Detailed implementation

E 2 6.7 6.878
3 6.7 6.898 Based on previous remarks to handle nonlinearities, the com-
1 6.2 6.340 plete heuristic algorithm is developed as depicted in Fig. 10.
F 2 6.2 6.226 Although global optimality cannot be ensured, the dual gap can
3 6.2 6.293 be estimated using the upper bound provided by v(P̃R (∗ )). Local
1 7.5 6.930 NLP solvers, such as CONOPT, are used for the heuristic steps as
G 2 7.5 7.360 the solution time is more critical than global optimality for these
3 7.5 7.360 problems. The hybrid restricted dual problem (P̂DK ) is solved using
the best solution of the modified heuristic problems (P̃H (yBK )) to
estimate the optimal dual solution v(PD ). The stopping criterion is
constraints gS (yC ) ≤ 0 are also relaxed (i.e. dropped). The modi- based on relative gaps but can also be expressed in terms of abso-
fied relaxed full-space MINLP problem (NLP + MILP) is denoted by lute gaps. Converging on the Lagrangian gap means that no further
(P̃R ()): improvement of the upper bound can be achieved and the current
⎧ T Lagrange multipliers are optimal.
⎪ max VPT xS − T xF + ( − VC ) yF
⎪
⎪ s.t. fP (xF , xI , xS ) ≤ 0
⎪
⎨ 6.7. Handling infeasibilities
gP (xI ) ≤ 0
(P̃R ())
⎪
⎪ fS (yB , yC , yF ) ≤ 0
⎪
⎪ xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S|
In this section, we discuss potential infeasibility issues in the
⎩ three subproblems (planning, scheduling and heuristic). First, the
yB ∈ {0, 1}|B| , yC ∈ R|C| , yF ∈ R|F| NLP solver may fail to find a feasible solution to the non-convex
refinery planning problem. In this case, no rigorous conclusion can
The corresponding modified crude-oil scheduling MILP sub-
be drawn on the existence of feasible solution for this problem. To
problem is denoted by (P̃S ()):
solve this issue, one may use a multi-start strategy or any global-
 T ization techniques for non-convex NLPs. Another possible strategy
max ( − VC ) yF
(P̃S ()) s.t. fS (yB , yC , yF ) ≤ 0 consists in perturbing the Lagrange multipliers and try to achieve
NLP convergence with the modified objective function.
yB ∈ {0, 1}|B| , yC ∈ R|C| , yF ∈ R|F|
If the crude-oil scheduling MILP is infeasible, the algorithm
The decomposability property is preserved: should terminate as the global MINLP is then proved infeasible. This
is an important difference of Lagrangian decomposition over typ-
v(P̃R ()) = v(PP ()) + v(P̃S ()) ical master/slave decompositions for which dedicated procedures
must be employed to handle infeasibilities in the slave problem
Finally, the modified dual problem (P̃D ) can be defined as: (Kelly & Zyngier, 2008).
The NLP solver may also fail to find a feasible solution to the
(P̃D ) min v(P̃R ()) heuristic problem. Again, globalization techniques can be used.

Additionally, one may add no-good or integer cuts to subsequent
This modified dual problem still provides a valid upper bound for scheduling and heuristic problems in order to prevent exploration
the original full-space problem (P). The following modified heuris- of identical discrete solutions yBK .
tic problem (P̃H (yBK )) is also defined. It is obtained from the original
heuristic problem (PH (yBK )) by dropping the nonlinear scheduling 7. Numerical illustration
constraints. It is used as a first heuristic step to get a good initial
point before solving the original heuristic NLP problem (PH (yBK )). In this section, several computational results are presented for
⎧ three approaches: the direct MINLP approach, a basic sequen-
⎪ max VPT xS − VCT yF tial scheduling-planning procedure and the proposed Lagrangian
⎪
⎪ s.t. fP (xF , xI , xS ) ≤ 0
⎪
⎪ decomposition method. Experiments have been performed on an
⎪
⎨ gP (xI ) ≤ 0 Intel Xeon 1.86 GHz processor using GAMS as the modeling and
(P̃H (yBK )) fS (yBK , yC , yF ) ≤ 0 algorithmic language. A 1000 s time limit has been used for each
⎪
⎪ yF − xF = 0
⎪
⎪ run. The following local NLP solvers have been used: CONOPT,
⎪
⎪ xF ∈ R|F| , xI ∈ R|I| , xS ∈ R|S| SNOPT and IPOPT. In our experiments, the convergence tolerance ε
⎩
yC ∈ R|C| , yF ∈ R|F| is set to 0.0001, the maximum step size parameter ˛ ¯ is set to 1 and
the step bound parameter ı̄ to 0.05.
2762 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

Table 9
Comparative performance of several MINLP algorithms for larger refinery problem (NLP solver: CONOPT).

MINLP solver 6 priority-slots 7 priority-slots

Objective value CPU time Objective value CPU time

Proposed 250.989 1045 s 250.128 4451 s

Sequential 116.814 46 s – 79 s
DICOPT – +3600 s – +3600 s
AlphaECP – +3600 s – +3600 s
SBB – +3600 s – +3600 s

strictly lower than 1, which corresponds to cases where the pure

subgradient multiplier update would violate some cutting planes
generated at previous iterations. The proposed approach automati-
cally overcomes this issue. The increase of CPU time between 6 and
7 priority-slots is mostly explained by the increase in size of the
MILP scheduling model. Figs. 11 and 12 display the evolution of the
objective value of various problems solved during the Lagrangian
iterations.
The optimal value of the Lagrange multipliers is ∗(X,A) = ∗(X,B) =
(Y,A) = ∗(Y,B) = 7. Figs. 13 and 14 display the evolution of the
∗

Lagrange multipliers during the Lagrangian iterations. The pro-

posed approach demonstrates its efficiency through stable updates
and fast convergence of the Lagrange multipliers.
A basic sequential procedure is introduced to compare with the
Lagrangian decomposition approach. First, the modified crude-oil
Fig. 18. Lagrangian iteration objective values for larger refinery problem (6 priority-
slots, NLP = CONOPT).
scheduling problem (P̃S ) is solved and the binary variables are fixed
to their solution value. Then, the modified and original heuris-
tic problems (P̃H (yB0 )) and (PH (yB0 )) are successively solved. This
Table 10
Blend compositions in the optimal solution of larger refinery problem.
procedure is not computationally expensive as it requires solv-
ing only one MILP, solved with CPLEX, and two NLPs, both solved
Crude Blend 1 Blend 2 Blend 3 with BARON. However, the solution obtained, if any, might not be
A 0.164 49.480 optimal.
B 33.293 16.707 Additionally, the Lagrangian decomposition approach is com-
C 50.000
pared to the direct approach, which consists of solving the
D 0.066 19.792
E 3.500 55.943 2.982 full-space problem (P) with various MINLP solvers. Table 4 shows
F 50.000 the computational performance of these MINLP algorithms. The
G 9.480 sequential approach quickly provides a feasible solution which is
8.0% lower than the global optimum (545.000 against 592.368). The
solvers DICOPT, AlphaECP and SBB cannot guarantee global opti-
The number of priority-slots for the crude-oil scheduling model mality of the solution returned while LINDOGLOBAL and BARON
is set to 6 and 7 (see Mouret et al., 2011). Tables 2 and 3 show are actual global MINLP solvers. DICOPT is able to find the global
iteration statistics for the Lagrangian decomposition method using optimal solution in many cases in reasonable CPU times. AlphaECP
SNOPT as the heuristic NLP solver. In particular, the optimal value never finds the global optimum. SBB seems to return the best solu-
of each problem solved is given as well as the optimal step size tions when CONOPT is used as NLP solver but it is requires large
calculated by the hybrid restrict dual problem and the cumulative CPU times. Neither LINDOGLOBAL nor BARON have found the global
CPU time at the end of each iteration (e.g., for 6 priority-slots, the optimum in the specified time limit.
first iteration took 3 s). Dashes are used when the information is not In comparison to these solvers, the proposed Lagrangian decom-
available (Hybrid Dual and Step Size for the first iteration), when the position approach proves to be very effective for the following
problem is unbounded (Pure Dual during the first few iterations), or reasons:
when it is locally infeasible (Original Heuristic in some iterations).
In each case, the global optimal solution is found (see under- • it is computationally effective (although DICOPT is faster);
lined entries in the Original Heuristic column) and proved optimal. • it always returns the global optimum;
It can be noted that in some iterations the step size variable ˛ is • it is very robust with the choice of NLP solver (although IPOPT is
significantly slower).

Fig. 15 depicts the composition of each crude blend in four dif-

ferent solutions. Crude blend X is mostly composed of crude A while
crude blend Y is mostly composed of crude B. If all solutions except
the second one (with objective value 568.077) are considered, one
may conclude that blend Y is “more profitable” than blend X because
the objective value increases when processing larger amounts of
blend Y and smaller amounts of blend X. Additionally, one could
say that increasing the amount of distilled crude increases the over-
all profit (solution 1 processes 253.034 units of crude, solution 3,
Fig. A.19. Layered artificial neural network 224.49 units of crude, and solution 4, 212.867 units of crude). Inter-
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2763

estingly, the second solution does not follow these observations. In pure restricted dual problem bounded. In order to achieve conver-
this solution, the amount of blend X is larger the amount of blend Y. gence after a few iterations, the Lagrangian gap is calculated using
Besides, the total amount of crude processed is the largest (276.923 the objective value of the hybrid dual problem, which is always
units of crude). This shows how difficult it is to approximate the bounded, as follows:
economic behavior of refinery operations, even for such a small case
study. Therefore, it is possible that linear approximations of refin- v(P̃R (K )) − v(P̂DK ) ?
≤ε
ery operations as proposed by Robertson et al. (2010) might not v(P̃R (K ))
be able to correctly evaluate the economic value of some feasible
The proposed approach converges within the Lagrangian gap in
solutions.
15 iterations. The final dual gap, although it does not represent a
valid optimality gap, is 3.8%. 83% of the time is spent on solving
the crude-oil scheduling MILPs. The optimal value of the Lagrange
8. Larger refinery problem
multiplier is given in Table 8. With 7 priority-slots, the decomposi-
tion procedure converges in 18 iterations and 4451 s, the increase
In this section, the proposed Lagrangian decomposition
of CPU time being explained mostly by the increased number of
approach is applied on a larger refinery example and compared
priority-slots. The solution obtained is slightly lower than the pre-
to standard MINLP solvers and the basic sequential procedure
vious one because the algorithm did not fully converge within tight
presented in Section 7. Instead of using classical fixed-yield (see
tolerances.
Section 2) or swing cut formulations (for example, see Zhang et al.,
Table 9 presents computational performances of the different
2001), the refinery planning problem is based on a crude distil-
approaches. Clearly, the Lagrangian decomposition procedure is
lation simulation model developed by Gueddar and Dua (2011).
much more robust than the other algorithms. Only the sequential
This crude distillation model is based on a layered artificial neu-
approach was able to deliver a solution with 6 priority-slots, but
ral network (ANN). This ANN is generated by solving an MINLP
its objective value is much lower (50% less profit) than the best
that aims at fitting empirical data for atmospheric distillation of
known solution. Indeed, the scheduling solution obtain during the
several crudes (crude assays) while simplifying the calculations.
first stage of the sequential procedure does not take into account
The model obtained is able to predict cut yields and cut properties
the economic impact on the refinery planning problem, which leads
from the chosen cut points and properties of the inlet crude. Fig. 16
to poor decisions. The standard local MINLP solvers were not able
depicts the full nonlinear planning model. Three different types of
to find a solution within 1 h because the MINLP model is too large:
crude blends are processed in three different CDU operating modes
3645 variables (84 binary) and 4177 constraints.
and five crude cuts are produced: LPG (liquefied petroleum gas),
Table 10 shows the crude compositions of the optimal solution
naphta, kerosene, diesel, and residue. Each discrete CDU mode is
with objective value 250.989. The two CDUs are mostly operated in
defined by the distillation cut point between diesel and residue
mode 2 which corresponds to the average cut point for diesel and
fractions: 340, 360, or 380 ◦ C. The decision variables for each CDU
residue cuts.
mode are composed of individual crude flows and distillation cut
points between naphta, kerosene and diesel fractions. The three
streams produced for each fraction are then blended into cut pools. 9. Conclusion
The bilinear pooling constraints discussed in Section 2 are classi-
cally used to calculate the properties of each cut pool. Table 5 shows In this work, a novel approach towards the integration of plan-
market prices for each distillation cut and corresponding prop- ning and scheduling has been developed in the context of oil
erty specifications. Crude availability constraints are also included refining. In particular, a precise crude-oil operations scheduling
in the model in accordance with the description of the crude-oil model and a coarse refinery planning model were optimized simul-
scheduling problem. This refinery planning model is composed of taneously using a new Lagrangian decomposition scheme, which
1831 variables and 1817 constraints. The full mathematical model makes use of Lagrange multipliers as a way to communicate eco-
is described in Appendix A. nomic information between the two subsystems. The methodology
The crude-oil scheduling problem is based on example 3 from leads to a classical primal-dual iterative algorithm to solve the
Lee et al. (1996) (see Fig. 17) with the parameters described in Lagrangian dual problem. The critical multiplier update step is
Table 6. It consists of three crude arrivals, three storage tanks, three performed by solving a new hybrid restricted dual problem. This
charging tanks (one for each type of crude mixture), and two identi- approach combines the strengths of cutting planes and subgradient
cal CDUs whose respective scheduled feedstocks are merged when steps and does not require to define heuristic updates of parameters
linked to the refinery planning problem. Seven different crudes are during iterations.
available, and 14 transfer operations can be executed to prepare the Although it is not guaranteed, our results achieved a 0% dual
different crude blends. When 6 priority-slots are used, the crude- gap for the smaller case study. It is well-known that augmented
oil scheduling problem is composed of 1814 variables (84 binary) Lagrangian techniques (see Li & Ierapetritou, 2010) can ensure to
and 2338 constraints. close the dual gap for any instance. However, this would require to
Table 7 and Fig. 18 show the iteration statistics for the solve the refinery planning subproblem to global optimality, which
Lagrangian decomposition method using 6 priority-slots for the is not yet achievable in an industrial context. It is therefore more
crude-oil scheduling model. The maximum step size parameter ˛ ¯ practical to use standard duality in order to obtain feasible heuristic
is set to 1 and the step bound parameter ı̄ to 0.05. Because we were solutions for the integrated problem.
not able to solve the refinery planning problem to global optimal- The proposed Lagrangian decomposition procedure has been
ity, we used CONOPT as the NLP solver. Other local NLP solvers applied to a larger and more complex refinery problem. The refinery
as SNOPT and IPOPT did not perform as well (slow convergence, planning problem is based on a crude distillation simulation model
poor solutions or local infeasibilities). As the refinery planning that has been developed by Gueddar and Dua (2011) using an arti-
problem is not solved to global optimality, it is not possible to rig- ficial neural network to fit empirical data. The crude-oil scheduling
orously estimate the optimality of the solution obtained for the problem is the third example of Lee et al. (1996). The results have
full-space problem. Therefore, the convergence tolerance ε is set to shown that the proposed approach remains very competitive com-
0.01 instead of 0.0001 as in the case study. Furthermore, many iter- pared to other MINLP algorithms, namely the two-step MILP-NLP
ations are needed to generate enough cutting planes and make the sequential procedure and standard MINLP solvers.
2764 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

Acknowledgment complexity. This tuning step is performed by solving a training

MINLP. We denote CDUANN(x, u) the set of constraints defining the
The authors are grateful to Total Raffinage & Marketing for finan- relation between the ANN inputs x and outputs u. In particular, the
ip
cial support of this project. model inputs include crude properties qC and CDU cut points  k
(k ∈ {naphta, kerosene, diesel}), and the model outputs include cut
Appendix A. Mathematical model for the refinery planning ijkp
yields ˛ijk and crude cut properties q1 . All crude property inputs
problem are fixed while CDU cut points are variable. The cut point between
diesel and residue cuts can take three discrete values defining the
In this section, the ANN model developed in Gueddar and Dua three CDU operating modes. All the outputs are variable. The full
(2011) for CDU simulations is presented. It is based on a layered refinery planning model is expressed as follows:

Appendix B. Mathematical model for the crude-oil

directed graph which represents the model calculations (see Fig. operations scheduling problem
A.19). Each node in the input/output layers correspond to one
input/output variable. Each node j = 1, . . ., Nn in the intermediate In this section, the mathematical model for crude-oil operations
layer l = 1, . . ., Nl corresponds to an activation variable alj and a trans- scheduling problems corresponding to the MOS time representa-
formed variable hlj . The activation variables are calculated from the tion (see Mouret et al., 2011) is presented.
transformed variables of the previous layer using an affine expres- It makes use of the following sets and parameters.
sion while the transformed variables are calculated by applying the
hyperbolic tangent to their associated activation variable. The ANN • T = {1, . . ., n} is a totally ordered set of priority-slots (indices i, j,
equations are expressed as follows: i1 , i2 )
• W is the set of all operations (indices v, v , v1 , v2 )

Nx
• Pv1 v2 = 1 denotes a precedence constraint between operations v1
a1j = wji1 xi + b1j , j = 1, . . . , Nn (A.1)
and v2
i=1
• PW W = 1 denotes a precedence constraint between set of oper-
1 2
hlj = tanh alj , l = 1, . . . , Nh , j = 1, . . . , Nn (A.2) ations W1 and W2
• WU ⊂ W is the set of unloading operations (WU = {1, 2, 3} for

Nn
COSP2)
alj = wjil hl−1
i
+ blj , l = 2, . . . , Nh , j = 1, . . . , Nn (A.3) • WT ⊂ W is the set of tank-to-tank transfer operations (WT = {4, . . .,
i=1 10} for COSP2)
• WD ⊂ W is the set of distillation operations (WD = {11, . . ., 14} for

Nn
COSP2)
uk = Wki hN
i
h +B ,
k k = 1, . . . , No (A.4)
• R is the set of resources (i.e. tanks, units): R = RV ∪ RS ∪ RC ∪ RD
i=1
• RV ⊂ R is the set of vessels
The model uses the following parameters: • RS ⊂ R is the set of storage tanks
• RC ⊂ R is the set of charging tanks
• Nx is the number of inputs • RD ⊂ R is the set of distillation units
• No is the number of outputs • Ir ⊂ W is the set of inlet transfer operations on resource r
• Nh is the number of intermediate layers • Or ⊂ W is the set of outlet transfer operations on resource r
• Nn is the number of nodes in each intermediate layer • C is the set of products (i.e. crudes)
• wjil , blj , Wki , Bk are parameters specific to the ANN • K is the set of product properties (e.g. crude sulfur concentration)
• H is the scheduling horizon
Dua (2010) demonstrates how to tune the ANN parameters • S v ∈ [0, H] is a lower bound on the start time of unloading opera-
by minimizing the total prediction error as well as the ANN tions v ∈ WU
 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766 2765

• NOv1 v2 is 1 if operations v1 and v2 must not overlap, 0 if they are Ziv = 1 if operation v is assigned to priority-slot i, Ziv = 0 other-
allowed to overlap wise
• GNO is the non-overlapping graph (see Mouret et al., 2011) whose • Time variables Siv ≥ 0, Div ≥ 0, Eiv ≥ 0 i ∈ T, v ∈ W
adjacency matrix is NO Siv is the start time of operation v if it is assigned to priority-slot
• [Vvt , Vvt ] are bounds on the total volume transferred during trans- i, Siv = 0 otherwise
fer operation v; in all instances, Vvt = 0 for all operations except Div is the duration of operation v if it is assigned to priority-slot
i, Div = 0 otherwise
unloadings for which Vvt = Vvt is the volume of crude in the marine
Eiv is the end time of operation v if it is assigned to priority-slot
vessel
i, Eiv = 0 otherwise
• [N D , ND ] are the bounds on the number of distillations • Operation variables V t ≥ 0 and Vivc ≥ 0 i ∈ T, v ∈ W, c ∈ C
• [FRv , FRv ] are flowrate limitations for transfer operation v iv
Vitv is the total volume of crude transferred during operation v
• [xvk , xvk ] are the limits of property k of the blended products
if it is assigned to priority-slot i, Vitv = 0 otherwise
transferred during operation v
Vivc is the volume of crude c transferred during operation v if it
• xck is the value of the property k of crude c
is assigned to priority-slot i, Vivc = 0 otherwise
• [Ltr , Lrt ] are the capacity limits of tank r • Resource variables Lt and Lirc i ∈ T, r ∈ R, c ∈ C
ir
• [Dr , Dr ] are the bounds of the demand on products to be trans- t is the total accumulated level of crude in tank r ∈ R ∪ R
Lir S C
ferred out of the charging tank r during the scheduling horizon before the operation assigned to priority-slot i
Lirc is the accumulated level of crude c in tank r ∈ RS ∪ RC before
The following variables are introduced in the model. the operation assigned to priority-slot i

• Assignment variables Ziv ∈ {0, 1} i ∈ T, v ∈ W The crude-oil scheduling MOS model is expressed as follows:
2766 S. Mouret et al. / Computers and Chemical Engineering 35 (2011) 2750–2766

References Kelley, J. E. (1960). The cutting-plane method for solving convex programs. Journal
of the Society for Industrial and Applied Mathematics, 8(4), 703–712.
Adhya, N., Tawarmalani, M., & Sahinidis, N. V. (1999). A Lagrangian approach Kelly, J. D., & Mann, J. L. (2003). Crude oil blend scheduling optimization: An appli-
to the pooling problem. Industrial and Engineering Chemistry Research, 38(5), cation with multimillion dollar benefits. Hydrocarbon Processing, 82(6), 47–53.
1956–1972. Kelly, J. D., & Zyngier, D. (2008). Hierarchical decomposition heuristic for schedul-
Alhajri, I., Elkamel, A., Albahri, T., & Douglas, P. L. (2008). A nonlinear programming ing: Coordinated reasoning for decentralized and distributed decision-making
model for refinery planning and optimisation with rigorous process models and problems. Computers and Chemical Engineering, 32(11), 2684–2705.
product quality specifications. International Journal of Oil, Gas and Coal Technol- Lee, H., Pinto, J. M., Grossmann, I. E., & Park, S. (1996). Mixed-integer linear pro-
ogy, 1(3), 283–307. gramming model for refinery short-term scheduling of crude oil unloading with
Audet, C., Brimberg, J., Hansen, P., Le Digabel, S., & Mladenović, N. (2004). Pooling inventory management. Industrial and Engineering Chemistry Research, 35(5),
problem: Alternate formulations and solution methods. Management Science, 1630–1641.
50(6), 761–776. Lemaréchal, C. (1974). An algorithm for minimizing convex functions. In Proceedings
Barahona, F., & Anbil, R. (2000). The volume algorithm: Producing primal solutions IFIP’74 congress (pp. 552–556).
with a subgradient method. Mathematical Programming, 87(3), 385–399. Li, W., Hui, C., & Li, A. (2005). Integrating CDU, FCC and product blending models
Bazaraa, M. S., & Sherali, H. D. (1981). On the choice of step size in subgradient into refinery planning. Computers and Chemical Engineering, 29(9), 2010–2028.
optimization. European Journal of Operational Research, 7(4), 380–388. Li, Z., & Ierapetritou, M. G. (2010). Production planning and scheduling integration
Bodington, C. E., & Baker, T. E. (1990). A history of mathematical programming in through augmented Lagrangian optimization. Computers and Chemical Engineer-
the petroleum industry. Interfaces, 20(4), 117–127. ing, 34(6), 996–1006.
Camerini, P. M., Fratta, L., & Maffioli, F. (1975). On improving relaxation methods by Maravelias, C. T., & Sung, C. (2009). Integration of production planning and
modified gradient techniques. Mathematical Programming Studies, 3, 26–34. scheduling: Overview, challenges and opportunities. Computers and Chemical
Cheney, E. W., & Goldstein, A. A. (1959). Newton’s method for convex programming Engineering, 33(12), 1919–1930.
and Tchebycheff approximation. Numerische Mathematik, 1(1), 253–268. Marsten, R. E., Hogan, W. W., & Blankenship, J. W. (1975). The boxstep method for
Dua, V. (2010). A mixed-integer programming approach for optimal configuration large-scale optimization. Operations Research, 23(3), 389–405.
of artificial neural networks. Chemical Engineering Research and Design, 88(1), Misener, R., & Floudas, C. A. (2009). Advances for the pooling problem: Model-
55–60. ing, global optimization, and computational studies. Applied and Computational
Erdirik-Dogan, M. E., & Grossmann, I. E. (2008). Simultaneous planning and Mathematics, 8(1), 3–22.
scheduling of single-stage multi-product continuous plants with parallel lines. Moro, L. F. L., & Pinto, J. M. (2004). Mixed-integer programming approach for short-
Computers and Chemical Engineering, 32(11), 2664–2683. term crude oil scheduling. Industrial and Engineering Chemistry Research, 43,
Fisher, M. L. (1981). The Lagrangian relaxation method for solving integer program- 85–94.
ming problems. Management Science, 27(1), 1–18. Mouret, S., Grossmann, I. E., & Pestiaux, P. (2009). A novel priority-slot based
Fisher, M. L. (1985). An application oriented guide to Lagrangian relaxation. Inter- continuous-time formulation for crude-oil scheduling problems. Industrial and
faces, 15, 10–21. Engineering Chemistry Research, 48(18), 8515–8528.
Floudas, C. A., & Visweswaran, V. (1993). A primal-relaxed dual global optimization Mouret, S., Grossmann, I. E., & Pestiaux, P. (2011). Time representations and math-
approach. Journal of Optimization Theory and Applications, 78(2), 187–225. ematical models for process scheduling problems. Computers and Chemical
Foulds, L. R., Haugland, D., & Jörnsten, K. (1992). A bilinear approach to the pooling Engineering, doi:10.1016/j.compchemeng.2010.07.007
problem. Optimization, 24(1–2), 165–180. Neiro, S. M. S., & Pinto, J. M. (2006). Lagrangean decomposition applied to mul-
Frangioni, A. (2005). About Lagrangian methods in integer optimization. Annals of tiperiod planning of petroleum refineries under uncertainty. Latin American
Operations Research, 139(1), 163–193. Applied Research, 36(4), 213–220.
Goffin, J. L., Haurie, A., & Vial, J. P. (1992). Decomposition and nondifferentiable Pinto, J. M., & Moro, L. F. L. (2000). A planning model for petroleum refineries.
optimization with the projective algorithm. Management Science, 38(2), 284– Brazilian Journal of Chemical Engineering, 17(4–7), 575–586.
302. Quesada, I., & Grossmann, I. E. (1995). A global optimization algorithm for linear
Gueddar, T., & Dua, V. (2011). Novel model reduction techniques for refinery-wide fractional and bilinear programs. Journal of Global Optimization, 6(1), 39–76.
energy optimization. Applied Energy Journal, in press. Robertson, G., Palazoglu, A., & Romagnoli, J. A. (2010). Refinery scheduling of crude
Guignard, M. (2003). Lagrangean relaxation. Top, 11(2), 151–228. oil unloading, storing, and processing considering production level cost. In 20th
Held, M., & Karp, R. M. (1971). The traveling-salesman problem and minimum span- European symposium on computer aided process engineering – ESCAPE20. Vol. 28
ning trees. Part II. Mathematical Programming, 1(1), 6–25. of Computer Aided Chemical Engineering (pp. 1159–1164).
Held, M., & Karp, R. M. (1974). Validation of subgradient optimization. Mathematical Zhang, J., Kim, N., & Lasdon, L. (1985). An improved successive linear programming
Programming, 6(1), 62–88. algorithm. Management Science, 31(10), 1312–1331.
Jia, Z., Ierapetritou, M. G., & Kelly, J. D. (2003). Refinery short-term scheduling using Zhang, J., Zhu, X. X., & Towler, G. P. (2001). A level-by-level debottlenecking
continuous time formulation: Crude-oil operations. Industrial and Engineering approach in refinery operation. Industrial and Engineering Chemistry Research,
Chemistry Research, 42(13), 3085–3097. 40(6), 1528–1540.