Optimization of the Petroleum Product Supply Chain under Uncertainty

Article pubs.acs.org/IECR

Optimization of the Petroleum Product Supply Chain under Uncertainty: A Case Study in Northern Brazil Fabrício Oliveira and Silvio Hamacher* Department of Industrial Engineering, Pontifícia Universidade Católica, Rio de Janeiro, Brazil S Supporting Information *

ABSTRACT: This Article addresses the problem of optimizing the investment planning process of a logistics infrastructure for the distribution of petroleum products under uncertainty. For this purpose, a two-stage stochastic model was developed. Because of the large number of possible scenarios for the future realization of demands, the Sample Average Approximation (SAA) methodology was used to produce approximations of the optimal solution. The proposed model and solution methodology were applied to a real case study on the distribution of petroleum products in northern Brazil. The results show that the methodology allows us to obtain solutions that are relatively close to the real optimal solution of the problem with optimality gaps estimated at up to 1%, which are reasonable for most practical purposes.

1. INTRODUCTION For almost 50 years, companies in the oil and chemical industries have led the development and use of linear mixed integer programming to support decision making at all levels of planning.1 Nevertheless, it is noteworthy that most of the applications in the field have not considered uncertainty in their evaluations.2 An overriding feature in the oil industry is its wide range of uncertainties, which are typically related to the unpredictable levels of demand for refined products, fluctuations in prices in domestic and international markets, and inaccuracies in the forecasted production of oil and gas. For this reason, some works have used techniques based on mathematical programming to support decision making under uncertainty in the oil industry (e.g., stochastic optimization).3−7 However, the incorporation of uncertainty into these models using stochastic optimization remains a challenge due to the large computational requirements involved.2 The predominant focus of the aforementioned articles was on the tactical and operational planning of refining activities. When it comes to the particular aspects related to the strategic planning of logistics of the oil industry, as investments in pipelines, distributions bases, and ports, it can be noticed that they are usually modeled in a simplified fashion. We believe that this represents a gap in the context of optimizing the industry operation in a broader point of view, especially because there is an inherent dependence between the continuity of the refining process and the distribution capability of the network. Aiming to fulfill this gap, we propose a mathematical model for the optimization of oil product distribution that considers issues of tactical planning (i.e., aggregate flows of products, inventory levels, and production levels among others) to evaluate decisions of a strategic nature (i.e., investments in logistics infrastructure). The proposed model is based on two-stage stochastic programming, where the first-stage decisions represent the decisions regarding investments (on terminals, bases, and pipelines) and the secondstage decisions are related to the distribution planning. To obtain solutions, we used a sampling strategy that was proposed in the seminal work of Shapiro and Homem-de-Mello,8 which is known in the literature as Sample Average Approximation © 2012 American Chemical Society

(SAA). The technique attempts to obtain approximate solutions for the total set of scenarios using a simulation-based approach. This approach has been used in the literature as a method of avoiding the difficulty of dealing with a large number of scenarios. Linderoth et al.9 applied the SAA technique to several linear programming models. Kleywegt et al.10 presented important theoretical considerations regarding the method for combinatorial problems and illustrated them with numerical examples of some applications. Verweij et al.11 demonstrated the application of SAA to routing problems with large numbers of scenarios (up to 21694) and obtained solutions with optimality gaps of approximately 1.0%. Santoso et al.12 proposed an application that was specifically aimed at supply chain design and applied it to a real study of the beverage industry. More recently, Schütz et al.13 applied the SAA methodology to a supply chain design problem for the supply chain of food products. The first significant contribution of this Article involves the development of a mathematical model for the optimal logistics planning of oil product distribution considering issues of tactical planning to evaluate decisions of a strategic nature. For this purpose, we detail the specific features of the logistics of petroleum products distribution and the possible investments that could be made in its infrastructure. Another important contribution is our consideration of the uncertain nature of the problem in a more comprehensive way than is usually found in the literature for planning the supply chain of petroleum and its byproduct through the use of SAA. The authors are not aware of an existing application of such a technique to the optimization of the petroleum products chain. The final contribution of this Article is its focus on a case study that is oriented toward a region of Brazil with very particular characteristics regarding the distribution of petroleum products. The implementation of the proposed approach led to important conclusions regarding the planning of logistics Received: Revised: Accepted: Published: 4279

June 22, 2011 January 25, 2012 February 8, 2012 February 8, 2012 dx.doi.org/10.1021/ie2013339 | Ind. Eng. Chem. Res. 2012, 51, 4279−4287

Industrial & Engineering Chemistry Research

Article

refinery availability, and the sources of production. Transportation and investment decisions are defined in conjunction, and the latter are chosen from a predefined portfolio of possibilities and allocated across the planning horizon. The uncertainties in the model are related to the levels of demand for petroleum products in the distribution bases, which are modeled as random variables. The demurrage cost is approximated by a piecewise linear curve such that the linearity of the model is preserved. The precision level of such approximation is directly proportional to the number of linear segments used for the representation of the original curve. However, there is a trade-off that must be considered because when one uses a greater number of these segments, the number of variables and constraints in the mathematical model increases. Given that the focus of the model is geared toward the commercialization, logistics, and allocation of petroleum products, aspects related to the petroleum refining process are not considered as objects of this decision model but as aspects that have already been decided in advance. Therefore, only the bases and the terminals of the model in question are subject to investment decisions. Each investment depends on the following factors individually or in combination: • Changes in operational costs, justified by the introduction of more modern and/or more energy-efficient equipment • Changes in the rotation capacity of the tanks, investments in more powerful pumps that allow faster pumping (i.e., tanks are filled and emptied more quickly) • Changes in demurrage cost curves, investing in ship management technology and/or adding new piers for berthing or modernizing existing berths • Changes in storage capacity, the installation of additional tanks for storage The most typical cases of investments in transportation arcs consist of improvements or expansions of the existing pipeline network, increasing pumping capacity, reducing operating costs, and/or changing the types of products that can be transported by that arc.

