Formulations For A Problem of Petroleum Transportation
Formulations For A Problem of Petroleum Transportation
Formulations For A Problem of Petroleum Transportation
Discrete Optimization
a r t i c l e
i n f o
Article history:
Received 11 March 2013
Accepted 19 January 2014
Available online 27 January 2014
Keywords:
Mixed Integer Programming
Transportation problem
Column generation-based heuristic
a b s t r a c t
Oil tankers play a fundamental role in every offshore petroleum supply chain and due to its high price, it
is essential to optimize its use. Since this optimization requires handling detailed operational aspects,
complete optimization models are typically intractable. Thus, a usual approach is to solve a tactical level
model prior to optimize the operational details. In this case, it is desirable that tactical models are as precise as possible to avoid too severe adjustments in the next optimization level. In this paper, we study
tactical models for a crude oil transportation problem by tankers. We did our work on the top of a previous paper found in the literature. The previous model considers inventory capacities and discrete lot
sizes to be transported, aiming to meet given demands over a nite time horizon. We compare several
formulations for this model using 50 instances from the literature and proposing 25 new harder ones.
A column generation-based heuristic is also proposed to nd good feasible solutions with less computational burden than the heuristics of the commercial solver used.
2014 Elsevier B.V. All rights reserved.
1. Introduction
Oil tankers play a fundamental role in every offshore petroleum
supply chain. Its optimization is usually divided in three levels:
strategic, tactical and operational. Strategic decisions deal with
eet sizing, facility location and capacity sizing. Tactical decisions
deal with production and distribution planning, transportation
mode selection, storage allocation and order picking strategies.
Finally, operational decisions deal with shipment and vehicle
dispatching (Ghiani, Laporte, & Musmanno, 2004). When dealing
with a planning level, all decisions made in the previous levels
are assumed to be known. Since each level considers more details
than the previous ones, past decisions usually have to be adjusted.
Hence, it is important to make both the strategic and tactical models as detailed as possible to guarantee low adjustments in the
operational level.
In this paper, we study a tactical optimization model for crude
oil distribution by tankers. The problem consists of scheduling the
shipments through routes linking platforms (offshore production
sites) and terminals (onshore consumer sites). The objective is to
ship the products from the platforms to supply the terminals with
minimum transportation cost for a nite planning horizon. For
each site, the inventory levels must lie between a lower and an
upper bound to avoid the lack or excess of product. Also, for each
site, the demand or the production is given for the whole planning
Corresponding author. Tel.: +55 2126295708.
E-mail addresses: luizaizemberg@gmail.com (L. Aizemberg), hugoharry
@gmail.com (H.H. Kramer), artur@producao.uff.br (A.A. Pessoa), uchoa@prod
ucao.uff.br (E. Uchoa).
0377-2217/$ - see front matter 2014 Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.ejor.2014.01.036
83
p: Platform, p 2 P;
t: Terminal, t 2 T;
c: Class of tanker, c 2 C;
d: day, d 2 f1; . . . ; Dg.
Subsets:
Tp # T: Terminals allowed to send tankers to a platform
p 2 P;
Pt # P: Platforms allowed to receive tankers from a terminal t 2 T;
Cp # C: Classes of tankers allowed to pick up oil from platform p 2 P.
Input data:
NF min
subject to
2. Mathematical formulations
X X
c;d
V c zp;t
8p 2 P; d 2 f1;...; Dg;
t2Tpc2Cp
This section presents the mathematical formulations considered for the problem in study.
The optimization problem dened in Rocha et al. (2011) is the
following. Let P be the set of platforms and T the set of terminals.
The planning horizon consists of D days. A unique product should
be shipped from platforms to terminals to satisfy given demands.
The inventory levels of platforms and terminals are limited by a
lower and an upper bound, i.e., the lack or surplus of the product
is not allowed at any point of the network. Oil tankers are grouped
p2P t2Tpc2Cp d1
X X
c;dDp;t
V c zp;t
8t 2 T; d 2 f1;...;Dg; 3
p2Ptc2Cp
0 6 ep;d 6 CAPp;d
0 6 et;d 6 CAP t;d
zc;d
p;t 2 f0; 1g
8p 2 P; d 2 f1;...;Dg;
8t 2 T; d 2 f1;.. .;Dg;
5
6
The objective function (1) aims to minimize the total transportation costs. Constraints (2) calculate the inventory balance at each
platform, while constraints (3) calculate the inventory balance at
each terminal. Constraints (4) and (5) assure that inventory levels
84
at each platform and terminal at each day will not be less than
their minimum storage capacity nor greater than their maximum
storage capacity.
Some of the problem assumptions deserve discussion:
no more than one tanker of each class can be used in each route
and day: this assumption is not realistic, we just followed the
previous paper. But it is easy to drop this assumption, just setting the binary variables as integer;
unlimited number of tankers in each class: we believe this
assumption is realistic because it is based on a real problem
where the oil company rents the tankers;
the instances use no realistic transportation and (un)loading
times: we could use realistic values, but this will not change
the formulation performance.
The NF is the starting point for a second formulation, which
needs some additional denitions listed on the following. Dene
Lp as the greatest common divisor (GCD) of all the shipping capacities of all classes of tankers that are allowed to pick up oil from
platform p. Let kmin p; d be minimum number of lots of size Lp that
should be shipped before day d from platform p to avoid its maximum inventory capacity to be exceeded. This number is given
by the formula:
KC min
D
XX XX
2F c Dp;t zc;d
p;t ;
p2P t2Tpc2Cp d1
subject to
d X
XX
d X
XX
V c zc;p;ts 6 AC p;d
8p 2 P; d 2 f1;...;Dg;
10
V c zc;p;ts P DAp;d
8p 2 P; d 2 f1;...;Dg;
11
t2Tp s1 c2Cp
d X
XX
d
X
ep;0
Pp;s :
c;sDp;t
V c zp;t
P AC t;d
8t 2 T; d 2 f1;...;Dg; 12
p2Pt s1 c2Cp
s1
AC p;d
;
kmax p; d
Lp
With the denitions given above, constraints (4) are replaced
with:
8p 2 P;
7
AC t;d
;
Lt
AC t;d et;0
t2Tp s1 c2Cp
kmin t; d
8t 2 T; d 2 f1;. ..;Dg:
AC p;d CAPp;d
;
kmin p; d
Lp
AC p;d
d
X
C t;s :
s1
CAPt;d AC t;d
;
kmax t; d
Lt
Given this, constraints (5) are replaced with:
d X
XX
c;sDp;t
V c zp;t
6 DAt;d
8t 2 T; d 2 f1;...;Dg;
13
p2Pt s1 c2Cp
zc;d
p;t 2 f0;1g
8p 2 P; d 2 f1;... ;Dg
15
8p 2 P; d 2 f1;...;Dg
16
t2Tp s1 c2Cp
d X
XX
t2Tp s1 c2Cp
c;sDp;t
V c zp;t
6 kmax t;dLt
8t 2 T; d 2 f1;...;Dg
17
P kmin t;dLt
8t 2 T; d 2 f1;...;Dg:
18
p2Pt s1 c2Cp
d X
XX
c;sDp;t
V c zp;t
p2Pt s1 c2Cp
XX X
2F c Dp;t xc;D
p;t ;
p2P t2Tpc2Cp
85
19
subject to
c;d
c;d1
xp;t
xp;t
zc;d
8p 2 P; t 2 Tp;c 2 Cp; d 2 f1;...;Dg;
p;t
X X
c;d
V c xp;t
6 kmax p;dLp 8p 2 P; d 2 f1;...;Dg;
20
21
t2Tpc2Cp
X X
c;d
V c xp;t
P kmin p;dLp
8p 2 P; d 2 f1;...;Dg;
22
t2Tpc2Cp
X X
c;dDp;t
6 kmax t;dLt
8t 2 T; d 2 f1;...;Dg;
23
c;dDp;t
P kmin t;dLt
8t 2 T; d 2 f1;...;Dg;
24
V c xp;t
p2Ptc2Cp
X X
V c xp;t
p2Ptc2Cp
c;d
2 Z
xp;t
c;d
zp;t
2 f0;1g or 0;1
25
26
The objective function (19) is to minimize the total transportation cost. Constraints (20) dene the x variables in terms of the z
variables. Constraints (21) and (22) ensure that the inventory level
in each platform p at day d is not less than the minimum inventory
capacity and does not exceed the maximum inventory capacity.
Constraints (23) and (24) avoid the inventory level in each terminal
t at day d to exceed the maximum inventory capacity and to be less
than zero. At last, constraints (25) dene the x variables as nonnegative integers and constraints (26) dene the z variables as binaries
or continuous variables in the range zero and one. Note that constraint (20) ensures the integrality of variables z regardless of their
denitions.
Proposition 2.1.
relaxation.
ep;d AC p;d
d X
XX
V c zc;p;ts
8p 2 P;
d 2 f1; . . . ; Dg
t2Tp s1 c2Cp
XX
C p;s kp;s ;
27
p2P s2Sp
subject to
X
kp;s 1 8p 2 P;
28
s2Sp
AC t;d
Lt 8t 2 T; d 2 f1;...;Dg;
L
t
p2P s2Sp s1
d
XXX
CAP t AC t;d
Lt 8t 2 T; d 2 f1;...;Dg;
V t;p;ssDp;t kp;s 6
L
t
p2P s2S s1
d
XXX
V t;p;ssDp;t kp;s P
29
30
kp;s 2 f0;1g
8p 2 P; s 2 Sp :
31
The objective function (27) aims to minimize the sum of transportation costs. Constraints (28) ensure that exactly one solution
s 2 Sp variable is chosen for each platform p 2 P. Constraints (29)
avoid the inventory level at each terminal to be less than zero. Constraints (30) ensure that the maximum inventory capacity of each
terminal will not be exceeded. Finally, constraints (31) dene the k
variables as binaries.
From (7),
d X
XX
V c zc;p;ts
t2Tp s1 c2Cp
8p 2 P; d 2 f1; .. .; Dg
which is equivalent to (15) and (16). The same procedure can be
used to prove that (3) and (8) have the same linear relaxation solution set of (17) and (18) over the z variables. Thus, RKC has the
same linear relaxation solution set over the z variables as BF .
Now, consider the RKC. Dene the x variables as follows:
xc;d
p;t
d
X
zc;p;ts
s1
86
SP p min
D
X XX
c;d
2F c Dp;t V c ht;dDp;t rt;dDp;t zp;t
pp ; 32
t2Tpc2Cp d1
subject to
d X
XX
33
t2Tp s1 c2Cp
d X
XX
zc;p;ts V c P DAp;d
8d 2 f1; . . . ;Dg;
34
t2Tp s1 c2Cp
zc;d
p;t 2 f0; 1g
35
The purpose of the Pricing Subproblem is to generate a k variable to be added to the RMP with the lowest reduced cost, which
is represented by the objective function (32). Constraints (33)
avoid the inventory level at each platform to be less than zero. Constraints (34) ensure that the maximum inventory capacity at each
platform will not be exceeded. At last, z variables are dened as
binary by (35). Even though the Pricing Subproblem is presented
as an Integer Linear Program, it can be solved more efciently by
dynamic programming.
3. Column generation-based heuristic
This section presents the column generation-based heuristic
(CGH). First, we show the dynamic programming algorithm to
solve the subproblems faster than the MIP of subSection 2.4.1.
Next, we show the heuristic procedure.
3.1. Pricing algorithm
In the following denitions we assume that the elements of the
set C are arbitrarily numbered from 1 to jCj. We dene our dynamic
programming cost table for platform p with Costp d; c; e being the
smallest partial reduced cost from day d to the end of the time
horizon, considering only the tanker classes c; . . . ; jCj at day d
(and all classes for the next days), and inventory level e at the
beginning of day d. Let Costp d; c; e 0 for d D 1. The dynamic
programming recursion is the following:
(
!
)
8
X
X c;d
>
>
>
min
Cost
d;c
1;e
P
V
a
C
a
p
c t
p;d
>
p;t t ; c 1
>
a2P 1;p;c;d;e
>
>
t2T
t2T
>
>
(
!
)
>
<
X
X c;d
min Cost p d;c 1;e
Cost p d;c; e
V c at
C p;t at ; 1 < c < jCj
a2P 2;p;c;d;e
>
>
t2T
t2T
>
>
(
!
)
>
>
>
X
X c;d
>
>
>
V c at
C p;t at ; c jCj;
: min Cost p d 1; 1; e
a2P 3;p;c;d;e
t2T
t2T
where
X
V c at g;
capacity of the platforms. For the harder class of instances, its performance is quite superior to when a MIP is used. Fig. 1 shows an
example with two lot sizes C f1; 2g. The daily production is always integer. There are two terminals with route to all platforms.
Each lot size can be sent only once per day, for each terminal. In
this example, the inventory capacity of the platform is CAPp 4
units, the daily production is P p;d 2 units, d 2 f1; 2; 3g, which is
added between the rst and the second lot size, and the initial
inventory is 2 units. The rst lot has size V 1 1 unit and costs
1.5. The second lot has size V 2 2 units and costs 2. Both can be
sent to all terminals, with the same cost. If there are two or more
terminals to send a lot size with different costs, the algorithm
chooses always the cheapest one according to the reduced costs.
In Fig. 1, the costs for each inventory level, day and lot size are
represented by the circles. At d 1 and c 1, there is an initial
inventory of 2 units and the lowest cost is Costp 1; 1; 2 4, which
is also the optimal cost. From d 1 and c 1 to d 1 and c 2,
there are three possibilities: (i) to send one lot (at 1, for some t),
resulting in inventory of 3 units and cost 3 1:5; (ii) to send two lots
(at 1; 8t 2 T), resulting in inventory of 2 units and cost 2 3; or
(ii) do not send a lot at 0; 8t 2 T, resulting in inventory of 4
units and cost 4. The third option is chosen. From d 1 and
c 2 to d 2 and c 2, no lot is sent. Then, one lot is sent from
d 2 and c 2 to d 3 and c 1, and from d 3 and c 2,
resulting in an inventory of 4 units at the end of the time horizon.
The main advantage of the dynamic programming procedure is
that the economy of scale on the lot sizes does not impact on the
runtime, differently of the MIP formulation where we observe a
major impact. This happens because the magnitude of the economy of scale only affects the reduced cost values, not the number
of calculations made by the algorithm.
For the case where the daily production is fractional in at least
one day, we need to do a preprocessing step. We accumulate the
daily productions, and if the accumulated value is fractional, we
decrease the inventory capacity by one unit to have space for the
fractional part. The adjusted daily production is the integer part
of the accumulated production, not added yet to the inventory.
Table 1 shows an example with fractional production where the
inventory capacity has 6 units. The rst column shows the daily
original production, which is fractional and equal 2.75 for all days.
The second column shows the accumulated production. Note that
there is only one integer value. Third column shows the original
capacity, 6 units. Fourth column shows the adjusted capacity,
t2T
P 2;p;c;d;e
X
fa 2 f0; 1g j0 6 e
V c at g;
jTj
Table 1
Dealing with fractional production.
t2T
X
V c at 6 CAPp g;
t2T
and C c;d
p;t 2F c Dp;t V c ht;dDp;t rt;dDp;t . The optimal reduced cost
for platform p is given by Costp 1; 1; ep;0 .
The Dynamic Programming algorithm is pseudopolynomial. Its
computational complexity depends on the maximum inventory
Original
production
Accumulated
production
Original
capacity
Adjusted
capacity
Adjusted
production
2.75
2.75
2.75
2.75
2.75
2.75
5.5
8.25
11
13.75
6
6
6
6
6
5
5
5
6
5
2
3
3
3
2
87
4. Computational results
In this section we show the results obtained by the described
formulations using the instances proposed in Rocha et al. (2011)
and a new set of harder instances. All the tests, including those
for the formulation with the best performance proposed by Rocha
et al. (2011) (HullRel), were run in a PC with an Intel Core 2 Quad
Q8300 2.50 gigahertz CPU, 2 gigabytes RAM under Windows 7 OS
using ILOG CPLEX 12.1. All instances have 9 platforms and 2 terminals, and are divided into 4 classes of 25 instances each with the
following characteristics:
In this work, the results obtained for the rst class of instances
(easy) will not be reported, due to the fact that such instances are
easily solved by BF combined with a commercial solver. The results for each instance class are shown separately. The tables contain, for each formulation, the average results obtained. Among the
results, two performance measures are shown, given by:
gapLB 100
UBb LBf
;
UBb
36
gapUB 100
UBf UBb
;
UBb
37
where UBb is the best integer feasible solution found among all formulations, LBf is the lower bound for a given formulation, and UBf is
the best integer feasible solution found by a given formulation.
The results obtained for the medium and hard classes of instances are shown in Tables 2 and 3, respectively. In the tables,
the column Formulation shows the formulations tested. The formulation ARCbin is the ARC formulation with z variables as binaries, ARCcont is the ARC formulation with z variables as
continuous between 0 and 1 and AC is the ARC formulation without rounding the inventory. Column Time (seconds) shows the
average execution time in seconds for each class of instances. For
all tests performed, the maximum execution time allowed for each
instance is 720 seconds. Column gapLB rootlp (%) shows the average gapLB at the root node before CPLEX cuts, calculated using the
gapLB of each instance given by (36) considering LBf as the lower
bound at the root node before CPLEX cuts. Column gapLB root
(%) shows the average gapLB at the root node after CPLEX cuts, also
calculated using the gapLB of each instance given by (36) considering LBf as the lower bound at the root node after CPLEX cuts.
Column gapLB nal (%) contains the average gapLB at the end of
the execution time of each instance and it is also calculated by
(36), where LBf is the lower bound at the end of the execution time
of each instance. Column gapUB (%) shows the average gapUB obtained at the end of executions and calculated by (37) for each instance. Column No. of nodes contains the average number of
nodes solved at the branch-and-bound algorithm. Finally, column
Inst. solved shows the number of instances of each class solved
to optimality.
Concerning the gapLB root lp , the results include the CPLEX
preprocessing routines. If they were not used, NF ; KC and AC
(see Proposition 2.2) would provide the same relaxation at root
node of the branch-and-bound tree, and the same would happen
with BF ; RKC and ARC (see Proposition 2.1). Indeed, for the
medium and hard instance classes, the average values of the linear
programming relaxation of the rounded formulations are 0.42 and
3.04, respectively, and for the formulations without rounding are
88
Table 2
Average results and number of instances solved for the mathematical formulations and results taken from Rocha et al. (2011) (HullRel) for the medium class of instances.
Formulation
Time (seconds)
gapLB root lp
gapUB (%)
No. of nodes
Inst. solved
NF
BF
KC
RKC
AC
ARCbin
ARCcont
DWM
HullRel
720
685.47
405.50
405.50
30.65
31.76
30.61
351.53
1.41
0.36
0.35
0.35
0.34
0.33
0.33
0.33
1.19
0.30
0.23
0.23
0.19
0.20
0.18
0.83
0.24
0.12
0.12
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.04
2,595,996
3,944,016
569,080
785,543
17,288
12,319
14,663
0
2
11
11
24
24
24
13
Table 3
Average results and number of instances solved for the mathematical formulations and results taken from Rocha et al. (2011) (HullRel) for the hard class of instances.
Formulation
Time (seconds)
gapLB root lp
gapUB (%)
No. of nodes
Inst. solved
NF
BF
KC
RKC
AC
ARCbin
ARCcont
DWM
HullRel
720.00
720.00
687.51
664.28
103.93
111.13
102.94
691.40
3.14
2.98
2.84
2.84
2.38
2.38
2.38
1.24
2.74
2.62
1.76
1.76
1.32
1.30
1.30
1.98
1.95
1.11
1.11
0.77
0.77
0.77
1.07
1.07
0.04
0.04
0.04
1.24
1,294,071
1,326,605
157,950
156,079
55,632
52,975
57,962
0
0
2
2
22
22
22
Table 4
Size of all formulations for an instance of the harder class.
Number of variables
Number of constraints
Number of nonzero coefs
a
b
NF
BF
KC
RKC
ARC
3030
330
6060
3030
330
6060
2700
660
167,400
2700
660
167,400
5400
3360a
18,900b
89
Table 6
Average results and number of feasible solutions found for the column generationbased heuristic with and without violation allowed in terminals for the harder class of
instances generated.
Formulation
Time (seconds)
gapUB (%)
Feasible found
11.5
13.2
0.12
0.40
16
14
Min
XX X
2F c Dp;t xc;D
p;t
p2P t2Tpc2Cp
D
XX
min
max
C v iol v iolt;d v iolt;d
38
t2T d1
max
X X
c;dDp;t
V c xp;t
min
8t 2 T; d 2 f1;.. .;Dg 40
p2Ptc2Cp
Fig. 2. (a) Representation of a set of z variables corresponding to the same route and
tanker class. (b) Branching on z variables. (c) Branching on x variables.
branching. This observation is depicted in Fig. 2(b). This characteristic, when replicated in several parts of the solution, generates an undesirable symmetry that prevents nishing the
branch-and-bound enumeration even in instances where the root
gap is very small. As shown in Fig. 2(c), the use of x variables
breaks this symmetry. Because of that, the ARC formulation
was tested in two ways: with the z variables as binaries
(ARCbin ), and with these variables as continuous between zero
and one (ARCcont ). Comparing both ways in Tables 2 and 3, one
can see that the ARCcont formulation shows a slight improvement
regarding the average gapLB at the root node in medium and hard
instances, and a little improvement concerning the average nal
gapLB in the hard instances. Thus, although the difference between
the two assumptions for the z variables is not expressive, we consider that the ARCcont formulation obtained the best results for the
instances of the medium and hard classes.
The results obtained for all formulation for the harder class of
instances are shown in Table 5. For this class of instances, we test
the original ARC and a modied version of this formulation. The
modication consists of allowing both inventory shortages and
maximum inventory capacities violations in the consumption sites,
keeping the violation cost high enough to avoid changing the optimal solution. For this modication, we include two new sets of
continuous variables:
max
For the harder class of instances, the GCD is always one. Rounding the inventory just drops the fractional part of the accumulating
production or consumption. So, we show only the results of the
rounding formulations. To nd optimal or near optimal solutions
for all instances of this class (to be used as the UBb value), we
run each instance for many hours, using as initial feasible solution
the best solution found in previous tests. Again the ARC has the
best performance. Comparing this formulation with and without
violation, we can see a good improvement in the quality of the
feasible solutions. This occurs because, when allowing violation,
the number of feasible solutions increases, facilitating the CPLEX
heuristic procedures to construct and improve feasible solutions.
Only the knapsack cascading formulations nd feasible solutions
for all instances, but the average gapUB of the ARC is better because
when this formulation nds a feasible solution for an instance, it
usually improves the solution faster than the knapsack cascading
formulations.
Concerning the DWM, it obtained the best average root relaxation, even when the CPLEX preprocessing routines and cuts are
considered for the other formulations. Without CPLEX preprocessing and cuts, the average value of the linear programming relaxation for all formulations is 8.82.
The results obtained for the column generation-based heuristic
are shown in Table 6 for the harder class of instances. Column
Feasible Found shows the number of instances where the heuristic nds a feasible solution, and column gapUB (%) shows the
average gapUB for these solutions. We also test the heuristic including violation in the ARC. Comparing the CGH with and without
violation, we can see a good improvement when the violation is
allowed. This occurs because the number of feasible solutions increases. The average gapUB for both cases shows that, as expected,
the heuristic nds good feasible solutions. For the instances where
the heuristic nds feasible solutions, the gapUB of the ARCv iol formulation is 0.69. Moreover, this formulation reaches feasibility
for 6 out of the 9 instances where the heuristic fails to nd feasible
solutions. However, the average time to obtain a feasible solution
Table 5
Average results and number of instances solved for the mathematical formulations for the harder class of instances generated.
Formulation
Time (seconds)
gapLB root lp
gapUB (%)
No. of nodes
Inst. solved
BF
RKC
ARC
ARCv iol
DWM
720.00
720.00
637.64
642.36
8.49
8.26
7.81
7.81
4.74
7.50
6.53
5.88
5.89
6.37
5.34
5.16
5.19
10.9
6.20
1.85
562,356
107,204
137,416
110,018
0
0
3
3
90