![Formulations For A Problem of Petroleum Transportation [PDF]](https://pdfs.asia/img/200x200/formulations-for-a-problem-of-petroleum-transportation.jpg)
11 0 582 KB
European Journal of Operational Research 237 (2014) 82–90
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Discrete Optimization
Formulations for a problem of petroleum transportation Luiz Aizemberg ⇑, Hugo Harry Kramer, Artur Alves Pessoa, Eduardo Uchoa Production Engineering Department, Fluminense Federal University, Rua Passo da Pátria 156, 24210-240 Niterói, RJ, Brazil
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 finite 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 find 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 fleet 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 finite 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: [email protected] (L. Aizemberg), hugoharry @gmail.com (H.H. Kramer), [email protected] (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
horizon. The products are shipped only by oil tankers and we assume that the global capacity of the fleet is unlimited. For each transportation, one has to decide the route, the lot size and the delivery day. This problem has already been addressed in Rocha (2010) and Rocha, Grossmann, and Poggi de Aragão (2011). This paper has two aims. The first aim is to compare mathematical formulations for the transportation problem proposed in Rocha et al. (2011). Firstly, a basic formulation is given, where the inventory level for each day in each site is represented by a continuous variable and the shipments are represented by binary variables. This formulation is similar to many others found in the literature. We refer to Sawik (2011) for an overview of such papers. We also test several variations of the formulation proposed in Rocha et al. (2011), referred to as Knapsack Cascading. Lastly, we devise a new formulation that outperforms all previous ones when given to a Mixed Integer Programming (MIP) solver. The second aim of this paper is to propose a column generationbased heuristic approach. We note that, for difficult instances, the commercial solver used is not capable of finding good feasible solutions at the beginning of the branch-and-bound algorithm tree. Finding a good feasible solution at the very beginning of the tree is crucial to reduce the search by cutting off nodes of such tree, increasing the probability of proving the optimality in a reasonable computational time. We did our work on the top of a previous paper found in the literature (Rocha et al., 2011). The instances used to test the formulations are those proposed by Rocha et al. Our paper shows a clear evolution concerning the results as the formulations are improved, leading to the generation of new and harder instances (based on the existent ones) to evaluate the potential of such formulations.
83
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
1.1. Literature review Several models for optimizing the transportation of petroleum products are found in literature, considering different characteristics. In Colvin and Maravelias (2011), the authors consider task selection, scheduling, and resource planning decisions, exploring risk management methods, and develop a multi-stage stochastic programming model. To optimize such models, Mixed Integer Programming is the most used technique. In Magatão, Arruda, and Neves (2004) and Boschetto et al., 2010, MIP is applied to schedule the activities in a real oil pipeline network. The former uses illustrative instances with four products, while the latter uses large instances, where more than 14 oil derivatives and ethanol are transported. In Moro and Pinto (2004), MIP is used to optimize the short-term crude oil scheduling. The authors deal with the separation of brine from the oil, interface separation between different types of oil, and energy consumption in the crude distillation units. MIP is also used in a real-time approach to dispatch petroleum tank trucks (Brown & Graves, 1981) and to optimize the long-term oil derivative scheduling under uncertainty (MirHassani, 2008). Finally, in Banaszewski, Tacla, Pereira, Arruda, and Enembreck (2010) and Aizemberg et al. (2011), the authors deal with the mediumterm planning of oil derivatives transportation in the Brazilian oil industry multimodal supply chain network. The former applied an auction based multiagent algorithm, while the latter applied MIP. Column generation based heuristics (CGH) are found in literature to deal with many different problems: multi-item scheduling (Bahl, 1983), vehicle routing (Mourgaya & Vanderbeck, 2007), cutting stock (Furini, Malaguti, Durain, Persiani, & Toth, 2012), sensor placement (Yavuz, Aras, Kuban, & Ersoy, 2010), generalized assignment (Moccia, Cordeau, Monaco, & Sammarra, 2009) and many others. In general, such approaches share similar main characteristics. First, the linear programming relaxation is solved by column generation. Then, a feasible solution is constructed for the original problem using the obtained columns. For example, in Choi and Tcha (2007) the authors apply a column generation heuristic for the heterogeneous fleet vehicle routing problem which consists of three steps: (i) solve the continuous relaxation; (ii) set the variables of the final restricted master as binaries; and (iii) solve the resulting MIP problem. In Joncour, Michel, Sadykov, Sverdlov, and Vanderbeck (2010), the authors review generic classes of column generation-based heuristics. 1.2. Paper organization The remainder of the paper is organized as follows. Section 2 presents the problem definitions and the mathematical formulations, including the Dantzig–Wolfe decomposition formulation used by the proposed heuristic. Section 3 shows the column generation-based heuristic procedure, including dynamic programming algorithm for the pricing subproblem. Section 4 discusses the obtained results. At last, Section 5 contains the conclusions of this work.
in classes according to their capacities. Let C be the set of all tanker classes. It is assumed that only oil tankers with fulfilled capacities are used in the transportation. The objective is to minimize the total transportation costs. The formulations described in the next subsections consider the following notation: Indices:
p: Platform, p 2 P; t: Terminal, t 2 T; c: Class of tanker, c 2 C; d: day, d 2 f1; . . . ; Dg.
Subsets: TðpÞ # T: Terminals allowed to send tankers to a platform p 2 P; PðtÞ # P: Platforms allowed to receive tankers from a terminal t 2 T; CðpÞ # C: Classes of tankers allowed to pick up oil from platform p 2 P. Input data:
ep;0 : Initial inventory at platform p 2 P; et;0 : Initial inventory at terminal t 2 T; Pp;d P 0: Production at platform p 2 P, in day d; C t;d P 0: Consumption at terminal t 2 T, in day d; CAPp;d : Maximum inventory capacity at platform p 2 P, in day d; CAPt;d : Maximum inventory capacity at terminal t 2 T, in day d; V c : Capacity of a tanker of class c 2 C; F c : Transportation cost per day for a tanker of class c 2 C (tankers must travel back from the platform to the terminal after each transportation); Dp;t : Transportation time between a platform p 2 P and a terminal t 2 T in days.
Other definitions will be added as needed by each formulation. 2.1. Natural formulation For this model the decision variables are defined as follows. zc;d p;t 2 f0; 1g is a binary variable indicating if a tanker of class c 2 C is assigned to ship the product from platform p 2 P to terminal t 2 T at day d. The variables ep;d 2 Rþ and et;d 2 Rþ represent the inventory levels at platform p 2 P and terminal t 2 T at day d, respectively. Thus, the natural formulation (NF ) is given by the following: D XX XX c;d 2F c Dp;t zp;t ;
ðNF Þ min subject to
ep;d ¼ ep;d1 þ P p;d
2. Mathematical formulations
X X
c;d V c zp;t
8p 2 P; d 2 f1;...; Dg;
ð2Þ
t2TðpÞc2CðpÞ
et;d ¼ et;d1 C t;d þ
This section presents the mathematical formulations considered for the problem in study. The optimization problem defined 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
ð1Þ
p2P t2TðpÞc2CðpÞ d¼1
X X
c;dDp;t
V c zp;t
8t 2 T; d 2 f1;...;Dg; ð3Þ
p2PðtÞc2CðpÞ
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;
8p 2 P; t 2 TðpÞ; c 2 CðpÞ; d 2 f1;. ..;Dg:
ð4Þ
ð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
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
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 definitions listed on the following. Define 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:
2.2. Knapsack cascading formulation For this third formulation, which is adapted from Rocha et al. (2011), some of the previous definitions are considered plus others that follow.
DAp;d ¼ AC p;d CAPp;d ; DAt;d ¼ CAPt;d AC t;d : Thus, the Knapsack Cascading (KC) formulation is given by:
ðKCÞ min
D XX XX 2F c Dp;t zc;d p;t ;
ð9Þ
p2P t2TðpÞc2CðpÞ d¼1
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Þ
t2TðpÞ s¼1 c2CðpÞ 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Þ
p2PðtÞ s¼1 c2CðpÞ
s¼1
Also, let kmax ðp; dÞ be the maximum number of lots of size Lp allowed to be shipped until day d from platform p avoiding its inventory level to be less than zero. Such number is calculated by:
AC p;d ; kmax ðp; dÞ ¼ Lp With the definitions given above, constraints (4) are replaced with:
AC p;d kmax ðp; dÞLp 6 ep;d 6 AC p;d kmin ðp; dÞLp d 2 f1; . . . ; Dg:
8p 2 P; ð7Þ
Now, define Lt as the GCD of all the shipping capacities of all classes of tankers that are allowed to supply terminal t. Let kmin ðt; dÞ be the minimum number of lots of size Lt that must be received by terminal t until day d such that its inventory level is not negative. This number is given by:
AC t;d ; Lt
where AC t;d is the accumulated amount of product at terminal t from day 0 to day d and is given by:
AC t;d ¼ et;0
ð8Þ
t2TðpÞ s¼1 c2CðpÞ
where AC p;d is the accumulated amount of product at platform p from day 0 to day d and is given by:
kmin ðt; dÞ ¼
8t 2 T; d 2 f1;. ..;Dg:
With the substitutions suggested above, the second formulation for the problem, which will be called Basic Formulation (BF ), is given by the objective function (1) subject to constraints (2), (3), (7), (8) and (6).
AC p;d CAPp;d ; kmin ðp; dÞ ¼ Lp
AC p;d
AC t;d þ kmin ðt;dÞLt 6 et;d 6 AC t;d þ kmax ðt;dÞLt
d X C t;s :
s¼1
Let kmax ðt; dÞ be the maximum number of lots of size Lt allowed to be received by terminal t until day d without exceeding its inventory capacity. Such quantity is calculated by the following formula:
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Þ
p2PðtÞ s¼1 c2CðpÞ
zc;d p;t 2 f0;1g
8p 2 P; t 2 TðpÞ;c 2 CðpÞ; d 2 f1;... ;Dg: ð14Þ
The objective function (9) minimizes the total transportation cost. Constraints (10) avoid the inventory levels at each platform to be less than the minimum capacity. Constraints (11) ensure that the inventory levels at each platform never exceed their maximum inventory capacity. Constraints (12) avoid the inventory levels at each terminal to be less than the minimum capacity and constraints (13) ensure that the inventory levels at each terminal never exceed their maximum inventory capacity. Finally, the z variables are defined by (14) as binaries. From the KC formulation, it is possible to derive a fourth formulation replacing constraints (10) and (11) by the new constraints (15) and (16) that follow: d X XX
V c zc;p;ts 6 kmax ðp;dÞLp
8p 2 P; d 2 f1;... ;Dg
ð15Þ
V c zc;p;ts P kmin ðp;dÞLp
8p 2 P; d 2 f1;...;Dg
ð16Þ
t2TðpÞ s¼1 c2CðpÞ d X XX t2TðpÞ s¼1 c2CðpÞ
and substituting (12) and (13) by (17) and (18) below: d X XX
c;sDp;t
V c zp;t
6 kmax ðt;dÞLt
8t 2 T; d 2 f1;...;Dg
ð17Þ
P kmin ðt;dÞLt
8t 2 T; d 2 f1;...;Dg:
ð18Þ
p2PðtÞ s¼1 c2CðpÞ d X XX
c;sDp;t
V c zp;t
p2PðtÞ s¼1 c2CðpÞ
The resulting formulation will be called Rounded Knapsack Cascading (RKC) formulation.
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
formulation, proving that it has the same linear relaxation solution set as the RKC over the z variables. h
2.3. Accumulated rounded cascading formulation The variables of this newly proposed formulation are defined as follows: xc;d p;t is an integer variable representing the number of tankers of class c 2 CðpÞ leaving platform p 2 P to terminal t 2 TðpÞ from day 1 to day d. The zc;d p;t variables are also considered as previously defined. Let xc;0 p;t ¼ 0. The Accumulated Rounded Cascading (ARC) formulation is given by: ðARCÞ min
XX X
2F c Dp;t xc;D p;t ; p2P t2TðpÞc2CðpÞ
85
The same proof can be applied for NF ; KC and AC, where AC is the accumulated cascading formulation without rounding, leading to the following proposition: Proposition 2.2. NF ; KC and AC have the same linear relaxation.
ð19Þ
2.4. Dantzig–Wolfe decomposition
subject to c;d c;d1 xp;t ¼ xp;t þ zc;d 8p 2 P; t 2 TðpÞ;c 2 CðpÞ; d 2 f1;...;Dg; p;t X X c;d V c xp;t 6 kmax ðp;dÞLp 8p 2 P; d 2 f1;...;Dg;
ð20Þ ð21Þ
t2TðpÞc2CðpÞ
X X
c;d V c xp;t P kmin ðp;dÞLp
8p 2 P; d 2 f1;...;Dg;
ð22Þ
t2TðpÞc2CðpÞ
X X
c;dDp;t
6 kmax ðt;dÞLt
8t 2 T; d 2 f1;...;Dg;
ð23Þ
c;dDp;t
P kmin ðt;dÞLt
8t 2 T; d 2 f1;...;Dg;
ð24Þ
V c xp;t
p2PðtÞc2CðpÞ
X X
V c xp;t
p2PðtÞc2CðpÞ c;d 2 Zþ xp;t
8p 2 P; t 2 TðpÞ; c 2 CðpÞ; d 2 f1;...;Dg;
c;d zp;t 2 f0;1g or ½0;1
8p 2 P; t 2 TðpÞ; c 2 CðpÞ; d 2 f1;...;Dg:
ð25Þ
ð26Þ
The objective function (19) is to minimize the total transportation cost. Constraints (20) define 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) define the x variables as nonnegative integers and constraints (26) define 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 definitions. Proposition 2.1. relaxation.
BF ; RKC and ARC have the same linear
Proof. We will prove the proposition in two steps. First, we will show that RKC is obtained from BF without changing its linear relaxation solution set over the z variables. Second, we will show that ARC is obtained from RKC without changing its linear relaxation solution set over the z variables. Now, consider the BF . Given the constraints (2) and expanding ep;d1 , we will have:
ep;d ¼ AC p;d
d X XX
V c zc;p;ts
8p 2 P;
d 2 f1; . . . ; Dg
t2TðpÞ s¼1 c2CðpÞ
Here, we present a Dantzig–Wolfe decomposition (DWD) (Dantzig & Wolfe, 1960) for the problem using the (RKC) formulation as starting point. The resulting reformulation is primarily used to derive improved lower bounds for the original problem via Column Generation (CG). The results obtained with this formulation are used as a benchmark to analyze the linear relaxation of the implemented formulations of this paper. Furthermore, we observe that the stronger linear relaxation provided by DWD helps to select z variables that are suitable for building good feasible solutions. Motivated by this observation, we propose a CG based heuristic in Section 3 which demonstrates the importance of the DWD formulation. The following notation is considered for the reformulation: let Sp be the set of all feasible solutions for a platform p 2 P, let V t;d p;s be the amount of product shipped from a platform p 2 P to terminal t 2 TðpÞ at a day d in a solution s 2 Sp , and let C p;s be the cost of a solution s 2 Sp for a platform p 2 P. Binary variables kp;s indicate if a solution s 2 Sp for a platform p 2 P is used. The reformulation obtained after applying the Dantzig–Wolfe decomposition, the so called Dantzig–Wolfe Master ðDWMÞ, is the following: ðDWMÞ min
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 s¼1 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 s¼1 d XXX
V t;p;ssDp;t kp;s P
ð29Þ
ð30Þ
p
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) define the k variables as binaries.
From (7),
AC p;d kmax ðp;dÞLp 6 AC p;d
d X XX
V c zc;p;ts
6 AC p;d kmin ðp;dÞLp
t2TðpÞ s¼1 c2CðpÞ
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. Define the x variables as follows:
xc;d p;t ¼
d X zc;p;ts
s¼1
and replace the right-hand side of the previous equality by xc;d p;t in the RKC. Including the constraints (20), we obtain the ARC
2.4.1. Pricing subproblem Given the large number of variables in DWM formulation, its direct use is impractical for large instances. To overcome this issue, a Column Generation (GC) procedure is often used in order to solve the linear programming relaxation of DWM. This procedure is started by considering only a subset of its variables. The problem with such subset of variables is called the Restricted Master Problem (RMP). At each iteration of the CG procedure, the RMP is solved and the dual variables of its optimal solution are used to generate new k variables by solving the Pricing Subproblem. A new k variable is added to RMP if its reduced cost, given by the objective function value of the Pricing Subproblem, is negative. If, at a given iteration, it is no longer possible to generate a k variable with negative reduced cost, the procedure is finished and the current solution for RMP is the optimal solution for the linear relaxation of DWM.
86
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
Let pp be the vector containing the dual variables associated to constraints (28), ht;d be the vector containing the dual variables associated to constraints (29), and rt;d be the vector containing the dual variables associated to constraints (30). The binary variable zc;d p;t indicates if the platform p sends a tanker of class c to terminal t at day d. Thus, for each platform p, a Pricing Subproblem is defined as follows.
ðSP p Þ min
D X XX c;d ð2F c Dp;t V c ðht;dþDp;t þ rt;dþDp;t ÞÞzp;t pp ; ð32Þ t2TðpÞc2CðpÞ d¼1
subject to d X XX
Fig. 1. Dynamic programming with two lot sizes.
zc;p;ts V c 6 AC p;d 8d 2 f1; .. . ; Dg;
ð33Þ
t2TðpÞ s¼1 c2CðpÞ d X XX
zc;p;ts V c P DAp;d
8d 2 f1; . . . ;Dg;
ð34Þ
t2TðpÞ s¼1 c2CðpÞ
zc;d p;t 2 f0; 1g
8t 2 TðpÞ; c 2 CðpÞ; d 2 f1; . . . ;Dg:
ð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 defined as binary by (35). Even though the Pricing Subproblem is presented as an Integer Linear Program, it can be solved more efficiently 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 definitions we assume that the elements of the set C are arbitrarily numbered from 1 to jCj. We define 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
P 1;p;c;d;e ¼ fa 2 f0; 1gjTj j0 6 e þ Pp;d
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 first and the second lot size, and the initial inventory is 2 units. The first 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 first 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
P 3;p;c;d;e ¼ fa 2 f0; 1gjTj j0 6 e
X V c at 6 CAPp g; t2T
and C c;d p;t ¼ 2F c Dp;t V c ðht;dþDp;t þ rt;dþDp;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
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
which is equal the original capacity if the accumulated production is integer, or equal the original capacity minus one, if the accumulated capacity is fractional. The fifth column shows the adjusted production, which is equal the integer part of the accumulated production not yet added to the inventory. 3.2. Column generation-based heuristic procedure The proposed heuristic uses the column generation procedure described in SubSection 2.4 with the pricing algorithm presented in the previous subsection to choose only some z variables of ARC with good chance to appear in the optimal solution. This reduction in the number of variables allows the branch-and-bound algorithm to find fast a feasible solution, or prove that no solution can be found with these variables. Such Column Generation-based heuristic is shown in Algorithm 1. First, the linear relaxation of the DWM formulation is solved by column generation (line 1). Next, the reduced costs of all columns present in the final RMP are calculated using the values of the dual variables in the last iteration of the column generation algorithm (line 2). As these columns must have non-negative reduced costs, we then sort such columns in ascending order according to their reduced costs (line 3). The parameters min and max are initialized in line 4. In line 5, the parameter k is initialized. This parameter indicates that the kth first columns of the final RMP must be converted into variables of ARC. At each iteration of the main loop (lines 6–16), the ARC formulation is populated with only the z variables corresponding to the k columns with least positive reduced costs in the final RMP (line 7). The ARC formulation run for up to five seconds (line 8). If a new feasible solution is obtained (line 9), we keep the solution if it improves the best found so far. If the branch-and-bound enumeration was not completed within the given five seconds (line 12), we finish the loop. Otherwise, we increase k by ten (line 15). It may happen that some platform has no chosen z variable in the reduced ARC formulation. If this is the case, the resulting MIP will be trivially infeasible and a larger value for k will be tried. When k is larger than the number of RMP columns with reduced cost zero, every platform will have at least one converted column. The time limit used in line 8 has been chosen based on preliminary experiments where we observed that most best feasible solutions found by the algorithm were attained when solving a MIP exactly in less than five seconds. Thus, this time limit seems to give a good compromise between solution quality and runtime. Algorithm 1. Column Generation based heuristic 1: Run column generation to solve linear relaxation of DWM 2: Calculate the reduced cost of the columns using the dual variables of the last iteration 3: Sort the columns according to their reduced costs 4: min 30; max 150 5: k min 6: while k 6 max do 7: Convert k columns of DWM with least positive reduced costs into variables of ARC 8: Run ARC for up to five seconds 9: if feasible solution found then 10: Keep the solution if the best found 11: end if 12: if enumeration incomplete within five seconds then 13: Exit while 14: end if 15: k k þ 10 16: end while
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:
only one lot size (easy); one lot size per route (medium); two lot sizes per route (hard); five lot sizes per route and all platforms are allowed to send product to all terminals (harder new instances).
In this work, the results obtained for the first 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 final (%) 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
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
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
gapLB root (%)
gapLB final (%)
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
gapLB root (%)
gapLB final (%)
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 – 1
1.47 and 3.19. Without CPLEX preprocessing, DWM would have the best gapLB rootlp for all instance classes. Yet, for the hard class of instances, DWM found the best average gapLB at the root node even when the CPLEX preprocessing routines and cuts are used for the other formulations. In Tables 2 and 3 one can observe that the NF and BF have the worst performance in comparison with the others. The gapUB of them is not shown because these formulations found few feasible solutions, so average value is not appropriate. The BF was capable of solving only two instances to optimality, both of the medium class. This fact impacts in the average execution time for this formulation, which is close to the maximum allowed for each instance. For the hard class of instances, as this formulation was not able to solve any instance to optimality, the average time is exactly equal to the maximum allowed for each instance. The NF did not solve any instance, so the average time is exactly 720 for both classes of instances. On the other hand, since these are more compact formulations with respect to the remaining ones in the sense that it has less nonzero coefficients in the matrix, a significantly larger number of nodes were solved at the branch-and-bound algorithm. Comparing KC with RKC and AC with ARCcont , the performance are similar. We can conclude that formulation NF is dominated by BF but the effect of the inventory rounding in KC and AC is not noticeable, even for GCD of lot sizes greater than one. In comparison with BF , the RKC formulation achieved a better performance, being able to solve 13 out of 50 instances, where 11 instances belong to the medium class and 2 are hard instances. This occurs because in the RKC formulation the constraints (10)–(13) are recognized by CPLEX as knapsack constraints (for (11) and (12) all binary variables are replaced with their complements to derive equivalent knapsack constraints). As a result, CPLEX is able to add specialized cuts, e.g. cover cuts. This effect can be observed by comparing the root lower bounds of BF and RKC in Tables 2 and 3. The best performance was obtained by the ARC formulation, outperforming the literature results reported in Rocha et al. (2011). This formulation was able to solve 24 instances of the medium class and 22 hard instances in a significantly small average
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
80% of these constraints are in constraints (20). 43% of these variables are in constraints (20).
execution time. For the 4 instances not solved by this formulation, it found 3 of the best known feasible solutions. For the 49 remaining instances, the gapUB is equal to zero. Thus, the average gapUB is equal to zero for the medium class of instances and close to zero for the hard instances. The success of this formulation can be explained by three main reasons: (i) the number of nonzero coefficients in ARC is much smaller than in RKC (and comparable to that of BF ) because each x variable in the constraints (21)–(24) replaces a sum of many z variables; (ii) as in the RKC, CPLEX recognizes the constraints (21)–(24) as knapsack constraints; and (iii) the x variables are more suitable for branching. To explain the item (i), we observe that the number of nonzero coefficients of NF ; BF and ARC is OðjPkTkCjDÞ while for KC and RKC is OðjPkTkCjD2 Þ. Table 4 compares the sizes of all formulations for an instance of the harder class. To explain the item (iii), we refer to Fig. 2. The graph of Fig. 2(a) represents the z variables associated to one route for the whole time horizon: the arcs represent the variables, the first column of nodes represents the source platform and the second column of nodes represents the destination terminal. Now, assume that z15 has a fractional value in a relaxed solution for the RKC formulation. If z15 is selected as a branching variable and forced to assume the value one, the relaxed solution will probably change significantly. However, when z15 is forced to be zero, it is likely that the relaxation gives the same lower bound if there is sufficient slack in the inventories of both p1 and t 1 to send the same amount of product one day latter, i.e. to set z16 to the same value as z15 had before
89
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
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
CGH with violation CGH no violation
11.5 13.2
0.12 0.40
16 14
Min
XX X
2F c Dp;t xc;D p;t þ
p2P t2TðpÞc2CðpÞ
D XX min max C v iol v iolt;d þ v iolt;d
ð38Þ
t2T d¼1 max
We also include the new variables v iolt;d in the constraints (23) min and v iolt;d in the constraints (24), as follows: X X c;dD max V c xp;t p;t v iolt;d 6 kmax ðt;dÞLt 8t 2 T; d 2 f1;...;Dg ð39Þ p2PðtÞc2CðpÞ
X X
c;dDp;t
V c xp;t
min
þ v iolt;d P kmin ðt;dÞLt
8t 2 T; d 2 f1;.. .;Dg ð40Þ
p2PðtÞc2CðpÞ
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 finishing 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 final 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 modified version of this formulation. The modification 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 modification, we include two new sets of continuous variables: max
v iolt;d : maximum capacity inventory violation at terminal t and day d, min v iolt;d : inventory shortage at terminal t and day d.
Let C v iol be the cost of inventory violations. As we do not want to change the optimal solution, we use a large enough value. The new objective function will be as follows:
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 find 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 find feasible solutions for all instances, but the average gapUB of the ARC is better because when this formulation finds 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 finds 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 finds good feasible solutions. For the instances where the heuristic finds 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 find 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
gapLB root (%)
gapLB final (%)
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
L. Aizemberg et al. / European Journal of Operational Research 237 (2014) 82–90
in these cases is 245.9 seconds. The average run time of both variants of the heuristic is low, with a minimum of 7 seconds and a maximum of 16.78 seconds. One may think that, as the knapsack cascading formulation finds feasible solutions for all instances, the heuristic should be used with this formulation. But our tests show that the accumulated cascading formulation has a better performance. This happens because the heuristic decreases the variable set, so it is easier to find a feasible solution. Comparing both formulations, when ARC finds a feasible solution, usually the solution is better than RKC. 5. Conclusions This work presented a modeling study for the oil transportation problem. Three classes of instances taken from the literature were used to test the formulations. Given the good results obtained, a new class was generated with harder instances. Six formulations were devised: the Natural Formulation (NF ), the Basic Formulation (BF ), the Knapsack Cascading Formulation (KC), the Rounded Knapsack Cascading Formulation (RKC), the Accumulated Rounded Cascading Formulation (ARC), and the Dantzig–Wolfe Master Formulation (DWM). The first five formulations were tested with the CPLEX branch-and-cut algorithm, and the best results found used the ARC formulation and its variants. Without the CPLEX cuts, the DWM formulation achieved the best lower bounds, which motivated the development of a column generation-base heuristic. In comparison with the results reported by Rocha et al. (2011), the results obtained in this work show a good improvement. This achievement is due only to a better understanding of the problem, incorporated in the mathematical formulation, reducing the computational effort. Rounding the inventories improves the formulation only for NF , not being noticeable when applied in RKC and ARC. The harder class of instances has been used to evaluate the ability of the best formulations to find good feasible solutions when the optimality could not be proved. In this case, the ARC formulation achieved the best performance when allowing inventory capacity violation in the demand sites, but keeping the violation cost high enough to ensure that the cost of the optimal solution will not change. The column generation-based heuristic used for the harder class of instances found good feasible solutions in a small computational time. We showed that it is also possible to improve the quality of the feasible solutions found by this heuristic allowing inventory capacity violations. References Aizemberg, L., Uchoa, E., Pessoa, A., Rocha, R., Coutinho, R., & Paula, U. C., Jr. (2011). Optimisation Model for the problem of oil transportation with local search by
MIP and simulation. In Proceedings of XLIII SBPO 2011, Ubatuba–SP, Brazil, in Portuguese. Bahl, H. C. (1983). Column generation based heuristic algorithm for multi-item scheduling. IIE Transactions, 15(2), 136–141. http://dx.doi.org/10.1080/ 05695558308974624. Banaszewski, R. F., Tacla, C. A., Pereira, F. R., Arruda, L. V., & Enembreck, F. (2010). Planning transport of crude oil derivatives with simultaneous auctions. In 2010 IEEE international conference on systems man and cybernetics (SMC) (pp. 2820–2827). IEEE. Boschetto, S. N., Magatão, L., Brondani, W. M., Neves, F., Jr., Arruda, L. V. R., BarbosaPóvoa, A. P. F. D., et al. (2010). An operational scheduling model to product distribution through a pipeline network. Industrial & Engineering Chemistry Research, 49(12), 5661–5682. Brown, G., & Graves, G. (1981). Real-time dispatch of petroleum tank trucks. Management Science, 27(1), 19–32. Choi, E., & Tcha, D.-W. (2007). A column generation approach to the heterogeneous fleet vehicle routing problem. Computers & Operations Research, 34(7), 2080–2095. Colvin, M., & Maravelias, C. T. (2011). R&D pipeline management: Task interdependencies and risk management. European Journal of Operational Research, 215(3), 616–628. http://dx.doi.org/10.1016/j.ejor.2011.06.023. ISSN 0377-2217. Dantzig, G. B., & Wolfe, P. (1960). Decomposition principle for linear programs. Operations Research, 8(1), 101–111. Furini, F., Malaguti, E., Durain, R. M., Persiani, A., & Toth, P. (2012). A column generation heuristic for the two-dimensional two-staged guillotine cutting stock problem with multiple stock size. European Journal of Operational Research, 218(1), 251–260. http://dx.doi.org/10.1016/j.ejor.2011.10.018. ISSN 0377-2217. Ghiani, G., Laporte, G., & Musmanno, R. (2004). Introduction to logistics systems planning and control (Vol. 14). Wiley. Joncour, C., Michel, S., Sadykov, R., Sverdlov, D., & Vanderbeck, F. (2010). Column generation based primal heuristics. Electronic Notes in Discrete Mathematics, 36(0), 695–702. http://dx.doi.org/10.1016/j.endm.2010.05.088. ISSN 15710653. Magatão, L., Arruda, L. V. R., Neves, F., et al. (2004). A mixed integer programming approach for scheduling commodities in a pipeline. Computers & Chemical Engineering, 28(1–2), 171–185. MirHassani, S. (2008). An operational planning model for petroleum products logistics under uncertainty. Applied Mathematics and Computation, 196(2), 744–751. Moccia, L., Cordeau, J.-F., Monaco, M. F., & Sammarra, M. (2009). A column generation heuristic for a dynamic generalized assignment problem. Computers & Operations Research, 36(9), 2670–2681. http://dx.doi.org/10.1016/ j.cor.2008.11.022. ISSN 0305-0548. Moro, L., & Pinto, J. (2004). Mixed-integer programming approach for short-term crude oil scheduling. Industrial & Engineering Chemistry Research, 43(1), 85–94. Mourgaya, M., & Vanderbeck, F. (2007). Column generation based heuristic for tactical planning in multi-period vehicle routing. European Journal of Operational Research, 183(3), 1028–1041. http://dx.doi.org/10.1016/ j.ejor.2006.02.030. ISSN 0377-2217. Rocha, R. (2010). Petroleum supply planning: models, reformulations and algorithms, Ph.D. thesis, Departamento de Informática, PUC-Rio. Rocha, R., Grossmann, I. E., & Poggi de Aragão, M. V. S. (2011). Cascading knapsack inequalities: Reformulation of a crude oil distribution problem. Annals of Operations Research, 1–18. Sawik, T. (2011). Scheduling in supply chains using mixed integer programming. Wiley. ISBN 9781118029091. . Yavuz, B., Aras, N., Kuban, I., & Ersoy, C. (2010). A column generation based heuristic for sensor placement, activity scheduling and data routing in wireless sensor networks. European Journal of Operational Research, 207(2), 1014–1026. http:// dx.doi.org/10.1016/j.ejor.2010.05.020. ISSN 0377-2217.