investments in the region. Again, we are not aware of the existence of previous academic papers that are concerned with this area of study. This Article is organized as follows: Section 2 provides a detailed description of the problem to be addressed. The developed mathematical model is presented in section 3. Section 4 gives an overview of the case study. The SAA methodology used and the setup of the experiments are explained in section 5. The computational results and an analysis of these results are presented in section 6. Finally, in section 7, some conclusions are drawn.

2. PROBLEM DESCRIPTION The problem in question can be defined as the strategic planning of petroleum products distribution, determining the levels of investments in logistics infrastructure, the distribution of flows, inventory policies, and the level of the external commercialization of refined products. These decisions are made in a multiperiod setting and are discretized into uniform intervals. Typically, a logistics network for oil products is composed of a set of nodes (international markets, refineries, terminals, and bases) that are connected by transportation arcs. A transportation arc is defined by its origin node, destination node, the mode of transport, and product group. The product group defines which products are compatible with the operational specifications of the arc and, therefore, can be transported by it. Some specific arc modes, such as pipelines, are reversible, which means that it can operate in both directions using the same structure. After they are produced in refineries, products are stored in tanks to be directed to distribution bases. These bases serve as negotiation points with distributors and are considered as aggregation points of demand for such products. The tanks of these bases are constantly being loaded and unloaded. This process is known as the tank rotation and is subject to the physical limitations that are inherent in the hardware associated with the tanks of the distribution base. The rotating capacity refers to the number of times a tank can be filled and emptied over a certain period of time. Many of the products are transported by sea through marine terminals. Therefore, it is important to consider the issue of the time that the ships spend in port while waiting for permission to dock and unload, which increases the logistical costs due to the demurrage of the vessel at the site. This demurrage cost can be modeled by curves obtained from simulation models. The demurrage curves are drawn proportionally to the amounts that are handled at the port. These curves are typically nonlinear and depend mostly on the operational characteristics of each marine terminal. Normally, an increase in port workload leads to the augmented congestion of vessels in the harbor, which increases logistical costs. To address this issue, we propose a two-stage stochastic mixedinteger linear programming model. The first-stage decisions are which projects to implement and when they are to be implemented; the second-stage decisions are those related to the flow of products, inventory, levels of supply provided to each demand site, amounts to be internally and externally commercialized, and the supply levels at the sources. The purpose of the model is to provide the optimal distribution of refined products to meet the demand of distribution bases, thereby minimizing the logistical costs of this operation and maximizing the revenue from commercialization with external markets. Meeting internal demand is mandatory in the present context. Moreover, supplying internal demand depends on the characteristics of the network operations,

3. MATHEMATICAL FORMULATION The objective of the mathematical model is to minimize the costs related to investment, freight, inventory, operations, demurrage in marine terminals, and commercialization and penalties for nonfulfillment of demand (backlog costs), which are subject to constraints relating to network operational conditions, the demand for oil products, and supply limits. To provide greater clarity in the problem’s notation, the domains of summations are omitted except when the summation is evaluated only on a subset of the natural domain. When there is no mention of this fact, its domain should be considered as the original set to which the index refers. 3.1. First-Stage Problem. The first-stage model comprises the decisions regarding when and where the investments should be made prior to the realization of the uncertainty. The firststage model can be written as follows: minimize ∑ CKAatyat + y,w

a,t

∑ CKLlt wlt + ΕΩ[Q (y , w , ξ)] i,t (1)

∑ yat t 4280

≤ 1 ∀ a ∈ AK (2) dx.doi.org/10.1021/ie2013339 | Ind. Eng. Chem. Res. 2012, 51, 4279−4287


∑ wlt ≤ 1 ∀ l ∈ LK

Article

olt,,pξ ≤ OIlt, p ∀ l ∈ L , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

(3)

t

t,ξ 3.2.3. Bound on Domestic Supply. The total supplied ol,p t must comply with the supply limit OIl,p .

Equation 1 refers to the expenses related to investments in infrastructure (arcs and places, respectively) that the model decides to deploy. The term EΩ[Q(y,w,ξ)] represents the expected value among the evaluation of the second-stage costs under each possible realization ξ ∈ Ω of the uncertain parameters. Constraints (2) and (3) indicate that only a single investment can be made in each arc or location along the planning horizon. 3.2. Second-Stage Problem. The second-stage problems represent the costs of supplying the demand given an investment decision (y,w) for a certain realization ξ of the uncertain parameters. The second-stage problem Q(y,w,ξ) is formulated as follows.

∑

minimize

CFDatxdat ,, pξ+

∑

⎛ xiat ,, ξp ⎞⎟ t , ξ ⎜ ∑ Va , p⎜xda , p + ≤ Q a(1 − ∑ yat′) Ia ⎟ ⎝ ⎠ p t ′≤ t + QKa

−

∑

PElt, pelt,,pξ + θ

∑

vlt,,pξ ≤ R l , p(1 −

3.2.1. Objective Function. Terms + ∑a,p,t t,ξ CFIat xia,p compute the total freight cost, that is, the cost of moving products over the arcs in the direct and reverse t,ξ directions, respectively. Term ∑l,g,p∈g,t CIgvl,p stands for the costs of maintaining the inventories of the product considering t,ξ the level at the end of each period t. Terms ∑l,g,t COlt zl,g + t,ξ ∑l,g,t COKlt zkl,g denote the operations cost of each site considering the current costs and the possible changes in t t,ξ operating costs due to investments. Terms ∑S,l,t CSS,l f S,l + t t,ξ ∑S,l,t CSKS,l f kS,l calculate the demurrage cost to be incurred for locations that have port access also considering the current values and possible changes due to investments in infrastructure t t,ξ t or operational improvements. Terms ∑l,p,t PIl,p il,p − ∑l,p,t PEl,p t,ξ el,p figure the net importation and exportation. Finally, the last t,ξ term θ∑l,p,t Sl,p represents the costs incurred by not meeting the demand.

∑ a|j=l

=

∑ a|i=l

∑

t,ξ xda,p

zlt,,pξ ≤ GR lt, pR l , p(1 −

∀ l ∈ LB ∪ LS , ∀ t ∈ T , ∀ ξ ∈ Ω zklt,,pξ ≤ GRKlt, pRKl , p

wlt ′

∀ l ∈ LB ∪ LS , ∀ t ∈ T , ∀ ξ ∈ Ω

zlt,,pξ + zklt,,pξ ≥

∑

(10)

(xdat ,, ξp + xiat ,, ξp)

a|j=l

∑

+


a|i=l

∀ l ∈ LB ∪ LS , ∀ t ∈ T , ∀ ξ ∈ Ω

(11)

3.2.6. Movement of Products through Each Location. t,ξ Constraint (9) limits the total throughput zl,p by the product of t the tank referential rotation GRl,p and the total available tankage Rl,p. In the event that the model chooses an investment wlt at time t, the throughput, which from the moment t becomes t,ξ represented by zkl,p , is subject to the new reference rotation t GRKl,p and the new tankage RKl,p as explained in constraint (10). Constraint (11) consolidates the flow of inputs and outputs (direct and reverse flows) in the respective variables. The t,ξ binary characteristic of investment variables implies that zl,p =0 t,ξ ∨ zkl,p = 0,∀l ∈ LB ∪ LS, ∀p ∈ P, ∀t ∈ T, ∀ξ ∈ Ω.

xiat ,, ξp + Dlt,,pξ + vlt,,pξ − slt,,pξ ∀

a|j=l

l ∈ LB ∪ LR ∪ LS , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

∑

(9)

t ′≤ t

1, ξ xiat ,, ξp + olt,,pξ + vlt,− p

∑

wlt ′)

∑ t ′≤ t

a|i=l

xdat ,, ξp +

t ′≤ t

3.2.5. Storage Capacity. The total tankage used for the storage of product p must not exceed the limit Rl,p. If an investment wlt is made at that location at some previous time t′ ≤ t, the site operates at its new storage capacity RKl,p.

(4)

xdat ,, ξp +

wlt ′

(8)

slt,,pξ

∑a,p,t CFDat

∑

∀ l ∈ LB ∪ LR ∪ LS , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

l ,p,t

l ,p,t

wlt ′) + RKl , p

∑ t ′≤ t

PIlt, pilt,,pξ

l ,p,t

s,l ,t

s,l ,t

(7)

3.2.4. Limit on the Capacity of Arcs. The total capacity Q a of movement along arc a is shared by the direct and inverse flows of all of the products that pass along it. The flow volume is adjusted by a viscosity factor Va,p and a reversal factor Ia in the case of inverse flows. In the event that the model chose an investment yat to be made in arc a at a prior time t′ ≤ t in the planning horizon, the arc operates based on its new capacity QKa.

CFIatxiat ,, pξ

∑ CSst, l fst,,lξ + ∑ CSKst, lfkst,,lξ + ∑

yat′ ∀ a ∈ A , ∀ t ∈ T , ∀ ξ ∈ Ω

∑ t ′≤ t

xd , xi , v , z , zk , f , fk , i , e , s a , p , t a,p,t + ∑ CIg vlt,,pξ + ∑ COlt zlt,,gξ + ∑ COKlt zklt,,gξ l ,g ,t l ,g ,p∈g ,t l ,g ,t

+

(6)

(5)

t,ξ 3.2.2. Material Balance. All direct flows ∑a|j=l xda,p and t,ξ inverse flows ∑a|i=l xia,p along arcs into location l plus the total t,ξ amount ol,p supplied by location l and the total remaining t−1,ξ inventory from the previous period vl,p should be equal to the t,ξ t,ξ and ∑a|j=l xia,p , sum of all direct and inverse flows (∑a|i=l xda,p t,ξ respectively) out of location l plus the local demand Dl,p and t,ξ the total inventory vl,p while discounting the eventual total of t,ξ t−1,ξ unmet demand sl,p in the period t in question. For t = 1, vl,p = 0 Vl,p . This constraint applies only to bases, refineries, and places with a marine terminal.

∑ a|i=l

xdat ,, ξp +

∑

xiat ,, ξp = elt,,pξ

a|j=l

∀ l ∈ LE , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω 4281

(12)

dx.doi.org/10.1021/ie2013339 | Ind. Eng. Chem. Res. 2012, 51, 4279−4287


Article

Figure 1. Distribution network.

elt,,pξ ≤ LElt, p ∀ l ∈ LE , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

linearization was performed along the demurrage curve dividing it into |S| segments s each with its respective linear cost. Constraints (16) and (17) bound the amount fs,lt,ξ handled for each segment. Given an investment wlt at time t, the amount handled, which t,ξ from the moment t becomes represented by fkS,l , is subject to the new limit FSKs,lt . Constraint (18) consolidates the terminal input flows undertaken by the maritime modal m ∈ MM ⊆ M in the respective handling variables. The binary characteristic of investment variables implies that fs,lt,ξ = 0 ∨ fks,lt,ξ = 0,∀s ∈ S,∀l ∈ LS, ∀t ∈ T, ∀ξ ∈ Ω.

(13)

∑

xdat ,, ξp +

a|j=l

xiat ,, ξp = ilt,,pξ

∑ a|i=l

∀ l ∈ LE , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

ilt,,pξ ≤ LIlt, p ∀ l ∈ LE , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

(14) (15)

3.2.7. External Commercialization. For foreign trade partners, the amounts traded are represented by specific variables: t,ξ el,pt,ξ represents exports, and il,p represents imports. Constraints (12) and (14) define the flows to and from foreign locations l as these variables. Constraints (13) and (15) bound the amounts of exports and imports, respectively. f st,,lξ ≤ FSst, l(1 −

∑

4. CASE STUDY The model was applied to a real case study on refined products distribution in northern Brazil. The transport in the region is primarily performed using waterway modals, which are strongly affected by seasonality issues regarding the navigability of rivers. For this study, four different products were considered, diesel, gasoline, aviation fuel, and fuel oil, to be distributed over 13 bases, three of which have sea terminals. Three supply sources were considered including one refinery and two external supply ́ (SP) and locations. The external supply, coming from Paulinia São Luiz (MA), represents the connection of the regional logistics network under study with the rest of the country. The case study does not include international commercialization. Four distinct modes of transportation are considered including waterways (using large ferries and small boats), roadways, and pipelines. Waterway transportation is generally performed by large ferries, which are typically used during periods of river flooding, and by smaller boats, which are able to navigate the rivers during droughts (i.e., low water level seasons); however, this form of transportation is costly. Figure 1 schematically represents the network under study. The region considered comprises approximately 3.7 MMkm2, which represents nearly 43% of Brazil’s national territory. As shown in this figure, the bases of Manaus (AM), Itacoatiara (AM), Santarém (PA), Macapá (AP), and Belém (PA) are particularly

wlt ′)

t ′≤ t

∀ s ∈ S , ∀ l ∈ LS , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω fkst,,lξ ≤ FSKst, l

∑

(16)

wlt ′

t ′≤ t

∀ s ∈ S , ∀ l ∈ LS , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

∑ f st,,lξ s

+

∑ fkst,,lξ = s

∑

(17)


a | j = l ∧ m ∈ MM , p

∀ l ∈ LS , ∀ p ∈ P , ∀ t ∈ T , ∀ ξ ∈ Ω

(18)

3.2.8. Demurrage. The amount of product handled in marine terminals, which is subject to demurrage costs, is modeled by constraints (16)−(18). The demurrage costs increase nonlinearly as a function of the amount handled. Therefore, a piecewise 4282



Article

site where each investment will be conducted and the type of project. Three distinct types of investments are considered at each location including investments in storage capacity that increase the location’s capacity for processing and storing a given product, investments in pumps and substations that reduce the operating costs and increase the ability to rotate the tanks, and investments in the construction of a new pier, which make the demurrage cost curves per handled volume smoother. The investment portfolio also has an investment available for the implementation of a multiproduct pipeline that connects the bases of Porto Velho and Rio Branco. The planning horizon considered was 8 years, which are divided into a total of 32 quarterly periods. All of the costs considered in the model are discounted to a present value under a yearly interest rate of 6.8%.

relevant because they act also as distribution points of the supply coming from São Luiz (MA). Depending on the season, these arcs may or may not be available for navigation. Table 1 shows how the seaworthiness Table 1. Seaworthiness between Locations origin Itacoatiara

destiny Itaituba

Itacoatiara

Porto Velho

Manaus

Caracarai ́

Manaus Base

Cruzeiro do Sul

Manaus

Itaituba

Manaus

Santarém

Santarém

Porto Velho Itaituba

Porto Velho

mode ferries small boats ferries small boats ferries small boats ferries small boats ferries small boats ferries small boats ferries small boats ferries small boats

first quarter

second quarter

X

X

third quarter

fourth quarter

X

X X

X X

X X X

X

X X X X

X

X

X

5. SOLUTION METHODS To take into account the uncertainty in demand levels for petroleum products, scenarios were generated by the following first-order autoregressive model:

X

X

1 Dlt, p = Dlt,− p [1 + ωp + σpε], t = 2, ..., |T |

X X

X

X X

where ωp represents the expected average growth rate for the consumption of product p over the planning horizon, σp represents the estimated maximum deviation of the product p consumption, and ε follows a standard normal distribution. The estimate of the maximum deviation used was simplified as being identical for each product due to the lack of data regarding the historical consumption of the products in the region studied. This estimation was made on the basis of an analysis of the annual Brazilian petroleum products consumption series over the last 40 years.14 Each scenario represents a possible demand curve for the entire time horizon considered and for each product and distribution base in the considered problem. It is not computationally feasible to consider all of the possible scenarios for demand levels due to the large number of possible realizations for the uncertain parameters. For instance, if discrete limits were defined for the levels of demand, the number of scenarios would be given by N|P|x|LB|, where N is the number of such levels, |P| is the total number of products, and |LB| is the number of distribution bases. In the case study considered, the definition of two possible values for demand levels would require the model to consider 252 scenarios, which might become computationally infeasible. To circumvent this difficulty, a scenario sampling scheme based on the Monte Carlo simulation was used. First proposed by Shapiro and Homem-de-Mello,12 this scheme is known in the literature as Sample Average Approximation (SAA). The idea of SAA is to repeatedly solve the problem (1−17) for a reduced set of scenarios. For M independent random samples of size N, an approximation of the value of the recourse costs is obtained. The original problem is then approximated by the following SAA problem:

X X

X X

was modeled in various parts of the region under study. The “X” cells represent periods in which the given stretch is unavailable for navigation by that mode. Observe that during certain times of the year, the Cruzeiro do Sul base remains completely isolated from communication with the system (third quarter), while the base of Caracarai ́ can only rely on supply via the roadway mode during the first two quarters of the year. Figure 2 shows the level of demand for each of the bases. Manaus (AM) is the main hub of the region’s demand, followed by Porto Velho (RO), Belém (PA), and Macapá (AP).

⎧ ⎪ minimize ⎨gN̂ (y , w) = ∑ CKAatyat + ∑ CKLlt wlt ⎪ a,t l ,t ⎩

Figure 2. Levels of demand.

The portfolio of projects considered in the study consists of 28 local projects and one arc project. Such projects are considered mutually independent and can therefore be combined. Table 2 represents the portfolio of investments considered, showing the

+

4283

1 N

⎫ ⎪ Q (y , w , ξn)⎬ ⎪ n = 1,..., N ⎭

∑

(19)



Article

Table 2. Investment Portfolio for Locations locations project type

Manaus

Macapá

Santarém

Belém

Cruzeiro do Sul

Itacoatiara

X X X X X X

X X X

X X X X X X

X X X X

X X X

X X X

diesel tankage gasoline tankage aviation fuel tankage fuel oil tankage pumps and substations pier

X

where Q(y,w,ξn) is the second-stage problem to be solved given a realization ξn of the uncertainty set. Linderoth et al.9 showed that, using this approach, it is possible to obtain upper and lower quotas for the optimal value of the real problem and, furthermore, that such boundaries converge to the true value with probability 1 as N increases. Kleywegt et al.10 showed that for combinatorial problems, assuming that the SAA problem is solved up to a optimality gap δ, the sample size required to ensure that the attainment of ε optimality with probability 1−α can be estimated as:

=

Therefore, one can estimate the required number of scenarios, once given a sampling estimator σ̂n, a desired confidence interval, and a confidence level 1 − α as: ⎛ ⎞2 z α /2 σ̂N ⎟⎟ N = ⎜⎜ ̂ (y , w) ⎠ ⎝ (β/2)gN

⎛ 2|AK |×|LK |×|T | ⎞ ⎜⎜ ⎟⎟ log α (ε − δ)2 ⎝ ⎠ 3σ2max (ε − δ)2

[(|LK | + |AK |)(|T |)log 2 − log α]

6. RESULTS The mathematical model and the scenario generation routines were implemented in AIMMS 3.10.16 The mixed-integer linear programming (MILP) model was solved using CPLEX 11.2. Table 3 describes the size of the models for the case study in

(20)

where ε ≥ δ, α ∈ [0,1], and 2|AK|×|LK|×|T| represent the total number of possible first-stage solutions. In estimate (20), the 2 term σmax is defined as the maximal variance of certain function differences in the optimal solution (y*,w*) (the reader should refer to Kleywegt et al.10 for further details regarding estimate 20). Even though estimate (20) can be too conservative for practical applications,10 it suggests that the sample size required to reach complete convergence grows at most linearly with the number of projects in the portfolio. Another possible way for estimating the number of scenarios is to define a desired accuracy level of the solution, which can be measured by the confidence interval of the expected total cost.8,15 Such a confidence interval can be defined as: ⎡ z α /2 σ̂N z α /2 σ̂N ⎤ , gN̂ (y , w) + ⎢gN̂ (y , w) − ⎥ ⎣ N N ⎦

Table 3. Summary of the Size of the Modela

a

N

σ̂N =

(gN̂ (y , w) − (∑ CKAatyat +

∑

n=1

a,t 2 n + Q (y , w , ξ )))

l ,t

N

total variables

total constraints

average solution time (s)

standard deviation of solution time (s)

20 30 40

250 400 455 504 606 864

304 144 374 880 499 360

532.83 971.42 1472.05

967.26 1500.20 864.41

50 replications.

question together with the mean and standard deviation of the solution time for solving each SAA problem. All of the experiments were performed using a Pentium Quad-Core 2.6 GHz with 8 GB RAM. To obtain estimates of the upper and lower limits, experiments were performed with N equal to 20, 30, and 40. These values for N were defined approximating the true values obtained using the estimate of Monte Carlo sampling standard deviation (eq 22) for N = 50(σ̂50) considering β = 0.1, α1 = 0.05, α2 = 0.025, and α3 = 0.01. The solution times ranged from 532.83 s for instances with 20 scenarios to 1472.05 s for those with 40 scenarios. To obtain the lower limits, we performed 50 replications (i.e., M = 50), with a time limit of 3600 s and a relative gap of 1% defined as the stop criteria. Thirty-six candidate solutions were generated for N = 20, 22 for N = 30, and 19 for N = 40. To obtain the upper limit, all of the candidate solutions were previously assessed with 50 replications. From this first evaluation, the three solutions that showed the best results in terms of solution gap were refined and subsequently evaluated with 1000 replications.

(21)

where N is the number of scenarios, σ̂n is the Monte Carlo sampling standard deviation estimator,8 and zα/2 is the standard normal deviate such that P(z ≤ zα/2) = 1 − α/2, where z follows a standard normal distribution. Assuming that one can solve problem (19) for N scenarios, the estimate σ̂N is given by:

∑

(23)

where β ∈ [0,1]. The term (β/2)ĝN(y,w) represents the fraction of the total cost one wishes to consider as the as the absolute confidence interval value for each side. For example, if one desires a confidence interval of ±5% of the total cost, then β = 0.1. In practical terms, the choice of the number of scenarios should take into account the trade-off between the computational effort to obtain a solution and the quality level required for the solution.

3σ2max

N≥

X

CKLlt wlt

N−1 (22) 4284



Article

Table 4 shows the best results for each experiment in terms of the lower and upper limits estimated for the solution of the

Table 6 provides the solutions with the lowest gap obtained from the three experiments including solution number 3 for

Table 4. Results of Experiments: Upper and Lower Statistical Limits

Table 6. Investment Profiles of Solution 3 for N = 20, Solution 3 for N = 30, and Solution 2 for N = 40

N 20

30

40

lower limit

upper limit

800.12 9.81 1.2% 801.25 10.22 1.2% 805.28 8.22 1.0%

818.76 109.66 13.40% 821.67 50.63 6.20% 817.12 40.03 4.90%

amount (106$) standard deviation (106$) percentage deviation amount (106$) standard deviation (106$) percentage deviation amount (106$) standard deviation (106$) percentage deviation

investment period project Inv. Inv. Inv. Inv. Inv. Inv.

Table 5. Results of Experiments: Estimative of the Optimality Gap and Its Statistical Upper Limit gap N

solution number

value (10 $)

%

standard deviation (106$)

20

1 2 3 1 2 3 1 2 3

19.61 26.4 18.64 23.18 26.95 20.42 17.38 11.83 14.85

2.40% 3.20% 2.30% 2.80% 3.30% 2.50% 2.10% 1.40% 1.80%

107.41 72.82 110.1 63.73 82.4 51.57 44.03 41.22 49.43

30

40

N = 20

N = 30

N = 40

7 21 24 1 29 23

7 16 22 1 24 27

6 17 16 1 27 26

N = 20 and N = 30, and solution number 2 for N = 40. As can be observed from Table 6, the profile of investments has little variability between experiments. An investment was made in fuel oil storage for Santarem in the first period in all runs, which indicates the attractiveness of this investment. This may be explained by the central position of Santarém in the petroleum products distribution network of the region and by the low level of tankage for fuel oil currently available at that location. Other investments also tend to have low variability in their positioning along the time horizon. The greatest variability was seen in the investment in pumps and substations in Macapá, which is directly related to the existence of an anticipated increase in demand for fuel oil in Belo Monte that is transported from São Luiz. The solutions also suggest that Santarém is a strategic location for the logistics of products other than fuel oil because the model suggests investing in three tanking projects in the region. Another relevant observation is related to the projects that did not constitute the optimal portfolio. None of the projects for the physical expansion of the marine terminals (piers) is selected for investment, which suggests that the terminal system as modeled is appropriate to the scenarios of demand for the products considered. The pipeline connecting the Porto Velho and Rio Branco bases turned out to be not economically attractive and was not included in the optimal portfolio of investments in any of the simulated scenarios. This is probably because of the high cost of that project and the ́ that, existence of an alternative road coming from Paulinia despite its high costs, is more economically efficient. All of the demand was satisfied in all experiments.

real problem. The results suggest that the configuration of the experiment with 50 replications (M = 50) for the lower bound was considered satisfactory given that the deviation obtained for the lower limit is approximately 1%. For the upper limit, it should be noted that its variability is related to the number of scenarios considered in obtaining the lower limit, and it is reduced from 13.4% (N = 20) to 4.9% (N = 40). Table 5 shows the statistics obtained on the estimate of the optimality gap and the upper limits for the three best solutions

6

Manaus Aviation Fuel Santarem Diesel Santarem Gasoline Santarem Fuel Oil Belém Diesel Macapa BS

obtained in each experiment. The optimality gaps obtained range from 3.3% to 1.4%, with no strong evidence of improvement for experiments with a larger number of scenarios. However, the experiments suggest a reduction of the estimated variability of the gap for the experiments with larger number of scenarios, which supports the hypothesis that these solutions are close to the real optimal solution of the problem. In practical terms, this optimality gap is considered acceptable, given the uncertainty inherent in the input data that are considered deterministic. This result is noteworthy, especially because this estimator is biased;8 thus, such an estimate always corresponds to an upper limit of the real gap. Upon analyzing the results of simulations to obtain the upper limit, the concept of solution resilience becomes clearer. Typically, a larger number of scenarios implies a more comprehensive investment profile in terms of its ability to cope with higher demands, which makes the system more robust with respect to variations in the total costs of meeting the demand (i.e., smaller fluctuations in second-stage costs), suggesting that there is indeed an increase in the accuracy of these solutions.

7. CONCLUSIONS This Article presented the problem of investment planning in the supply chain of petroleum product distribution in northern Brazil considering the uncertain demand for such products in the region. The mathematical model presented can be applied to any problem in the supply chain of crude oil and its byproducts and considers new and/or restructuring investments for the supply chain. Specific operational issues were modeled (e.g., navigability of transportation systems, operational capacity of rotating tanks, and demurrage costs, among others) that affect particular details of the logistics chain. Given the size of the possible set of scenarios for constructing the uncertainty of the demand, the SAA method was used. The results show that it was possible to obtain solutions with acceptable estimates of optimality gaps in practical terms (i.e., in terms of the solution quality and the computational time required to obtain the solutions) even with a relatively small number of scenarios. It is worth mentioning that if discrete limits were defined 4285



Article

ξ∈Ω

for the levels of demand, the number of scenarios would be given by N|P|x|LB|, where N is the number of such levels, |P| is the total number of products, and |LB| is the number of distribution bases. Even for a small number of levels N, the problem might become computationally infeasible. The proposed approach allows us to consider an infinite set (however limited) of realizations of demand and, therefore, of possible scenarios. In the proposed approach, it is possible to delineate reasonably acceptable confidence intervals and thereby define the total number of scenarios required to statistically guarantee that the solutions obtained. It is important to highlight that the amount of scenarios required ais strongly related to the variability of the recourse cost of the particular instance considered as can been seen in eq 20. The case study showed that from the proposed portfolio, only six projects comprise the optimal portfolio of investments. The results suggest that the Santarém region has a particular strategic importance for the planning as one-half of these investments were assigned to that region. Another important observation is the finding that many projects in the set of possible investments were not relevant to the optimization of the logistics in the region for the data set considered. The results seem to be in line with what has been observed in the literature9−13 when it comes to the successful application of the SAA methodology to solve practical large-scale problems. It can be observed that, even for a modest number of scenarios (ranging from 20 to 40, in this case), the method can provide high-quality solutions with relatively small optimality gaps. Therefore, the proposed methodology can support the decision making process, while identifying solutions and statistically ensuring its quality without the need for time-consuming discussions of the adequacy of possible methods for generating scenarios.

■

AI ⊆ A

AK ⊆ A LB ⊆ L LE ⊆ L LF ⊆ L LR ⊆ L LS ⊆ L LK ⊆ L MM ⊆ M Parameters

CFDa t CFIat CIg CKAat CKLlt COlt COKlt

ASSOCIATED CONTENT

S Supporting Information *

t CSS,l

Sample average approximation. This material is available free of charge via the Internet at http://pubs.acs.org.

■

t CSKS,l

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

t,ξ Dl,p

Notes

The authors declare no competing financial interest.

■

t FSS,l

ACKNOWLEDGMENTS We gratefully acknowledge the support and encouragement of the Brazilian National Council of Research-CNPq and Coordination for the Improvement of Higher Education Personnel in BrazilCAPES. We are also grateful to the reviewers whose insightful comments helped to improve the presentation of this Article.

■

set of scenarios

Subsets

t FSKS,l

t GRl,p

NOMENCLATURE

Sets i,j,l ∈ L set of locations p,p1,p2 ∈ P set of products g = {(p1,p2,...,pk)∈Pk} ∈ G set of product groups m∈M set of transportation modes a={(i,j,g,m) ∈ L × L set of arcs × G × M} ∈ A s∈S set of segments of the piecewise linearization of the demurrage cost function t∈T set of periods

t GRKl,p

Ia t LEl,p t LIl,p

4286

subset of arcs that can operate in reverse subset of arcs subject to project subset of bases subset of international markets subset of production sources subset of refineries subset of nodes with access via marine terminal subset of locals subject to projects subset of shipping modes ∀a ∈ A,∀t ∈ T, freight cost per unit of product transported by the arc during period t ∀a ∈ AI,∀t ∈ T, reverse freight cost per unit of product transported by the arc during period t ∀g ∈ G, holding cost of product group g ∀a ∈ AK,∀t ∈ T, cost of investment in the project for arc a in period t ∀l ∈ LK,∀t ∈ T, cost of investment in the project at location l for period t ∀l ∈ LF ∪ LS,∀t ∈ T, cost of operation of the site during period t ∀l ∈ LK,∀t ∈ T, cost of operation after completion of project at location l during period t ∀s ∈ S,∀l ∈ LS,∀t ∈ T, demurrage cost in segment s for period t in location l ∀s ∈ S,∀l ∈ LS ∩ LK,∀t ∈ T, demurrage cost in segment s after the completion of the project for location l during period t ∀l ∈ LB,∀p ∈ P,∀t ∈ T,∀ξ ∈ Ω, demand in location l for product p in period t for realization ξ ∀s ∈ S,∀l ∈ LS,∀t ∈ T, maximum volume of segment s for location l during period t ∀s ∈ S,∀l ∈ LS ∩ LK,∀t ∈ T, maximum volume of segment s after the project conclusion for location l during period t ∀l ∈ LB ∪ LS,∀p ∈ p, ∀t ∈ T, reference tank rotation at location l for product p during period t ∀l ∈ (LB ∪ LS) ∩ LK,∀p ∈ P,∀t ∈ T, reference tank rotation after the project conclusion for location l for product group g in period t ∀a ∈ AI,∀t ∈ T, inversion factor for arc a ∀l ∈ LE,∀p ∈ P,∀t ∈ T, export limit for product p for location l in period t ∀l ∈ LE,∀p ∈ P,∀t ∈ T, import limit for product p for location l in period t dx.doi.org/10.1021/ie2013339 | Ind. Eng. Chem. Res. 2012, 51, 4279−4287

Industrial & Engineering Chemistry Research t OIl,p

t PEl,p

t PIl,p

Qa QKa Rl,p RKl,p Va,p 0 Vl,p

θ Variables

fs,lt,ξ

f ks,lt,ξ

t,ξ el,p

t,ξ il,p

t,ξ ol,p

t,ξ sl,p

t,ξ vl,p

t,ξ xda,p

t,ξ xia,p

Article

∀l ∈ LF ∪ LR,∀p ∈ P,∀t ∈ T, supply of product p at location l in period t ∀l ∈ LE,∀p ∈ P,∀t ∈ T, sales price in location l of product p in period t ∀l ∈ LE,∀p ∈ P,∀t ∈ T, purchase price at location l of product p in period t ∀a ∈ A, current capacity of arc a ∀a ∈ AK, capacity after project completion for arc a ∀l ∈ LB ∪ LR ∪ LS,∀p ∈ P, storage capacity of location l for product p ∀l ∈ (LB ∪ LR ∪ LS) ∩ LK,∀p ∈ P, storage capacity after project conclusion in location l for product p ∀a ∈ A,∀p ∈ P, viscosity factor in arc a for product p ∀l ∈ LB ∪ LR ∪ LS,∀p ∈ P, initial inventory of product p in location l penalty for not meeting the demand

yat ∈ {0,1}|a ∈ AK, t ∈ T, decision to implement a project for arc a in period t wlt ∈ {0,1}|l ∈ LK, t ∈ T, decision to implement a project for location l in period t t,ξ zl,p ∈ 9 + |l ∈ LB, ∪ LS, P ∈ P, t ∈ T, ξ ∈ Ω, amount handled at location l of products p in period t under realization ξ t,ξ zkl,p ∈ 9 + |l ∈ LB ∪ LS, P ∈ P, t ∈ T, ξ ∈ Ω, amount handled after the completion of the project for location l of products in group g for period t under realization ξ

yat wlt t,ξ zl,p

t,ξ zkl,p

■

REFERENCES

(1) Shapiro, J. Challenges of strategic supply chain planning and modeling. Comput. Chem. Eng. 2004, 28, 855−861. (2) Sahinidis, N. Optimization under uncertainty: state-of-the-art and opportunities. Comput. Chem. Eng. 2004, 28, 971−983. (3) Escudero, L.; Quintana, F.; Salmerón, J. CORO, a modeling and an algorithmic framework for oil supply, transformation and distribution optimization under uncertainty. Eur. J. Oper. Res. 1999, 114, 638−656. (4) Dempster, M.; Hicks Pedron, N.; Medova, E.; Scott, J.; Sembos, A. Planning logistics operations in the oil industry. J. Oper. Res. Soc. 2000, 51, 1271−1288. (5) Al-Othman, W.; Lababidi, H.; Alatiqi, I.; Al-Shayji, K. Supply chain optimization of petroleum organization under uncertainty in market demands and prices. Eur. J. Oper. Res. 2008, 189, 822−840. (6) Khor, C.; Elkamel, A.; Ponnambalam, K.; Douglas, P. Two-stage stochastic programming with fixed recourse via scenario planning with economic and operational risk management for petroleum refinery planning under uncertainty. Chem. Eng. Process. 2008, 47, 1744−1764. (7) Ribas, G.; Hamacher, S.; Street, A. Optimization under uncertainty of the integrated oil supply chain using stochastic and robust programming. Int. Trans. Oper. Res. 2010, 17, 777−796. (8) Shapiro, A.; Homem-de-Mello, T. A simulation-based approach to two-stage stochastic programming with recourse. Math. Program. 1998, 81, 301−325. (9) Linderoth, J.; Shapiro, A.; Wright, S. The empirical behavior of sampling methods for stochastic programming. Ann. Oper. Res. 2006, 142, 215−241. (10) Kleywegt, A.; Shapiro, A.; Homem-de-Melo, T. The sample average approximation method for stochastic discrete optimization. SIAM J. Optim. 2002, 12, 479−502. (11) Verweij, B.; Ahmed, S.; Kleywegt, A.; Nemhauser, G.; Shapiro, A. The sample average approximation method applied to stochastic routing problems: A computational study. Comput. Optim. Appl. 2003, 24, 289−333. (12) Santoso, T.; Ahmed, S.; Goetschalckx, M.; Shapiro, A. A stochastic programming approach for supply chain network design under uncertainty. Eur. J. Oper. Res. 2005, 167, 96−115. (13) Schütz, P.; Tomasgard, A.; Ahmed, S. Supply chain design under uncertainty using sample average approximation and dual decomposition. Eur. J. Oper. Res. 2009, 199, 409−419. (14) Empresa de Pesquisa Energética (EPE) - Balanço Energético Nacional (BEN). Ministério de Minas e Energia, 2010; available at: https://ben.epe.gov.br/downloads/Relatorio_Final_BEN_2010.pdf. (15) You, F.; Wassick, J. M.; Grossmann, I. E. Risk management for a global supply chain planning under uncertainty: Models and algorithms. AIChE J. 2009, 55, 931−946. (16) Bisschop, J.; Entriken, R. AIMMS The Modeling System; Paragon Decision Technology: Haarlem, The Netherlands, 1993.

fs,lt,ξ ∈ 9 + |s ∈ S,l ∈ LS,t ∈ T,ξ ∈ Ω, amount handled in segment s subject to demurrage at location l in period t under realization ξ f ks,lt,ξ ∈ 9 + |s ∈ S,l ∈ LS ∪ LK,t ∈ T, ξ ∈ Ω, amount handled in segment s subject to demurrage after project completion for location l, product p, period t under realization ξ t,ξ el,p ∈ 9 + |l ∈ LE, p ∈ P, t ∈ T, ξ ∈ Ω, amount exported from location l of product p in period t under realization ξ t,ξ il,p ∈ 9 + |l ∈ LE, p ∈ P, t ∈ T, ξ ∈ Ω, amount imported to location l of product p in period t under realization ξ t,ξ ol,p ∈ 9 + |l ∈ LF ∪ LR, p ∈ P, t ∈ T, ξ ∈ Ω, quantity supplied at location l of product p in period t under realization ξ t,ξ sl,p ∈ 9 + |l ∈ LB, P ∈ P, t ∈ T, ξ ∈ Ω, shortage of product p at location l in period t under realization ξ t,ξ vl,p ∈ 9 + |l ∈ LB ∪ LR ∪ LS, P ∈ P, t ∈ T, ξ ∈ Ω, stock at location l of product p in period t under realization ξ t,ξ xda,p ∈ 9 + |a ∈ A, p ∈ P, t ∈ T, ξ ∈ Ω, flow in the arc a of product p in period t period under realization ξ t,ξ xia,p ∈ 9 + |a ∈ AI, p ∈ P, t ∈ T, ξ ∈ Ω, inverse flow in the arc a for product p in period t under realization ξ

4287


Optimization of the Petroleum Product Supply Chain under Uncertainty

Recommend Documents