Planning under Demand and Yield Uncertainties in an Oil Supply Chain

ARTICLE pubs.acs.org/IECR

Planning under Demand and Yield Uncertainties in an Oil Supply Chain Kailiang Tong, Yiping Feng, and Gang Rong* Institute of Cyber-System and Control, Department of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China ABSTRACT: In this article, a stochastic programming approach for an optimal refinery planning problem under uncertainties is proposed. The nominal planning model is extended by means of conditional value-at-risk theory. Demand amount uncertainty and product yield fluctuation are taken into account simultaneously. The risk of customer dissatisfaction and inventory violation are considered as constraints according to the decision maker’s risk tolerance. Sample average approximation approach is employed to test the robustness of the model and determine the suitable risk aversion value. A more accurate product yield distribution based upon a Markov chain is introduced and applied in this model. A problem with such endogenous uncertainty is solved using a heuristic iterative algorithm integrating stochastic programming and simulation framework. Also, the scenario number in the stochastic programming model is determined by a statistical analysis, which is a compromise of model accuracy and problem size. At the end of this article, a comprehensive analysis is presented to illustrate the effectiveness of the proposed model.

1. INTRODUCTION Supply chains in the process industry are usually large-scale systems with hundreds of components such as raw material suppliers, production facilities, distribution centers, and customers that perform the following functions: procurement of raw materials, transformation of raw materials into finished products, and distribution of finished products to customers.1 The goal of supply chain management is to achieve high revenues at a high customer satisfaction level and at a low risk level. Supply chain management has three vital components: strategic planning, tactical planning, and scheduling.2 Each of them plays a different role. Strategic planning takes consideration of changing the supply chain structure such as setting up a new site or eliminating an existing facility. The time horizon is usually 1 year or even longer. The objective is to achieve a long-term decision facing the changing market. Tactical planning is usually referred to as “operational planning”. It determines the production profile as well as raw material procurement and final production distribution under many uncertainties within a time horizon from 1 week to 3 months. Scheduling gives the detailed task allocation and unit utilization in the daily production profile with the information given by operational planning. The time horizon of scheduling is usually 1 day or hours. As a typical example of supply chain management, oil supply chains attract many researchers’ attention. Shah et al.3 proposed a comprehensive study about key issues, advances, and opportunities in petroleum refinery operations. Uncertainties in oil supply chains often make the problem more complicated. In the petroleum refining industry, demand, market price, transportation delay, unit capacity, product yield, unit breakdown, and order cancellation can be regarded as uncertain. Among these, the demand amount uncertainty was most discussed by many researchers.46 Only a few considered the product yield uncertainty. Lababidi7 presented a scenariobased two-stage stochastic approach to discuss the uncertainty of key parameters, such as demand, market price, raw material costs, and product yield. Khor et al.8 proposed a two-stage stochastic programming model with fixed recourse via a scenario analysis r 2011 American Chemical Society

with incorporation of risk management for an optimal midterm refinery planning that addresses three factors of uncertainties: prices of crude oil and saleable products, product demands, and product yield. These works assume that product yield obeys normal or uniform distribution. This assumption may lose the dynamic information about product yield. For instance, product yield will be dramatically influenced by operation mode changeover. Yang9 and co-workers assumed that the product yield is only influenced by its last state, which has the Markov property being time-independent and time-homogeneous. A constraint programming approach was employed to keep the inventory over safety stock level within the largest probability. According to this assumption, the distribution of product yield is dependent on the operation mode changeover list, which is called “endogenous uncertainty”.10 Based on the dependence of the stochastic process on the decision, Jonsbraten11 classified uncertainty into two categories: project exogenous uncertainty and project endogenous uncertainty. Problems where the stochastic process is independent of project decisions are said to be exogenous uncertain. For these problems, the scenario tree is fixed. Demand uncertainty in this model is an example of exogenous uncertainty. Problems where the project decisions influence the stochastic process are said to be endogenous uncertain. For problems with such uncertainty, the scenario tree and, hence, the nonanticipativity constraints are decision dependent. Product yield uncertainty obeying Markov chain character is included in this category. There is little literature dealing with problems having endogenous uncertainty. Jonsbraten et al.12 proposed an implicit enumeration algorithm addressing such uncertainty. In Yang’s model,9 product yield probability distribution was generated when all the available mode changeover list was enumerated. Goel and Grossmann10,13 Received: January 27, 2011 Accepted: November 17, 2011 Revised: November 10, 2011 Published: November 17, 2011 814

dx.doi.org/10.1021/ie200194w | Ind. Eng. Chem. Res. 2012, 51, 814–834

Industrial & Engineering Chemistry Research

ARTICLE

dealt with the gas field development problem under uncertainty in size and quality of reserves where decisions on the timing of field drilling yield an immediate resolution of the uncertainty. They first introduced a heuristic approach to deal with such a problem and then developed a more rigorous approach that is a branch-and-bound solution method based on a Lagrangean decomposition scheme relying on relaxing the disjunctions and logic constraints for the nonanticipativity constraints. Solak14 used the sample average approximation method for solving the project portfolio optimization problem that deals with the selection of research and development projects and determination of optimal resource allocations under decision dependent uncertainty where uncertainty was resolved gradually. Balasubramanian and Grossmann15 proposed an approximation strategy based on a series of two-stage stochastic programming approaches integrated with a shrinking horizon. Such an algorithm reduces the computational time significantly. Grossmann and co-workers13,16,17 developed many methods to reduce the computational time of the multistage stochastic programming problem with endogenous uncertainty these years. In this paper, because of the Markov property, the scenario is dependent on the time period t, which is different from the one in the literature aforementioned. The scenario number increases dramatically when the time period increases. We propose a novel iterative algorithm to address such uncertainty where the enumeration number can be greatly reduced. Many sophisticated approaches have been developed to deal with uncertainties, namely two-stage or multistage stochastic programming, parametric programming, chance constraint programming, fuzzy programming, robust programming, and conditional value-at-risk approaches. Each of them has advantages and disadvantages. Two-stage stochastic programming is often scenario based, in which the first-stage variable must be assigned before uncertainties are realized and the second stage is dependent on the realization of uncertainty. The first stage is often referred to as a “here-and-now” decision, while the second stage is referred to as a “wait-and-see” decision.18 Parametric programming is another alternative way dealing with uncertainty explicitly. Based on the sensitive analysis theory, the parametric programming approach provides an optimal solution of the entire uncertainty space.2 Chance constraint programming is a common approach addressing uncertainties. It explicitly takes into account the stochastic manner and reformulates it into a deterministic form using probabilistic theory. Fuzzy programming is often used in the situation where either information is incomplete or uncertainty distribution is not available. Fuzzy member and fuzzy function are the main parts of fuzzy programming. In the robust programming approach, the nominal model is reformulated into a deterministic robust counterpart.19 It guarantees the solution will be feasible even in the worst case. Risk management has received much attention in recent years. Standard stochastic methods usually do not provide any control on the solution variability over different scenarios. In other words, the decision makers are assumed to be risk-neutral.20 Li et al.5 defined a loss function and analyzed the decision maker’s service level, namely, the fill rate and the confidence level. Mulvey et al.21 proposed a robust optimization model to control the mean and the variance of the objective value. Ahmed and Sahinidis22 defined the variability index as a nonnegative continuous variable to manage the trade-off between objective and model variability. Barbaro and Bagajewicz 23 introduced the probabilistic financial risk as a metric of risk for the planning problem under

uncertainty. Eppen et al.24 proposed the downside risk as a risk measure and incorporated it into a two-stage stochastic model. Recently, You and Grossmann20 applied four popular risk measures including variance, variability index, probabilistic financial risk, and downside risk into the global chemical supply chain planning problem and presented a comprehensive comparison. Conditional value-at-risk (CVaR) theory is a newly developed risk management alternative dealing with uncertainties. Compared with the value-at-risk (VaR) theory, CVaR performs a more optimistic view as it is the expectation of value-at-risk. Conditional value-at-risk was first introduced by Rockafellar.25 He applied CVaR theory to a portfolio problem. Both value-atrisk and conditional value-at-risk can be used to constrict the level of risk according to the users’ degree of risk aversion.26 Value-at-risk can be defined as the maximum loss that exceeds a probability of 1 β, whereas conditional value-at-risk is the expected loss given that the loss is greater than the value-at-risk for a given confidence level of β. Compared with VaR, the explicit distribution of the stochastic vector is not needed.25 Carneiro27 applied CVaR theory to the oil supply chain under supply, demand, and price uncertainties. The objective is to maximize the gross profit while keeping risk under a certain value. Verderame and Floudas26 proposed a novel CVaR framework in midterm planning under demand and transportation time uncertainty. Later, they28,29 applied this framework with different uncertainties in other models. Also, they presented a comprehensive comparison between the CVaR results and the robust optimization results. In scenario-based stochastic programming approaches, the discrete representation of the uncertainty is one of the key issues. When the scenario number increases, the problem size increases exponentially and the problem becomes computationally intractable. On the other hand, the model will be inaccurate and the solution will be not robust enough if the scenario number is too small. Recently, many researchers applied a statistical approach to determine the scenario number to achieve a desired confidence level. Shapiro30 proposed a Monte Carlo sampling approach determining the sampling size in a stochastic programming problem. He states that for a problem with 51000 scenarios, a sample size of around 400 can find the optimal solution with probability 95%. You et al.20 applied this approach in a large industrial-size supply chain problem. He used a sampling number of 600 and got the solution within the desired confidence level. Mele et al.31 used another statistical model to determine the scenario number based on the work of Law and Kelton.32 This approach reduced the problem size dramatically. Karuppiah et al.33 proposed a novel heuristic approximation strategy to reduce the scenario number. This method can reduce the number of scenarios while keeping the worst scenario in the reduced set of scenarios. The objective of this work is to apply conditional value-at-risk theory to an oil supply chain planning framework under demand uncertainty and product yield fluctuation simultaneously, which has not been addressed in previous literature. In this article, Gu’s6 model is extended with conditional value-at-risk theory in a twostage stochastic programming framework. In the supply chain planning model, the objective is to minimize the total cost associated with purchasing cost, transportation cost, mode changeover cost, and the penalties related to demand dissatisfaction and inventory violation. Simple scheduling period and operation mode changeover are integrated into a short-term planning model in order to get a more precise solution. Demand uncertainty, which causes customer’s dissatisfaction, and product yield fluctuation, which affects the inventory level of a storage 815

dx.doi.org/10.1021/ie200194w |Ind. Eng. Chem. Res. 2012, 51, 814–834


ARTICLE

site, are taken into account simultaneously. The customer dissatisfaction and inventory violation are captured in a loss function in our work. The associated risks about loss being greater than a certain value at a predefined confidence level are considered as hard constraints in the model. The application of CVaR is based on the work of Verderame and Floudas.26 In their work, loss function about overproduction or underproduction was defined in an aggregation manner. It would lose the detailed information about each product in each period. In this article, each product at a different site in a different period has its own loss function. Every risk is captured in our model. To reduce the problem size caused by the scenario number, several statistical methods are employed to determine the scenario number and a sample average approximation approach is adopted to determine the proper risk aversion value. Just as aforementioned, the assumption about product yield obeying a normal or uniform distribution is unreasonable. We characterized the product yield fluctuation as a Markov chain process in our work. Such uncertainty depending on decision is called “endogenous uncertainty”.11 The problem with such uncertainty cannot be incorporated with CVaR theory directly. As a result, a novel heuristic iterative algorithm integrating approximation of two-stage stochastic programming and simulation framework is introduced to overcome this difficulty. Comprehensive analysis is presented in this work to illustrate the effectiveness of this model and to show how the CVaR parameter works in this framework. The remainder of this article takes the following form. Section 2 presents a detailed description of the oil supply chain problem. A mathematical deterministic model is given in section 3. In section 4, CVaR theory is applied to this model with demand uncertainty and product yield fluctuation. Statistical methods are employed to determine the sampling size and proper risk aversion. A more accurate yield distribution is adopted and a novel algorithm is introduced in section 5. In section 6, comprehensive analysis is covered to illustrate the effectiveness of the proposed model. A brief conclusion is made in section 7.

Figure 1. Network of oil supply chain.6

The inventory should be higher than the safety stock most of the time. Any violation of the inventory will be penalized in the objective function.

3. MODEL FORMULATION OF THE DETERMINISTIC SHORT-TERM PLANNING MODEL The short-term planning model aims at minimizing the total cost. This work is extended with Gu’s short-term planning model6 and Karimi’s work.34 A brief description is given below. The objective function in eq 1 is to minimize the cost associated with purchasing cost, transportation cost, inventory holding cost, changeover cost, and penalty with inventory violation and customer dissatisfaction, respectively. Equation 2 denotes that the production volume of product i equals the processing amount times the yield of this product. Equation 3 indicates that the flow of feedstock i equals the processing amount times the corresponding mixing ratio Xi,m. Equation 4 shows that only one operation mode should be executed in each scheduling period. Equation 5 limits the upper bound and lower bound for the processing amount in each scheduling period. If the operation mode is not executed, the processing amount is forced to be zero. Equations 6 and 7 indicate whether there exists an operation mode changeover. Equations 8 and 9 indicate the relationship between binary variables Zm,t,k and CHm,t,k and their corresponding auxiliary variables ZAm,h and CHAm,h. ZAm,h and CHAm,h indicate whether operational mode m is executed and whether there exists an operation mode changeover in aggregated period h. To reduce the problem size, these auxiliary variables are regarded as continuous variables. Equations 10 and 11 ensure that one operation mode at least continues for tn scheduling periods to guarantee steady production. Equation 12 states the material Δ in eq 13 captures the shortage of balance in storage sites. Dc,i,t product i at site c in period t. It is penalized in the objective function. Equation 14 ensures that products delivered to customers should never be more than their demand. Similar to Δ and IVΔ Dc,i,t s,i,t in eq 15 capture the violation of safe stock level in storage sites, which will be penalized in the objective function as well. Equation 16 provides both upper and lower bounds for material purchased. Equations 17 and 18 present the lower and upper bounds for decision variables.

2. PROBLEM STATEMENT We consider a simplified oil supply chain optimization problem. This supply consists of crude oil suppliers, a jetty tank area where crude oil is unloaded, a crude oil tank which stores crude oil before production, a refinery with an inplant and an outplant as its input and output interfaces, final product tanks where products are stored temporarily, distribution centers where products are stored for demand fulfillment, and customers where the demand comes from. The network of the oil supply chain is shown in Figure 1. A short-term planning model in the oil supply chain which extends Gu’s work6 will be presented. In the shortterm planning model, given predictions of product demand, crude oil supply capacity, market price in each period, production information such as process capacities, run time length limit, etc., the objective is to minimize the total cost with transportation, crude oil procurement, and penalties for customer dissatisfaction and inventory violation. Decisions provide a detailed production profile such as process volume and run time length in each scheduling period of a short-term period. In order to capture some scheduling information in the short-term period, scheduling information such as changeovers is integrated into this model. Uncertainties of demand amount and product yield fluctuation are simultaneously considered. In this model, each storage site has a minimum stock level that is called the “safety stock level”.

Objective function: minimize fS ¼

∑

s ∈ SOS , s0 , i, t

þ

∑ ips, i IVΔs, i, t

s, i, t

816

pricei Fs, s0 , i, t þ

þ

∑

hs, s0 vi ∈ Link, t

∑ dpiDΔc, i, t

c, i, t

þ

fcs, s0 , v, t Fs, s0 , v, t þ

∑ hcs, i IVs, i, t

s, i, t

∑ cc CHm, t, k

ð1Þ

m, t, k



ARTICLE

4. CONDITIONAL VALUE-AT-RISK FRAMEWORK Much attention has been paid to the study of uncertainty in the literature. Demand uncertainty is the most discussed one. On the other hand, product yield uncertainty is rarely studied for lack of stochastic information and its dependence on the operation. This article addresses the demand amount uncertainty and product yield fluctuation simultaneously with conditional value-at-risk theory. This framework is scenario-based where a number of scenarios with different demand amounts and product yield values are created. Before we apply CVaR constraints to the deterministic model, a brief introduction of conditional valueat-risk theory is presented. This work is based on those of Rockafella25,35 and Verderame.26 4.1. CVaR Theory Background. A loss function f(x,y) is defined first to introduce value-at-risk, where x is the decision vector and y is the stochastic vector. y has its own distribution p(y). The value-at-risk problem is defined to determine the minimal value-at-risk α, such that the probability of the loss function being less than α is greater than the confidence level of β. This formulation is shown in eqs 19 and 20.

subject to manufacturing constraints: Fs, s0 , i, v, t ¼

Fs, s0 , i, v, t ¼

∑ Qm, t, kYi, m m, k

"s ∈ Soutplant

ð2Þ

∑ Qm, t, kXi, m

"s0 ∈ Sinplant

ð3Þ

m, k

∑m Zm, t, k ¼ 1

"t, k

ð4Þ

Q L Zm, t, k e Qm, t, k e Q U Zm, t, k

"t, k

ð5Þ

CHm, t, k g Zm, t, k Zm, t, k1

"k g 2

ð6Þ

CHm, t, k g Zm, t, k Zm, t1, K

"k ¼ 1

ð7Þ

ZA m, ðt1ÞNk þ k ¼ Zm, t, k

"m, t, k

CHA m, ðt1ÞNk þ k ¼ CHm, t, k

ð8Þ

"m, t, k

s.t. ð10Þ

Z f ðx, yÞ e α

"1 < h0 e tn, m

ZA m, 1 e ZA m, h0

ð11Þ

∑ Fs , s, i, v, t s∑, v Fs, s , i, v, t 0

00

s0 , v

∑s, v Fs, c, i, v, t

∑s, v Fs, c, i, v, t e Dc, i, t

∑ Fs, s , i, v, t e OUs, i s ,v 0

y

ð15Þ

max½0, f ðx, yÞ α pðyÞ dy

Property 1: ϕβ ðxÞ ¼ min Fβ ðx, αÞ

ð23Þ

α

ð16Þ

Property 2:

lower bound constraints:

αβ ¼ left end point of arg min Fβ ðx, αÞ α

Fs, s0 , i, v, t ,

IV s, i, t ,

IV Δ s, i, t ,

CHm, t, k ,

DΔ c, i, t g 0

ð17Þ

ϕβ ðxÞ ¼ Fβ ½x, αβ ðxÞ

"Æs, iæ ∈ SI, t;

∑i IVs, i, t e Caps ,

ð24Þ

Property 3:

upper bound constraints: IV s, i, t e IV Us, i ,

ð21Þ

A function Fβ(x,α) in eq 22 is introduced to replace the former expectation ϕβ(x) for its solving difficulties. Verderame26 summarized four properties of Fβ(x,α) which are demonstrated by Rockafella.25,35 Because of these properties, it suffices to use Fβ(x,α) instead of ϕβ(x) when applying conditional value-at-risk theory.

"s ∈ SOS , i ∈ I RM , t

0

Z

f ðx, yÞ pðyÞ dy

ð22Þ

ð14Þ

"s ∈ SIV , i, t

f ðx, yÞ g αβ

Fβ ðx, αÞ ¼ α þ ð1 βÞ1

ð13Þ

"c, i, t

L IV Δ s, i, t g IV s, i IV s, i, t

OLs, i e

"c, i, t

ð20Þ

Z

ð12Þ

IV

DΔ c, i, t g Dc, i, t

ϕβ ðxÞ ¼ ð1 βÞ1

00

"s ∈ S , i, t

pðyÞ dy g β

where β is the user-specified confidence level having a value between 0 and 1. The expectation of the loss function being greater than α is expressed in eq 21.

supply chain constraints: IV s, i, t ¼ IV s, i, t1 þ

ð19Þ

α

ð9Þ

"h < h0 e h þ tn, m

CHA m, h e ZA m, h0

αβ ¼ min α

ð25Þ

Property 4: "Æs, iæ ∈ SI, t

ð18Þ

min ϕβ ðxÞ ¼ min Fβ ðx, αÞ x

817

x, α

ð26Þ



ARTICLE

Property 1 (eq 23) indicates that the minimization of Fβ(x,α) in terms of α transforms the special function Fβ(x,α) into the conditional expectation of loss function ϕβ(x). Property 2 (eq 24) states that the β value-at-risk can be ascertained by means of Fβ(x,α). Property 3 (eq 25) conveys the fact that Fβ(x,α) is equal to ϕβ(x) when the generic loss variable α is replaced by the β value-at-risk, [αβ(x)], which can be expressed as a function of x. Property 4 (eq 26) demonstrates that the minimization of ϕβ(x) in term of x is equivalent to the minimization of Fβ(x,α) in terms of x and α. For a more detailed analysis of these four properties and their transition to eq 27, please refer to Rockafellar’s work.25,35 When applying scenario-based approaches, eq 22 can be approximated by eq 27. e Fβ ðx, αÞ ¼ α þ ð1 βÞ1

∑sc psc max½0, f ðx, scÞ α

positive value, then the unsatisfied demand amount at site c is beyond the acceptable threshold in period t for scenario sc. The left-hand side in eq 33 represents the scenario-based conditional value-at-risk for the particular loss function in eq 31, while the right-hand side indicates the loss threshold of decision makers. The loss threshold is associated with original demand, product price, and the users’ risk aversion. θdem indicates the decision maker’s risk aversion and takes a value between 0 and 1. fc,dem i, t, sc ¼ ðDc, i, t, sc

dem dem πdem c, i, t, sc g fc, i, t, sc αc, i, t

ð27Þ

dem 1 αdem Þ c, i, t þ ð1 β

where sc is one of the scenarios and psc is the probability for each scenario sc. The original CVaR theory minimizes the risk F~β(x,α) in the objective function, which usually results in a multiobjective programming problem. One of the methods solving multiobjective programming problem is to treat one objective as the objective function while others are regarded as constraints with suitable tolerant values. As a result, minimizing eq 27 is transformed into e Fβ ðx, αÞ ¼ α þ ð1 βÞ1

where θ is a parameter which indicates the decision maker’s risk aversion. To remove the nonlinearity in F~β(x,α), the nonnegative variable of the adjusted loss function π(sc) is introduced and eq 28 is reformulated in eqs 29 and 30.

θβ ðx, αÞ ¼ α þ ð1 βÞ1

∑sc

ð29Þ psc πðscÞ e θ

"c, i, t, sc

ð32Þ

dem Dc, i, t pricei ∑sc psc πdem c, i, t, sc e θ

ð33Þ

4.3. Yield Uncertainty. The inventory in storage sites is directly influenced by the product yield fluctuation. When product yield fluctuates at the original yield level, the outflow from processing unit changes as well. In this model, we divide the storage sites into three types, namely, raw material tanks (rtank), product tanks (ptank), and distribution centers (DC). Only ptank is integrated into the CVaR framework for product yield uncertainty. It should be noted that there is only one ptank and it is in an aggregation manner. Because of the fluctuation of product yield, both the inventory of ptank and the flow feeding ptank are scenario dependent, whereas as other variables are scenario independent. In order to apply CVaR theory with product yield uncertainty, a loss function combined with the inventory of ptank is introduced. Before applying CVaR theory, some modification should be made. First, product yield should be combined with scheduling period. Equation 2 is reformulated in eq 34.

∑sc psc max½0, f ðx, scÞ α e θ

"sc

ð31Þ

"c, i, t, sc

"c, i, t

ð28Þ

πðscÞ g f ðx, scÞ α

∑s, v Fs, c, i, v, t Þpricei

Fs, s0 , i, v, t, sc ¼

ð30Þ

∑ Qm, t, k Yi, m, t, k, sc

"i, t, sc

s ∈ Soutplant

k, m

ð34Þ Because CVaR constraints should be combined with scenario information, only the inventory of ptank is taken into account. To get the direct relation between the inventory of ptank and product yield, eq 12 should be adjusted.

4.2. Demand Uncertainty. In an oil supply chain, demands are predicted based on historical data. The precision of the prediction may be the key issue to the robustness of the model. The demands become much more unpredictable in the changing environment. A robust model incorporated with CVaR theory that deals with fluctuating demands is adopted here. A loss function that captures the unsatisfied demand is introduced. Unlike the loss function in Verderame’s work,26 the loss function here is in a disaggregation manner, which is defined for each product in each site at each period. This loss function is more reliable and can capture more detailed information. In Verderame’s model, the loss function was aggregated with sites, materials, and time periods. It may lose some information. For instance, there are two products A and B in one site. The loss for product A is 100 and the loss for B is 100. This means that demand for product A is not met and product B is overproduced. This feature cannot be captured in Verderame’s model for its aggregation manner. In his model, the loss for demand is 0. Equation 31 defines the loss function with unsatisfied product i demand dem is the at site c in time period t for each scenario sc. In eq 32, πc,i,t,sc dem nonnegative variable of the adjusted loss function. If πc,i,t,sc takes a

IV s, i, t, sc

∑s Fs , s, i, v, t, sc ∑s Fs, s , i, v, t, sc ¼ IV s, i, t1, sc þ ∑ Qm, t, k Yi, m, t, k, sc ∑ Fs, s , i, v, t, sc m, k s ¼ IV s, i, 0 þ ∑ Qm, t , k Yi, m, t , k, sc ∑ Fs, s , i, v, t , sc m, k, t e t s ,t e t ¼ IV s, i, t1, sc þ

0

00

0

00

00

00

0

0

00

0

00

0

0

ð35Þ

"s ∈ Sptank , i, t, sc

Equation 35 reveals how product yield influences inventory level of ptank. Combined with safety stock level limits, eq 15 can be reformulated as eq 36. IV s, i, 0 þ

∑

m, k, t 0 e t

Qm, t 0 , k Yi, m, t 0 , k, sc

L þ IV Δ s, i, t g IV s, i

∑

s00 , t 0 e t

Fs, s00 , i, v, t0 , sc


ð36Þ

After the preparation for product yield uncertainty in CVaR theory, we get the loss function of inventory violation associated 818



ARTICLE

with yield uncertainty. fs,invi, t, sc

¼ IV Ls, i IV s, i, t, sc ¼ IV Ls, i IV s, i, 0 þ

∑

s00 , t 0 e t

Fs, s00 , i, v, t0 , sc

∑ Qm, t , k Yi, m, t , k, sc m, k, t e t 0

0

0

"s ∈ Sptank , i, t, sc ð37Þ

L inv πinv s, i, t, sc g ðIV s, i IV s, i, t, sc Þ αs, i, t

"s ∈ Sptank , i, t, sc ð38Þ

inv 1 αinv s, i, t þ ð1 β Þ

∑sc psc πinvs, i, t, sc e θinv IVLs, i


ð39Þ

The loss function is defined in eq 37 with inventory violation which is directly caused by product yield fluctuainv tion. In eq 38, π s,i,t,sc is the nonnegative variable of the inv takes a positive value, then adjusted loss function. If π s,i,t,sc the inventory violation of product i at site s is beyond the acceptable threshold in period t for scenario sc. The left-hand side of eq 39 represents the scenario-based conditional valueat-risk for the particular loss function in eq 37, whereas the right-hand side indicates the loss threshold of decision makers. The loss threshold is associated with safety stock level and users’ risk aversion θinv. Note that both yield and demand uncertainties have an effect on the inventory violation. The risk of inventory violation accounts for both uncertainties. Eqs 3739 show that the product yield uncertainty has more impact than the demand. In a real-world situation, the influence of inventory violation from demand uncertainty is so small that it can be ignored. Moreover, in this framework, only the inventory level of ptank and the flow feeding ptank are regarded as second-stage variables. Others are first-stage variables. This means that the first-stage variables should be determined before uncertainty is realized. As a result, we assume that all three kinds of storage sites are independent while coping with uncertainties. In other words, DCs have the ability to cope with demand uncertainty without influencing other storage tanks; ptank has the ability to deal with product yield fluctuation. This assumption is reasonable: the capacity of storage tanks in a plant is much smaller compared with the storage capacity in DCs. The influence from demand uncertainty on ptank is so slight that it can be ignored when compared to product yield fluctuation. In the same manner, the uncertainty of product yield will not impact the inventory of DCs. It is reasonable and the same as the situation in realworld application. From the case study in section 6.2, we can Δ "s ∈ S ptank , i, t in the see that there is no positive IV s,i,t deterministic model or under demand uncertainty. In other words, the demand uncertainty will not influence the ptank. Δ only occurs when product yield is considered. On the IVs,i,t other hand, inventory levels in DCs are regarded as first-stage variables that are decided before realization of uncertainties. They are not influenced by product yield fluctuations. Compared with the nominal deterministic model for shortterm planning, some modifications should be made to adopt

Figure 2. Flowchart of sample average approximation (SAA) approach.

the CVaR framework. IV s, i, t, ðscÞ ¼ IV s, i, t1, ðscÞ þ

∑ Fs , s, i, v, t, ðscÞ 0

s0 , v

∑ Fs, s , i, v, t, ðscÞ s ,v 00

00

DΔ c, i, t, sc g Dc, i, t, sc

∑s, v Fs, c, i, v, t

L IV Δ s, i, t, ðscÞ g IV s, i IV s, i, t, ðscÞ

"s ∈ SIV , i, t, ðscÞ

ð40Þ

"c, i, t, sc

ð41Þ

"s ∈ SIV , i, t, ðscÞ

ð42Þ

Equations 4042 imply that these constraints are adopted for both scenario-based and non-scenario-based constraints. The objective function should be reformulated as well. Decision variables with and without scenarios should be addressed separately. Decision variables with scenarios should be added to the objective function with each probability. fS ¼

∑

s ∈ SOS , s0 , i, t

þ

pricei Fs, s0 , i, t

∑

hs, s0 vi ∈ Link, t, sc

ðpsc Þfcs, s0 , v, t Fs, s0 , v, t, ðscÞ

þ

ðpsc Þhcs, i IV s, i, t, ðscÞ þ ∑ ðpsc Þips, i IV Δ ∑ s, i, t, ðscÞ s, i, t, ðscÞ s, i, t, ðscÞ

þ

ðpsc Þdpi DΔ ∑ c, i, t, ðscÞ þ ∑ cc CHm, t, k m, t, k c, i, t, ðscÞ

ð43Þ

Some additional constraints should be added to keep the model reasonable. Equation 44 guarantees the product amount 819



ARTICLE

Figure 3. Overview of heuristic iterative algorithm.

delivered to the customer will not be greater than the maximum demand. This constraint conflicts with eq 14. As a result, eq 14 will not be included in the CVaR model. Equations 45 and 46 provide a lower bound for the demand satisfactory and inventory level.

∑s, v Fs, c, i, v, t e

max Dc, i, t, sc

"c, i, t

sc

yield

Yi, m, t, k, sc ¼ Yi, m ð1 þ ξi, t, k Þ

"c, i, t

ð45Þ

∑s IVs, i, t, sc g 0:7IVLs, i

"c, i, t, sc

ð46Þ

∑

where n is the number of scenarios and Costsc is the total cost of scenario sc in simulation. Then the confidence level of 1 α is given by " # Zα=2 SðnÞ Zα=2 SðnÞ ð50Þ E½Cost pffiffiffi , E½Cost þ pffiffiffi n n

With the aforementioned information, the short-term planning model in the oil supply chain is comprised of objective function eq 43, constraint eqs 318, eqs 3134, eqs 3742, and eqs 4446. Decision variables IVs,i,t,sc, Fs0 ,s,i,v,t,sc " s ∈ Sptank, i, t, v, sc, are regarded as second-stage variables that make “wait-andsee” decisions, whereas the other variables are regarded as firststage variables that make “here-and-now” decisions. 4.4. Sampling Size Determination. In the scenario-based two-stage stochastic programming approach, the problem size increases exponentially as the number of scenarios increases. A toolarge scenario number makes the problem intractable within a considerable time. A too-small scenario number cannot represent the full scenario space. In this article, two Monte Carlo sampling techniques are adopted to obtain a reasonable sampling size. dem and We defined the uncertainties by eq 47 and 48, where ξc,i,t yield ξi,t,k are random seeds. They may obey a specified distribution. Each scenario can be obtained when generating random seeds. Dc, i, t, sc ¼ Dc, i, t ð1 þ ξdem c, i, t Þ

"c, i, t, sc

ð48Þ

The first approach determining the sample size is based on the work of You and Grossmann.20 The Monte Carlo sampling variance estimator of the result for a stochastic programming problem is given by vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi un u u ðE½Cost Costsc Þ t sc ð49Þ SðnÞ ¼ n1

ð44Þ

∑s, v Fs, c, i, v, t g 0:7Di, c, t

"i, m, t, k, sc

where zα/2is the standard normal deviate such that 1 α/2 satisfies for a standard normal distributed variable z ∼ N(0,1), Pr(z e z1α/2) = 1 α/2. On the other hand, if we are given the sampling estimator S(n) and the desired confidence interval H, the minimum number of scenarios required can be determined by !2 zα=2 SðnÞ ð51Þ N1 ¼ H As a result, we can first generate a small set of scenarios and estimate S(n) by eq 49 with the information of the total cost in the simulation. Then, with the user-specific confidence interval H, we can get the approximate scenario number by eq 51. Another approach is based on the work of Mele et al.31

ð47Þ 820



ARTICLE

Figure 4. Heuristic approach integrated with SAA approach.

Statistically speaking, N1 from approach 1 and N2 from approach 2 should be similar. However, they can be quite different in some cases. As a result, the scenario size can be determined by

The objective of this approach is to obtain a good number of sampling size with a simulation cost relative error of γ and a confidence level of 1 α. The algorithm is as follows: Step 1. Make n0 g 2 replications of simulations and set n = n0. Step 2. Compute the expectation of simulation cost E[Cost] and the confidence level half-length δ(n,α), where rffiffiffiffiffiffiffiffiffiffi S2 ðnÞ δðn, αÞ ¼ tn1, ð1αÞ=2 n

N ¼ maxðN1 , N2 Þ

ð53Þ

4.5. Sample Average Approximation Approach. The key issue in applying scenario-based CVaR framework is the solution feasibility and its optimum when considering the full scenario space. To maintain the computational tractability, only a few scenarios are considered in the planning optimization framework. As a result, the sample average approximation (SAA), which is introduced by Wang and Ahmed36 and developed by Verderame,26 is adopted here. The aim of the SAA method is to determine a suitable risk aversion value under small scenario space. The first step is to generate N scenarios, where N is determined by eq 53. Solve the CVaR planning model with this set of

ð52Þ

Step 3. If δ(n,α)/E[Cost] e γ/(1 γ), the sampling size is N2 = n and stop. Otherwise, increase n and go to step 2. It is noted that, in this article, expectation and variance of simulation cost are calculated instead of expectation and variance of uncertain parameters as it is in Mele’s work.31 Similar work can be found in Law and Kelton’s work.32 821



ARTICLE

Figure 5. Simplified supply chain network of a refinery.6

Table 1. Price of Both Raw Materials and Final Products material (I)

type

price (k$/kbbl)

RM

Table 3. Product Yield for Each Operation Mode final product

R1 R2

I IRM

120 124

mode (m)

G90

G93

D

FO

R3

IRM

116

RM

153

m1 m2

0.31 0.19

0.17 0.29

0.32 0.29

0.14 0.17

MTBE

I

G93

IFP

160

m3

0.23

0.15

0.43

0.13

G90

I

FP

153

m4

0.14

0.21

0.36

0.24

D

I

FP

155

FO

IFP

145

dem being the value λ which indicates the percentage of LHSc,i,t,sc dem larger than RHSc,i,t , which means violation of the risk constraint eq 56. If λ is larger than a certain confidence level χ, this means that risk aversion parameter θdem is too large to get a robust solution. We choose a suitable risk aversion that is the compromise of objective value and the violation of constraint. As will be seen in section 6.2.4, the violation percentage does not monotonically decrease as the risk aversion value decreases. With the information about the objective value and violation percentage, the decision maker can choose a suitable risk aversion value. A similar method is applied to inventory risk management as well. Figure 2 presents an overview of the proposed SAA algorithm.

Table 2. Mix Ratio of Crude Oil for Each Operation Mode crude mode (m)

R1

R2

R3

MTBE

m1

0.65

0.3

0

0.05

m2

0.3

0.5

0.1

0.1

m3

0.3

0.2

0.5

0

m4

0.15

0.25

0.6

0

scenarios. Then generate a larger, independent set of scenarios dem dem and RHSc,i,t with each scenario in the and calculate LHSc,i,t,sc simulation according to eqs 54 and 55. dem 1

dem LHSdem c, i, t, sc ¼ αc, i, t þ ð1 β

∑s, v

Þ

inv 1 inv max½0, IV Ls, i LHSinv s, i, t, sc ¼ αs, i, t þ ð1 β Þ

IV s, i, 0 þ

max½0, ðDc, i, t, sc

Fs, c, i, v, t Þpricei αdem c, i, t

"c, i, t, sc

dem LHSdem c, i, t, sc e RHSc, i, t

"c, i, t

"c, i, t, sc

Fs, s00 , i, v, t 0 , sc

Qm, t , k Yi, m, t , k, sc αinv ∑ s, i, t m, k, t e t 0

0

0


ð54Þ dem RHSdem Dc, i, t pricei c, i, t ¼ θ

∑

s00 , t 0 e t

inv L RHSinv s, i, t, sc ¼ θ IV s, i

ð55Þ

inv LHSinv s, i, t, sc e RHSs, i, t

ð56Þ

Theoretically, eq 56 should not be violated in simulation. However, with some large risk aversion θdem, this is not the case. While simulating a large set of scenarios, there really exist some violations of eq 56. The reason is that the optimization algorithm will result in some extreme solutions with a large risk aversion that is not so robust. As a result, we simulate from large to small risk aversion values and record the information. Get the objective value and

ð57Þ

"s, i, t, sc

ð58Þ

"s, i, t, sc

ð59Þ

5. PRODUCT YIELD UNCERTAINTY In most articles, product yield is regarded as deterministic or obeying a normal or uniform distribution.7,8 This assumption may lose some dynamic characters of production yield. In realworld applications, product yield fluctuates dramatically when the operation model is changing. In this article, we consider the 822



ARTICLE

Table 4. Configuration of Site and Corresponding Inventory site (s)

material (i)

IVU s,i (kbbl)

IVSs,i (kbbl)

initial inventory (kbbl)

holding cost (k$/kbbl)

inventory penalty (k$/kbbl)

JTK

R1

10 000

800

1000

3

230

JTK

R2

10 000

800

1000

3

230

JTK

R3

10 000

800

1000

3

230

JTK

MTBE

4 000

400

600

4

230

MTK

R1

20 000

3000

4500

3

250

MTK

R2

20 000

3000

4500

3

250

MTK

R3

20 000

3000

4500

3

250

MTK PTK

MTBE G90

8 000 5 000

800 500

1200 1000

4 12

300 400

PTK

G93

5 000

500

1000

14

400

PTK

D

5 000

500

1000

12

400

PTK

FO

5 000

500

1000

12

400

DC1

G90

15 000

200

1500

8

300

DC1

G93

15 000

200

1500

9

300

DC1

D

15 000

200

1500

8

300

DC1 DC2

FO G90

15 000 15 000

200 200

1500 1500

8 8

300 300

DC2

G93

15 000

200

1500

9

300

DC2

D

15 000

200

1500

8

300

DC2

FO

15 000

200

1500

8

300

DC3

G90

15 000

200

1500

8

300

DC3

G93

15 000

200

1500

9

300

DC3

D

15 000

200

1500

8

300

DC3

FO

15 000

200

1500

8

300

fluctuation of product yield from a dynamic point of view. Markov chain is applied to describe yield fluctuation. This character is closely related to the changeover of operation modes. This work is an extension of Yang, Gu, and Rong’s work.9 5.1. Markov Chain Property. When applying the Markov yield for product i in scheduling chain property, the range of ξi,t,k period k of period t in eq 48 is a series of discrete values Λi = {λi,1,λi,2,...,λi,Ni}, which is called the Markov state space. If there is no changeover of operation mode and the current yield fluctuation is known, the product yield fluctuation in the next scheduling period is independent of historical data. Thus, the Markov chain can be applied to describe the state transition of yield fluctuation ξi,t,kwithin the same operation mode. We can get

Based on the work of Grinstead,37 the n-step transition prob(n) equals the abth entry of matrix Pni , where Pni is the nth ability pi,a,b power matrix of the state transition matrix Pi. If there is no operation mode changeover, the probability of ξi,t,k can be calculated by the state transition matrix Pi. If there is a changeover of operation mode occurring in scheduling period k of period t, we define the initial probability distribution vector of product i’s yield fluctuation Pbi(1). This initial probability distribution is adopted when it is the first scheduling period in the whole horizon or there is a mode changeover occurring. Pbi ð1Þ ¼ ½pbi, 1 ð1Þ, pbi, 2 ð1Þ, :::, pbi, Ni ð1Þ

Prfξi, t, kþ1 ¼ λi, b jξi, t, k ¼ λi, a , :::, ξi, 1, 1 ¼ λi, x g ¼ Prfξi, t, kþ1 ¼ λi, b jξi, t, k ¼ λi, a g

After a new operation mode has been executed, there exists a state path formed by yield fluctuation until the τth scheduling period when a new changeover occurs. We define the state transition path λi,x1 f λi,x2 f ... f λi,xτ (λi,x1, λi,x2, ..., λi,xτ ∈ Λi) as e; then Q the probability of path e can be calculated as pbi,x1(1) τ1 l=1 pi,xl ,xl +1. Given parameters Pi and Pbi(1), if the operation mode changeover list is known, the product fluctuation ξi,t,k at each time period can be obtained by sampling techniques. Each state of ξi,t,k has a probability associated with its state in the previous scheduling period. Here is the procedure generating yield fluctuation uncertainty for product i. Step 1. Random generate ξi,1,1 with probability Pbi(1). Go to step 2. Step 2. Check whether there exists operation mode changeover in this period. If yes, go to step 3; else, go to step 4.

ð60Þ

where λi,a, λi,b, λi,x ∈ Λi. For simplicity, we define the one-step transition probability as pi, a, b ¼ Prðξi, t, kþ1 ¼ λi, b jξi, t, k ¼ λi, a Þ

ð61Þ

The state transition matrix of product i within the same operation mode can be indicated by λi, 1 ::: λi, Ni 0 1 pi, 1, 1 ::: pi, 1, Ni λi, 1 B C Bl C ⋱ l Ρi ¼ l @ A λi, Ni pi, Ni , 1 3 3 3 pi, Ni , Ni

ð63Þ

ð62Þ

where the abth entry of Pi stands for the probability of product yield fluctuation transitioning from state λi,a to λi,b in one step. 823



ARTICLE

Table 5. Forecasted Demand for Each Short Planning Period

Table 6. Unit Transportation Cost source

short-term period (t) customer product

1

2

3

C1

G90

0

C1

G93

0 800

C1

D

0

0

C1

FO

C2 C2

4

0 160

5

6

7

8

9

10

destination

vehicle

unit cost (k$/kbbl)

CrudeSR1

JTK

ship

4.6

CrudeSR2

JTK

ship

3

0

0

0

0 800

0

0

CrudeSR3

JTK

ship

3.5

0

0

0

0

0

0

0 100

MTBESR

JTK

train

4.8

0

0

0 100

0

0

0 600

JTK

MTK

pipe

0

0

0 180

0

0

0

0

0

0

MTK

inplant

pipe

0

G90

0

0

0

0 170

0 800

0

0

0

outplant

PTK

pipe

0

G93

0

0

0

0

0

0 160

0

0

0

C2

D

0

0 100

0

0

0

0 100

0

0

PTK PTK

DC1 DC2

train train

3 3.7

C2 C3

FO G90

150 0

0 130 0 0

0 0

0 0

0 0

0 0

0 0 0 900

PTK

DC3

train

3.9

DC1

C1

truck

2

C3

G93

0

0 100

0

0

0

0

0

0

0

DC1

C2

truck

3

C3

D

0

0

0

0

0

0

0 190

0

0

DC1

C3

truck

3.1

C3

FO

0

0 140

0

0

0

0

0 130

0

DC1

C4

truck

3

C4

G90

0 700

0

0

0

0

0

0

0

0

DC2

C4

truck

2.5

C4

G93

0 800

0

0

0

0

0

0 400

0

C4

D

0 100

0

0

0

0

0

0 100

0

DC2 DC2

C5 C6

truck truck

2.9 3.2

C4 C5

FO G90

0 0

0 0 120 0 200 0

0 0

0 0

0 0

0 100 0 0

0 0

DC3

C6

truck

2.7

DC3

C7

truck

2.8

C5

G93

400

0

0

0

0

0 130

0

0

0

DC3

C8

truck

2.9

C5

D

0

0

0

0 200

0

0

0

0

0

DC3

C1

truck

3

C5

FO

0

0

0 900

0

0

0

0 900

0

C6

G90

0

0

0

0

0

0 130

0

0 700

C6

G93

0

0 130

0

0

0

0

0 400

C6

D

0

0

0

0

0 140

0

0

0 800

C6 C7

FO G90

0 800

0 0

0 0

0 0

0 0

0 0

0 0

0 900 0 0

0 0

C7

G93

100 190

0

0

0

0

0

0

0

0

C7

D

0

0 100

0

0 110

0

0

0

0

C7

FO

0

0 110

0

0

0

0

0

0 800

C8

G90

0

0

0

0

0

0 170

C8

G93

0

0

0

0

0

0

0 900

C8

D

0 800

0

0

0

0

0

0

0 100

C8

FO

0

0

0

0 120

0

0

0

0 0

0 600 190

0 0

0

0

resolved completely and becomes determined after a certain action has been made. Their scenario is independent of the time period. As a result, they model the problem with separate scenario and time period indices. This is not the case in our paper. Because of the Markov chain property, the product yield is still uncertain even though the operation mode list is known to be fixed. In our model, the scenario index is combined with the time period index. The number of scenarios increases exponentially with the time period. As a result, the nonanticipativity constraints are hard to write as the scenario is related to the time period. To overcome the difficulty of solving such a problem with multistage stochastic programming, we proposed a novel heuristic iterative algorithm integrating approximation two-stage stochastic programming and a simulation framework. Just as mentioned in section 4.1, Fβ(x,α) in eq 22 can be approximated by F~β(x,α) in eq 27, and the accuracy of sampling is the key issue. However, the distribution of product yield is based on the operation mode changeovers. It cannot be obtained before the optimization model is solved. We fix the operation mode list first. Then the yield uncertainty can be regarded as exogenous uncertainty, as it is independent of other decision variables. It can be solved using the two-stage stochastic programming approach. Moreover, the chain is a time-homogeneous Markov chain if the operation mode list is fixed. The probability will converge to a steady state within a long time. We also can use the probability in steady state to generate scenarios as an approximation. In such a way, we can get one optimal value with one operation mode list e. Thus, enumerating all the possible operation mode changeover lists and choosing the best one is one of the methods. However, the number of operation mode lists increase exponentially as the time period increases. This makes enumerating all the lists not realistic. As a result, we proposed a heuristic approach based on enumerating the potential operation mode lists in which the number of enumerations is reduced.

0

0

Step 3. Random generate ξi,t,k with probability Pbi(1). Go to step 5. Step 4. Random generate ξi,t,k with probability in state transition matrix Pi with the state in the previous scheduling period. Go to step 5. Step 5. If it is the end of scheduling period, terminate; else, go to the next scheduling period. Go to step 2. 5.2. Heuristic Iterative Algorithm. Just as aforementioned, the distribution of product yield is dependent on the operation mode and the state in the last period. If the operation mode changes at period t, the distribution of the product yield at a period is determined by its initial probability. Otherwise, the distribution is determined by its last state and the transition matrix. However, the operation modes at each scheduling period are variables to be determined in the optimization model. This is a case of endogenous uncertainty. With endogenous uncertainty, the scenario tree is not fixed and is generated only after certain decisions have been made. Researchers10,13,14,17 usually use a multistage stochastic programming approach solving such problems with endogenous uncertainty. In their research, the uncertainty is 824



ARTICLE

Table 7. Other Model Parameters parameter description

notation

value

max refinery throughput (kbbl/schedule period)

QU

2550

min refinery throughput (kbbl/schedule period)

QL

1620

changeover fee (k$/per time)

cc

30

penalty for demand shortage (k$/kbbl)

dpi

[G90: 306, G93: 320, D: 310, FO: 290]

number of periods in short-term horizon

Nt

10

number of scheduling periods per short time period

Nk

3

number of scheduling periods operation mode continues

tn

5

sites’ capacity (kbbl) max/min available crude per long period (kbbl)

Caps OU s,i

[JTK: 30000, MTK: 60000, PTK:16000, DC1DC3: 50000] [CrudeSR1R1: 100000/0, CrudeSR2R2: 100000/0, CrudeSR3R3:

confidence level for loss function

βdem/βinv

0.99

100000/0, MTBESRMTBE: 50000/0]

Table 8. Sample Size Determination for Approach 1 E(Cost) (109$)

S(n)

α

H (109$)

N1

demand

8.8194

0.175

0.01

0.065

48.32

yield demandyield

8.8102 9.4027

0.346 0.391

0.01 0.01

0.125 0.115

51.06 77.06

model

simulation are similar. Therefore, in the proposed heuristic approach, it is not practical to wait until the operation mode list converges. Because in every sampling the scenario sets are quite different, they may lead different solutions with different operation mode changeover lists. Even their simulation results are similar. Our strategy is to enumerate the potential operation mode lists and compare their simulation results, choosing the best one. There are several termination criterions to guarantee getting a “good” solution. A simulation part is integrated into the algorithm. The total cost in simulation is a criterion to evaluate the decision made by the CVaR planning model. We enumerate the possible operation mode list to get the near-optimal solution. The quality of the solution depends on the quality of sampling scenarios. As a result, a good quality solution can be obtained in the simulation framework. Algorithm 2 Step 1. Set titer = 0 and E(Cost)min = +∞. Generate a set of yield scenarios where both yield product fluctuation ξi,t,k dem and demand uncertainty ξc,i,t obey normal or uniform distribution. Step 2. Solve the CVaR model and get the operation mode changeover list CHm,t,k and objective value Obj. Turn into the simulation step. Simulation Step 2.1. Simulate with the solution and a new larger set of scenarios based on the mode changeover list and Markov chain. Get the expectation of the total cost E(Cost). Simulation Step2.2. If E(Cost) e E(Cost)min , set E(Cost)min = E(Cost). Simulation Step 2.3. Fix the operation mode Zm,t,k. As Zm,t,k are the only binary variables in the CVaR framework, the model becomes a linear programming model. Generate a small set of scenarios based on the Markov chain and Zm,t,k. Resolve this model with newly generated scenarios. Record the objective value Obj. Simulation Step 2.4. Repeat steps 2.12.3 several times and record the smallest E(Cost) in step 2 and its corresponding decision. Go to step 3. Step 3. If the termination criterion is met, terminate; else, go to step 4. Step 4. Set titer = titer + 1. Using the information of CHm,t,k, Pi, and Pbi(1) and the theory proposed in section 5.1, dem generate a new set of scenarios for ξyield i,t,k , where ξc,i,t remain the same. Unfix the operation mode Zm,t,k. Go to step 2.

The idea of the proposed algorithm is to utilize the approximation of two-stage stochastic programming instead of a multistage one. Generate a set of yield scenarios based on the operation mode list which is from the last iteration and resolve the problem. Repeat this procedure until the operation mode list remains unchanged. We first introduce a simplified heuristic iterative algorithm to find the possible changeover list, and then it is integrated into a simulation framework to get a more reliable solution. Algorithm 1 Step 1. Set titer = 0. Generate a set of scenarios where both yield and demand unceryield product fluctuation ξi,t,k dem tainty ξc,i,t obey normal or uniform distribution. Step 2. Solve the CVaR model and get the operation mode changeover list CHm,t,k. Step 3. If the termination criterion is met, terminate; else, go to step 4. Step 4. Set titer = titer + 1. Using the information of CHm,t,k, Pi, and Pbi(1) and the theory proposed in section 5.1, dem remain the generate a new set of scenarios, where ξc,i,t same. Go to step 2. Only when both criterions are met can this algorithm stop: iter The percentage in the number of Zm,t,k = Ziter1 m,t,k is larger than the specified confidence level χ1. The deviation of the objective function |(objiter objiter1)/ objiter| is less than the specified confidence level χ2, In practice, the first termination criterion is hard to meet, as Zm,t,k is easily different with different scenario sampling sets, although the objective value changes little. Therefore, this algorithm should be improved. As in the scenario-based stochastic programming approach, the sampling quality of the scenario space is the most important thing. Due to limited computation capacity, we can only use a small set of scenarios in the stochastic programming. As a result, the optimization result only gets a “good” solution. The operation mode list may not be fixed with different sampling scenarios: even both the optimal value and the expected total cost in 825



ARTICLE

Table 9. Sample Size Determination for Approach 2 α

λ

λ/(1 λ)

demand

0.01

0.0074

0.0075

0.0083 (20), 0.0085 (25), 0.0080 (30), 0.0075 (35)

35

yield

0.01

0.0142

0.0144

0.0287 (20), 0.0236 (25), 0.0198 (30), 0.0171 (35),

50

demandyield

0.01

0.0122

0.0124

0.0291 (20), 0.0243 (25), 0.0208 (30), 0.0180 (35), 0.0164 (40),

model

iteration with N

N2

0.0156 (40), 0.0155 (45), 0.0147 (50), 0.0131 (60) 80

0.0162 (45), 0.0152 (50), 0.0137 (60), 0.0115 (80)

two mode lists is the time when the operation mode changes. In order to get a low-cost solution, the number of operation mode changeovers will not be so large, which means that the number of potential “good” operation mode lists is finite and is much smaller than the total operation mode lists. This proves that we can get a near-optimal solution within several steps. Termination criterion 1 indicates the algorithm will stop when the operation mode list converges, although it proves hard to converge. Termination criterions 2 and 3 are the criterions for expected cost convergence. If the expected cost in the simulation does not improve for a predefined number of iterations, the algorithm will terminate. In Algorithm 2, most of the computational part is the simulation part in step 2. As a result, some useless simulation should be abandoned. For example, if the objective value from step 2.3 is much higher than the minimal obj, it is unlikely to get a smaller E(Cost) than the minimal E(Cost). The simulation for such a decision is not needed. In certain iterations, if the operation mode changeover list is the same with the one got from a former iteration, at most one simulation is needed. It can be easily seen, compared with the uniform or normal distribution representation for product yield uncertainty, that the heuristic approach for Markov chain representation can achieve better simulation results, although it needs a little more computational time. It is worth doing so, as it captures the dynamic information of product yield and leads to less expected total cost. 5.3. Integrating Heuristic Iterative Algorithm with SAA approach. The SAA approach presented is to determine the suitable risk aversion value. If the product yield obeys the Markov chain process, the SAA approach should be modified. It is integrated with the heuristic approach dealing such endogenous uncertainties. With each risk aversion, the heuristic iterative algorithm is executed and the violation of CVaR constraints is recorded. We record the violation percentage, objective value, and simulated expected cost with each risk aversion, and then choose a suitable risk aversion value comparing the results. Figure 4 shows a flowchart of a heuristic iterative algorithm integrated with the SAA approach.

Table 10. Model Information for Four Cases scenario binary continuous constraint model

number variable

variable

deterministic demand

0 50

240 240

1 968 50 439

yield

50

240

12 311

demandyield

80

240

95 351

time

solution

number (CPU s) (109 $) 3 896 52 446

2.3 54

8.2263 8.6823

23 285

37

8.3510

111 815

420

8.9603

Figure 6. Operation modes in all 30 scheduling periods for four models.

The termination criterions are modified from the ones in Algorithm 1. If any one of the criterions is met, the algorithm will terminate. Termination Criterion 1. The percentage in the number of iter = Ziter1 Zm,t,k m,t,k is larger than the specified confidence level χ1. Termination Criterion 2. The deviation of total cost |(E(Cost)min E(Cost))/E(Cost)| is less than the specified confidence level χ2. Termination Criterion 3. It reaches the maximal iteration number. The improved algorithm is shown in Figure 3. The block with dashed lines on the right is the simulation part. The one on the left is the termination criterion. Although no rigorous theoretical analysis is possible, the quality of the proposed algorithm can be illustrated in a case study. From the case study in section 6.3, we find that the objective value may vary from different samples of scenarios. However, the simulation result changes little. Another interesting thing is that, for some changeover lists, even if they are different, they have similar simulation costs. This proves that we can get near-optimal solutions with different operation mode lists. The difference between

6. CASE STUDY 6.1. CVaR Model Framework. In order to demonstrate the effectiveness of the proposed framework, a moderate oil supply chain is presented. The motivated case is inherited from Gu’s work. Hence, the supply chain structure is the same as the case in ref 6, whereas some data are adjusted. The refinery processes three crude oils and the raw material methyl tert-butyl ether (MTBE) to make four products. The supply chain network is shown as Figure 5, where R1R3 are three crude oils purchased from three crude sources CrudeSR1 CrudeSR3; JTK, MTK, and PTK respectively denote tank areas 826



ARTICLE

Table 11. Simulation Results of Different Solutions under Different Uncertainties uncertainties demand decisions

yield

demand and yield

case 1

case 2

case 1

case 3

case 1

case 2

case 3

case 4

total cost (10 $)

8.825

8.707

8.796

8.360

9.395

9.392

8.977

8.926

purchase cost (109$)

5.612

6.258

5.612

5.744

5.612

6.258

5.744

6.323

changeover cost (k$)

90

90

90

90

90

90

90

90

transportation cost (109$)

0.494

0.551

0.494

0.499

0.494

0.551

0.499

0.551

9

holding cost (109$)

1.201

1.244

1.201

1.253

1.201

1.245

1.253

1.287

inventory penalty (109$)

0.209

0.134

0.778

0.197

0.778

0.818

0.197

0.195

demand penalty (109$)

1.310

0.520

0.711

0.667

1.310

0.520

1.285

0.570

Figure 7. Histogram results for different decisions under different uncertainties.

of jetty, refinery feedstock, and final products; DC1DC3 are three distribution centers; and C1C8 stand for eight customers or sales regions. The refinery produces four final products 90# gasoline (G90), 93# gasoline (G93), diesel (D), and fuel oil (FO). There are four operation modes considered in this motivated case. Each operation mode has a different pipeline connection between units in the refinery. From a holistic view of the refinery, each operation mode has its own feedstock consumption ratio and product yield. The planning horizon is divided into 10 short planning periods, whereas each short planning period consists of three scheduling periods. Each operation mode lasts in the whole scheduling period. In other words, changeover only occurs at the beginning of each scheduling period. To guarantee a

steady production, each operation mode continues for at least five scheduling periods. The price of both raw materials and finished products are shown in Table 1. Information about the operation mode for the mixing ratio and product yield is listed in Tables 2 and 3. Configurations of storage sites and inventory are shown in Table 4. Table 5 presents the forecasted demand in the whole short planning horizon. Unit transportation cost is shown in Table 6. Other configurations including CVaR parameters are presented in Table 7. 6.2. Comparative Result. 6.2.1. Case Description. In this subsection, we present four cases to illustrate the effectiveness of the proposed framework and to show the influence of different risk aversion values. The model in case 1 is deterministic, in 827



ARTICLE

Figure 10. Costs with different demand risk aversions in case 4. Lines in red, blue, green, and black denote inventory risk aversion θinv = 0, 0.08, 0.12, and 0.2, respectively.

Figure 8. Costs with different demand risk aversions in case 2.

Figure 9. Costs with different inventory risk aversions in case 3. Figure 11. Costs with different inventory risk aversions in case 4. Lines in red, blue, green, and black denote demand risk aversion θdem = 0.22, 0.25, 0.30, and 0.35, respectively.

which demand and yield are treated as fixed values. In case 2, demand is generated by eq 47, and ξdem c,i,t obeys a normal dis2 tribution, ξdem c,i,t ∼ N(0,0.1 ), whereas yield is still treated as deteryield ministic. In case 3, product yield is generated by eq 48 and ξi,t,k yield obeys a uniform distribution, ξi,t,k ∼ U( 0.1,0.1), whereas demand is deterministic. In case 4, demand and product yield are yield 2 both stochastic, where ξdem c,i,t ∼ N(0,0.1 ) and ξi,t,k ∼ U( 0.1,0.1). These four cases are referred as the deterministic model, demand model, yield model, and demandyield model, respectively. In the models where demand uncertainty is considered, the risk aversion of customer dissatisfaction θdem is set to 0.25. Meanwhile, in the models where yield uncertainty is considered, the risk aversion of inventory violation θinv is set to 0.1. On the basis of the theory presented in section 4.4, we evaluate the scenario sampling size. In approach 1, α is set at 0.01. Values of the confidence horizon H in the demand, yield, and demandyield cases are set at $0.065 109, $0.125 109, and $0.015 109, respectively. We got sample sizes of 48, 51, and 77 for the three cases. Detailed information is listed in Table 8.

In approach 2, α is the same as in approach 1. To keep pace with approach 1, λ is set according to H. We iterate the sample size with a step of 5. Approximated sample sizes of 35, 50, and 80 are determined. Detailed information is listed in Table 9. Both approaches get similar results. As a result, the sample sizes for these cases are set to be 48, 51, and 80, respectively. For simplicity, the sample sizes are set to be 50, 50, and 80, respectively. All the models are all implemented in CPLEX 11.0 on an Intel Core DUO 2.80 GHz CPU and 2 GB RAM. The optimal tolerances in all cases are set to 0. Table 10 lists the problem size, solving time, and objective value in four models. From Table 10, we can conclude that, when uncertainties are taken into account, the objective value increases as well. The problem size increases exponentially with the increasing scenario number. Figure 6 presents the operation modes in 30 scheduling periods 828



ARTICLE

Figure 14. SAA results for different demand risk aversions. Lines in red, blue, green, and black denote inventory risk aversions 0, 0.08, 0.12, and 0.2, respectively.

Figure 12. SAA result for different demand risk aversions in case 2.

Figure 13. SAA results for different inventory risk aversions in case 3. Figure 15. SAA results for different inventory risk aversions. Lines in red, blue, green, and black denote demand risk aversions 0.21, 0.25, 0.31, and 0.35, respectively.

for four models respectively. The points in pink, red, blue, and black represent four operation modes m1, m2, m3, and m4 executed in the corresponding scheduling periods. Operation mode m1 is not executed in whole scheduling periods with all four models. The reason is that the refinery consumes too much R1 and produces too much G90 in mode m1. That will cause a shortage of R1 in the storage tanks and an overproduction of G90. We also can find some interesting features. Compared with operation modes in the deterministic model, the ones in the demand model change a lot, whereas the operation modes in the yield model remain almost unchanged. A similar conclusion can be draw when comparing these two models with the demand yield model. 6.2.2. Stochastic Results. In this subsection, we test the solutions obtained from four models under different uncertainty scenarios. All the simulation results are presented in Table 11 and Figure 7. Table 11 presents the expected costs with different decisions under demand and yield uncertainties for 1000 scenarios. The cost consists of the purchase cost, changeover cost, transportation cost, holding cost, and penalties for inventory

violation and demand shortage. In the demand uncertainty scenario, we test the solution from the deterministic model and the demand model. The total cost in solution 1 is $8.825 109, which is a little higher than the one in solution 2. The main reason is that the solution obtained from the demand model can lead to less of a demand penalty compared with the solution obtained from the deterministic model. Also, purchasing more raw materials can hedge the risk of demand dissatisfaction. From the histogram in Figure 7, decision 2 leads to both less expectation and less variance of total cost. In the yield uncertainty scenario, we test the solution obtained from the deterministic model and the yield model. The cost with decision 1 is 5% higher than the one with decision 3. The main reason is that the decisions obtained from the inventory model can reduce the inventory penalty in a yield uncertain world. It should be noted that constraint 46 is relaxed in the simulation. Otherwise, most simulations with decision 1 will be infeasible as the inventory in 829



ARTICLE

Table 12. Parameters of Product Yield Fluctuation parameter

value

Markov state space: Λi for ξi,t,k

½ 0:2

initial probability distribution when operation mode is changed: Pbi(1)

½ 0:7 0:25 0:05 0

0:1

0:05 0 0:05 0:1 0 0

state transition matrix when operation modes executed continuously: Pi 20:2 0:2 0 6 0:1 6 0 6 0:05 6 0 6 0 6 60 0:05 6 40 0:1 0

0:1 0:05 0:5 0:3 0:15 0:2 0:4 0:25 0:2 0:3 0:3 0:1 0:25 0:3 0:05 0:15 0:3 0:05 0:1 0:25

0

0:05 0:13 0:05 0 7 0:1 0:05 7 7 0:15 0:05 7 7 0:25 0:1 7 7 0:3 0:2 7 5 0:4 0:2

less than 0.16. Figure 8 records the trajectories of corresponding costs with different risk aversions. All the costs increase as the demand risk aversion θdem decreases, but the demand penalty decreases. This is easy to explain. When θdem decreases, the demand shortage decreases directly; thus, the demand penalty decreases. Meanwhile, other operational costs may increase to guarantee low customer dissatisfaction. Inventory penalty occurs when θdem is less than 0.34. This illustrates that, with high demand risk aversion, the inventory of storage sites will be greater than the safety stock level. From Figure 8, we also can see that, when θdem decreases to 0.2, the inventory penalty increases significantly and the holding cost decreases. In case 3 where the yield uncertainty is considered, θinv is iterated from 0.3 to 0 with a step of 0.01. Figure 9 presents the costs with different inventory risk aversions. Costs increase as θinv decreases, except the penalties of inventory violation and demand shortage. However, all the changes are quite slight compared to corresponding costs. When θinv is larger than 0.28, the demand penalty increases as θinv decreases. In case 4 where the demand uncertainty and yield fluctuation are taken into account simultaneously, θdem is iterated from 0.35 to 0.15 with a step of 0.01; meanwhile, θinv is iterated from 0.2 to 0 with a step of 0.02. Iteration terminates when θdem = 0.21. Solution becomes infeasible when θdem is less than 0.21. Figures 10 and 11 present the costs with different demands and inventory risk aversions. The trends are similar to the ones in cases 2 and case 3. It can be concluded again that the objective value changes significantly as θdem changes, but it changes little as θinv changes. 6.2.4. Sample Average Approximation Approach. To determine a proper risk aversion value, the evaluation of the objective value is far beyond enough. In a two-stage stochastic programming model, model robustness is also a key issue. A solution to an optimization model is defined as solution robust if it remains “close” to optimal for all scenarios of the input data and model robust if it remains “almost” feasible for all data scenarios.21 In this subsection, an SAA approach described in section 4.5 is utilized to determine proper risk aversions. In case 2 where demand uncertainty is considered, θdem is iterated from 0.4 to 0.16 with a step of 0.01. In each iteration, a small set of 50 scenarios are generated and case 2 is solved. With such a decision and a new independent set of 1000 scenarios, simulations are done and the violation count of constraint eq 56 is recorded. Figure 12 presents the SAA results. The total demand violation percentage decreases from 2.2 to

Figure 16. Operation modes in heuristic iterative algorithm in case 5.

PTK changes dramatically in a yield fluctuating environment. Even when constraint 46 is relaxed, there are still some infeasible solutions in the simulation. Note that the inventory of PTK will be 0 in some periods with unwise decisions. It is easy to explain the infeasibility of the simulation. As shown in Figure 7, most simulation total costs lie between $8.3 109 and $8.4 109 with decision 3. On the other hand, with decision 1, the total cost ranges from $8.3 109 to $9.5 109, which is mainly caused by the inventory penalty. A similar result can be seen from the simulation under both demand and yield uncertainties. Decision from case 4 leads to the least inventory penalty and demand penalty, which leads to the least total cost, although the purchasing cost with decision 4 is a bit high. Clear results can be seen from Figure 7. From this subsection, a simple conclusion can be drawn that the proposed model can hedge the risk caused by the demand and yield uncertainties significantly. 6.2.3. Iteration with Different Risk Aversion. Just as noted in subsection 4.1, the optimization result will be influenced by different risk aversion values. In this subsection, we analyze the relation between risk aversion and objective value. In case 2 where only the demand uncertainty is considered, we iterate θdem from 0.5 to 0 with a step of 0.01. Iteration terminates when θdem equals 0.16. Solution becomes infeasible when θdem is 830



ARTICLE

Table 13. Information for Each Iteration in Case 5 information

iteration 1

iteration 2

iteration 3

iteration 4

iteration 5

iteration 6

minimum simulation cost E(Cost) (109 $)

9.1344

8.6147

8.5589

8.5408

8.5253

8.5193

corresponding objective value (109$)

8.3332

8.5631

8.5256

8.5298

8.5165

8.5003

infeasible solution number

43

1

0

0

0

0

of Figure 13. Similar to case 2, the figure about violation is not smooth. We also take two individual examples. Violations of G93 at PTK in period 6 and FO at PTK in period 10 are taken into account. Finally, we choose θinv = 0.09 as the proper inventory risk aversion. Similar to case 2 and case 3, in case 4 where the demand uncertainty and yield fluctuation are simultaneously taken into account, θdem is iterated from 0.35 to 0.21 with a step of 0.02. Meanwhile, θinv is iterated from 0.2 to 0 with a step of 0.02. As no feasible solution will be obtained when θdem is less than 0.21, iterations are not executed with smaller θdem. Figures 14 and 15 present the SAA results for case 4. Similar conclusions can be drawn compared with SAA in case 2 and case 3. Meanwhile, some new interesting features are also presented. As the demand risk aversion θdem decreases, both the violation percentage of inventory and demand decrease notably. So does the infeasible percentage. Compared with θdem, change of θinv makes little effect. It should be noted in the last subfigure in Figure 14that the violation of G90 at C6 in period 10 with θdem = 0.29 and θinv = 0 is relatively high compared to it with other θdem. The reason for this phenomenon is understandable. With low inventory risk aversion θinv = 0, the decision will change notably when the demand risk aversion changes a little. When the demand risk aversion shifts from 0.31 to 0.29, in order to get a minimum cost, the decision maker scarifies this product at this site for others. As a result, the violation increases a lot. Similar situations can be seen with demand risk aversion θ dem = 0.27 in the same subfigure. From Figure 15 we can see that both total cost and violation percentage are influenced little with different inventory risk aversion θinv. To trade off between low cost and model robustness, we choose θdem = 0.25 and θinv = 0.1 as proper risk aversions. 6.3. Accurate Product Yield Distribution. Just as mentioned in section 5, it is not quite reasonable for product yield fluctuation to obey normal or uniform distribution. In this subsection, yield fluctuation based on the Markov chain is described and the effectiveness of the proposed algorithm is demonstrated. To illustrate the efficiency of the proposed algorithm, problems with different sizes are considered. Case 5 is modified from case 3 where only product yield uncertainty is considered. The only difference between case 5 and case 3 is the distribution representation of product yield. For simplicity, we assume that each product has the same initial probability distribution and state transition matrix. Table 12 presents the parameters of product yield fluctuation based on the Markov chain. After the operation mode has been changed from one to another, it needs time to reach the steady condition. As a result, the product yield is likely to decrease for the first period after the operation mode changeyield with a value of 0.2 at a large over. That is why we assign ξi,t,k probability after the operation mode changeover. As mentioned in subsection 5.3, the heuristic iterative algorithm should iterate with each inventory risk aversion value in order to determine a suitable risk aversion value with fewer CVaR constraint violations. However, to simplify the problem, we use

Figure 17. Operation modes in heuristic iterative algorithm in case 6.

1.2% and then increases to 1.6% as θdem decreases from 0.4 to 0.25 and to 0.16. On the other hand, the maximum demand violation individual decreases from 13 to 5.8% and remains unchanged as θdem decreases from 0.4 to 0.26 and to 0.16. This is easy to understand. Constraint eq 56 reduces the risk of demand violation of each product in each customer zone at each period as θ dem decreases. To achieve a minimum cost, the inventory for other products in another site at different periods may reach the lower bound of constraint eq 56 and easily be violated. We take two individual examples. Violation of product G90 at C8 in period 9 decreases from 2.2 to 0.4% and remains unchanged as θdem decreases from 0.4 to 0.26 and to 0.16. Most situations are of this type. Violation of product G90 at C6 in period 6 increases from 0 to 2.8% as θdem decreases from 0.26 to 0.2, and remains unchanged with other θdem. Few situations are of this type. The reason is similar to the one aforementioned. The lines are not smooth enough. This denotes a slight change of θdem can lead to a large change of decision. From Figure 12, we can see that θdem = 0.27 may be the proper demand risk aversion as it yields low customer dissatisfaction and total cost. In case 3 where yield uncertainty is considered, inventory risk aversion θinv is iterated from 0.3 to 0 with a step of 0.01. In each iteration, a small set of 50 scenarios are generated and case 3 is solved. With a new independent large set of 1000 scenarios and the solution from case 3, simulations are done and numbers for violation of constraint eq 59 are recorded. Figure 13 presents the SAA results. Both the infeasible solution number and violation numbers decease as inventory risk aversion θinv decreases, but total cost increases. The increasing speed for the total cost is slow. Therefore, the total cost is acceptable even when θinv drops to 0. In this kind of simulation, constraint eq 46 is not relaxed, as with low inventory risk aversion most solutions will be feasible even if yield changes. This can be seen from the second subfigure 831



ARTICLE

Table 14. Information for Each Iteration in Case 6 information

iteration 1

iteration 2

iteration 3

iteration 4

iteration 5

iteration 6

minimum simulation cost E(Cost) (109$)

9.3162

8.9953

9.0089

9.0010

9.0127

9.0075

corresponding objective value (109$)

8.8062

8.9523

8.9653

8.9562

8.9673

8.9695

infeasible solution number

18

11

11

11

11

11

the inventory risk aversion value θinv of 0.1, which is determined from the SAA approach with the product yield with uniform distribution. The model contains 240 binary variables, 12 311 continuous variables, and 23 285 constraints for the first iteration with uniform distribution and is solved in 42 CPU s. In order to capture the cost information, the constraint eq 46 is ignored in the simulation framework. Here we give a simple description of the procedure of the heuristic iterative algorithm. In the first iteration, the product yield scenarios are generated according to the uniform distribution. After the problem has been solved, a large set of 1000 scenarios based on the Markov chain is generated and simulation is made to get the average total cost E(Cost) = $9.1344 109. This means, with the decision made from uniform distribution in iteration 1, the expected total cost is $9.1344 109 under the Markov chain distribution environment. Update E(Cost)min = $9.1344 109. In the next iteration, a small set of scenarios based on the Markov chain and the operation mode changeover list obtained from the last iteration is generated. The CVaR model is applied with this new set of scenarios. Then go to the simulation part. First, generate a large set of scenarios based on the current mode list. This is the simulation environment for the simulation part. Simulate with the decision and yield uncertainties and get the expected total cost. Fix the operation mode in the simulation part. Generate another small set of scenarios based on this mode list. Solve it and simulate in the same environment. Repeat this procedure several times (three times in this case) and get the smallest expected cost E(Cost) = $8.6147 109 and corresponding objective value E(Cost) = $8.5631 109. Compared with E(Cost)min, update E(Cost)min = $8.6147 109. Go to the next iteration; unfix the operation mode list. Get a new solution with the new mode list. Repeat the procedure described above until a maximum iteration number of 6 is reached. Figure 16 presents the operation mode in each iteration. Table 13 presents the minimal simulation cost, corresponding objective value, and infeasible count in each iteration. It should be noted that, in iteration 5, there is only one simulation needed, because the operation mode list is the same as the one in iteration 4. After six iterations, the best simulation result is in the sixth iteration where E(Cost) = $8.5193 109 and the corresponding objective value is $8.5003 109. In order to illustrate that the proposed algorithm is not only effective for small size problems, a more complex problem is considered. In case 6, both demand and yield uncertainties are considered. Each planning period consists of six scheduling periods instead of three as in previous cases. The distribution of product yield is based on the Markov chain and is shown in Table 12. Maximum and minimum process volumes are assigned as 1275 and 810, respectively. The minimum time that one operation mode should continue is assigned as nine scheduling yield in eq 48 is based on periods. The procedure of generating ξi,t,k the Markov chain whose procedure is described in section 5.1. The scenario number is set at 50. Both the demand and inventory risk aversion values θdem and θinv are set at 0.1. Other configurations and parameters are the same as for case 4. It should be

noted that the demand scenarios are remain unchanged in the heuristic iterative algorithm for simplicity. Case 6 contains 480 binary variables, 60 701 continuous variables, and 73 401 constraints and is solved in 356 CPU s for iteration 1. The operation mode list and other information for each iteration are presented in Figure 17 and Table 14. After six iterations, the best simulation result is in the second iteration where E(Cost) = $8.9953 109 and the corresponding objective value is $8.9523 109. There are 11 infeasible solutions in the simulation of the last five iterations. This results from the sampling of demand uncertainty. Demand uncertainty with normal distribution may violate constraint eq 45 at a low probability. It should be noted that in the simulation part for each iteration, if the objective value from the latest optimization is much larger than the minimum objective value, it is unlikely to get the smaller simulation cost than the one obtained from the minimum objective value. The simulation for this decision can be abandoned. For example, in iteration 4, the best objective value since then is $8.9523 109. In this iteration, if the objective value obtained from step 2.2 is much larger than $8.9523 109 (larger than $8.9650 109 in this example), there is no need to simulate with this decision. This saves much computational time. From these two cases, we find that with a fixed operation mode changeover list the objective value may vary from different samples of scenarios. However, the simulation result changes little. Another interesting thing is that, for some changeover lists, even if they are different, they have similar simulation costs. This shows that we can get a near-optimal solution with different operation mode changeover lists. It also proves that we can get a near-optimal solution within several steps. From the results, we can easily see the improvement from using the Markov chain instead of uniform distribution (iteration 1 in case 5 and case 6). It reduces much cost, because the Markov chain captures the dynamic characters of product yield. Product loss is captured with the Markov chain in the proposed algorithm. Also, from Tables 13 and 14, we can conclude that not only the cost in simulation but also the infeasible solution number is reduced.

7. CONCLUSION In this article, a two-stage stochastic programming approach based on the novel application of conditional value-at-risk to the problem of demand amount uncertainty and product yield fluctuation has been presented. A statistical analysis of sampling number is done to trade off between problem size and model accuracy. The sample average approximation approach is employed to determine a suitable risk aversion. A more reasonable distribution of yield uncertainty based on the Markov chain is introduced and incorporated into the CVaR framework. This model is solved using a novel heuristic iterative algorithm integrating simulation framework. Comprehensive analysis is presented to illustrate the effectiveness of the proposed model and algorithm. Results show the product yield with the Markov 832



ARTICLE

Fs,s0 ,i,v,t = flow quantity of product i from site s to s0 via vehicle v in time period t, where Æs,s0 ,i,væ ∈ Tran IVs,i,t = inventory level for product i ∈ I at the end of time period t at site s = IVΔ s,i,t deviation below safety stock level for product i ∈ I at site s ∈ SIV in time period t Qm,t,k = processing amount of total feedstock to the plant in schedule period k in short period t under mode m Zm,t,k = binary variable, denotes whether mode m is used in scheduling period k of short-term t ZAm,h = auxiliary variable, denotes Zm,t,k in whole horizon, where h = (t 1)Nk + k CHm,t,k = binary variable, denotes whether the production process is changed from other modes to mode m in scheduling period k of short-term t CHAm,h = auxiliary variable, denotes CHm,t,k in whole horizon, where h = (t 1)Nk + k

chain can achieve better performance compared with uniform distribution representation. As discussed, the proposed CVaR framework can efficiently reduce the total cost in an oil supply chain and the risk of customer dissatisfaction and inventory violation simultaneously under demand and yield uncertainties. Also, this framework can be extended to different industry areas under both exogenous and endogenous uncertainties.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT The authors acknowledge financial support from the National High Technology Research and Development Program of China (No. 2009AA04ZX1485932) and the National Basic Research Program of China (No. 2012CB720500).

Parameters

pricei = price of raw material and finished product i hcs,i = inventory cost for holding a unit of product i ∈ I at site s in a time interval of short-term time period dpi = penalty for demand shortage below order of final product i ∈ IFP ips,i = penalty for inventory shortage below safety stock of material i ∈ I at site s cc = changeover cost for production mode change one time Dc,i,t = predicted demand for finished product i ∈ IFP for customer c due at the end of time period t Caps = capacity of site s S IVU s,i, IVs,i = upper bound and safety stock bound of inventory level for product i ∈ I at site s ∈ SIV L U Os,i, Os,i = lower bound and upper bound of raw material i ∈ IRM an oil or MTBE supplier s ∈ SOS could supply in a long period QL, QU = lower bound and upper bound of the processing amount of total feedstock to the plant in a scheduling period Xi,m = consumption ratio of raw material i ∈ IRM under operation mode m Yi,m, Yi,m,t,k = yield of product i ∈ IFP under operation mode m (in scheduling period k of period t) tn = minimum scheduling periods that one operation mode must continue fcs,s0 ,v,t = unit transfer cost from site s to s0 with vehicle v in period t psc = probability of scenario sc dem = random seed for demand uncertainty of product i in ξc,i,t customer c in period t yield = random seed for product yield uncertainty of product i in ξi,t,k scheduling period k of period t

’ NOMENCLATURE Sets

I {i} = set of materials and products IRM ⊂ I = set of raw materials IFP ⊂ I = set of finished products S {s} = set of sites, including suppliers, jetty tank areas, crude oil tanks, manufacturing plants, product tanks, distribution centers, customers SOS ⊂ S = set of crude oil and MTBE suppliers C {c} ⊂ S = set of customs; they are considered to be one kind of site for the sake of convenience DC ⊂ S = set of distribution centers Sinplant ⊂ S = label of the plant, considering the consumption of raw materials Soutplant ⊂ S = label of the plant, considering the output of final products SIV S\(SOS ∪ Sinplant ∪ Soutplant ∪ C) = set of sites which are the inventory equipment of a refinery Sptank ⊂ S = storage site that receives products from Soutplant Srtank ⊂ S = storage site that feeds products to Sinplant SI {Æs,iæ} = product i that a site s possesses V {v} = set of vehicle modes for transportation Link {Æs,s0 ,væ} = two sites can be linked by transportation mode v Tran {Æs,s0 ,i,væ} = specific transportation tuple, denotes the feasible scenario of transportation path, product, and vehicle M {m} = set of operation modes T {t} 1, ..., Nt = range of time period of short-term models K {k} 1, ..., Nk = range of scheduling time period; its range equals one short time period. It is assumed that an operation mode at least continuously occupies a scheduling time period H {h} 1, ..., NtNk = range of aggregated scheduling time periods in the whole short-term planning horizon

Conditional Value-at-Risk Variables and Parameters

dem αc,i,t , αinv s,i,t = value-at-risk for customer dissatisfaction/inventory violation in period t βdem, βinv = user-specified confidence level for loss function less dem , αinv than αc,i,t s,i,t in period t θdem, θinv = user-specified risk aversion for customer dissatisfaction/ inventory violation inv f dem c,i,t,sc, f s,i,t,sc = loss function for customer dissatisfaction/inventory violation inv πdem c,i,t,sc, πs,i,t,sc = nonnegative adjusted loss function for customer dissatisfaction/inventory violation

Variables

Δ Dc,i,t = amount of shortage of finished product i ∈ IFP for customer c in time period t

833



ARTICLE

’ REFERENCES

(23) Barbaro, A.; Bagajewicz, M. J. Managing financial risk in planning under uncertainty. AIChE J. 2004, 50 (5), 963–989. (24) Eppen, G. D.; Martin, R. K.; Schrage, L. A Scenario Approach to Capacity Planning. Oper. Res. 1989, 37 (4), 517–527. (25) Rockafellar, R. T. Optimization of Conditional Value at Risk. J. Risk 2000, 2, 21–41. (26) Verderame, P. M.; Floudas, C. A. Operational Planning of Large-Scale Industrial Batch Plants under Demand Due Date and Amount Uncertainty: II. Conditional Value-at-Risk Framework. Ind. Eng. Chem. Res. 2010, 49 (1), 260–275. (27) Carneiro, M. C.; Ribas, G. P.; Hamacher, S. Risk Management in the Oil Supply Chain: A CVaR Approach. Ind. Eng. Chem. Res. 2010, 49 (7), 3286–3294. (28) Verderame, P. M.; Floudas, C. A. Multisite Planning under Demand and Transportation Time Uncertainty: Robust Optimization and Conditional Value-at-Risk Frameworks. Ind. Eng. Chem. Res. 2011, 50 (9), 4959–4982. (29) Verderame, P. M.; Floudas, C. A. Integration of Operational Planning and Medium-Term Scheduling for Large-Scale Industrial Batch Plants under Demand and Processing Time Uncertainty. Ind. Eng. Chem. Res. 2010, 49 (10), 4948–4965. (30) Shapiro, A.; Homem-de-Mello, T. A simulation-based approach to two-stage stochastic programming with recourse. Math. Program. 1998, 81 (3), 301–325. (31) Mele, F. D.; Guillen, G.; Espu~ na, A.; Puigjaner, L. A SimulationBased Optimization Framework for Parameter Optimization of SupplyChain Networks. Ind. Eng. Chem. Res. 2006, 45 (9), 3133–3148. (32) Law, A. M.; Kelton, W. D. Simulation Modeling & Analysis, 3rd ed.; McGraw-Hill: New York, 2000. (33) Karuppiah, R.; Martín, M.; Grossmann, I. E. A simple heuristic for reducing the number of scenarios in two-stage stochastic programming. Comput. Chem. Eng. 2010, 34 (8), 1246–1255. (34) McDonald, C. M.; Karimi, I. A. Planning and Scheduling of Parallel Semicontinuous Processes. 1. Production Planning. Ind. Eng. Chem. Res. 1997, 36 (7), 2691–2700. (35) Rockafellar, R. T.; Uryasev, S. Conditional value-at-risk for general loss distributions. J. Banking Finance 2002, 26 (7), 1443–1471. (36) Wang, W.; Ahmed, S. Sample average approximation of expected value constrained stochastic programs. Oper. Res. Lett. 2008, 36 (5), 515–519. (37) Grinstead, C. M.; Snell, J. L. Introduction to Probability, 2nd ed.; American Mathematical Society: Providence, RI, 1997.

(1) Maravelias, C. T.; Sung, C. Integration of production planning and scheduling: Overview, challenges and opportunities. Comput. Chem. Eng. 2009, 33 (12), 1919–1930. (2) Verderame, P. M.; Elia, J. A.; Li, J.; Floudas, C. A. Planning and Scheduling under Uncertainty: A Review Across Multiple Sectors. Ind. Eng. Chem. Res. 2010, 49 (9), 3993–4017. (3) Shah, N. K.; Li, Z.; Ierapetritou, M. G. Petroleum Refining Operations: Key Issues, Advances, and Opportunities. Ind. Eng. Chem. Res. 2011, 50 (3), 1161–1170. (4) Khor, C. S.; Elkamel, A.; Ponnambalam, K.; Douglas, P. L. Twostage 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 (910), 1744–1764. (5) Li, W.; Hui, C.-W.; Li, P.; Li, A.-X. Refinery Planning under Uncertainty. Ind. Eng. Chem. Res. 2004, 43 (21), 6742–6755. (6) Gu, H.; Rong, G. A Two-stage Discriminating Framework for Making Supply Chain Decisions under Uncertainties. Chem. Biochem. Eng. Q. 2010, 24 (1), P51–56. (7) Lababidi, H. M. S.; Ahmed, M. A.; Alatiqi, I. M.; Al-Enzi, A. F. Optimizing the Supply Chain of a Petrochemical Company under Uncertain Operating and Economic Conditions. Ind. Eng. Chem. Res. 2004, 43 (1), 63–73. (8) Khor, C. S.; Elkamel, A.; Douglas, P. L. Stochastic refinery planning with risk management. Pet. Sci. Technol. 2008, 26 (14), 1726–1740. (9) Yang, J.; Gu, H.; Rong, G. Supply Chain Optimization for Refinery with Considerations of Operation Mode Changeover and Yield Fluctuations. Ind. Eng. Chem. Res. 2010, 49 (1), 276–287. (10) Goel, V.; Grossmann, I. E. A stochastic programming approach to planning of offshore gas field developments under uncertainty in reserves. Comput. Chem. Eng. 2004, 28 (8), 1409–1429. (11) Jonsbraten, T. W. Optimization Models for Petroleum Field Exploitation. Thesis, Stavanger College, Stavanger, Norway, 1998. (12) Jonsbraten, T.; Wets, R.; Woodruff, D. A class of stochastic programs with decision dependent random elements. Ann. Oper. Res. 1998, 82, 83–106. (13) Goel, V.; Grossmann, I. E.; El-Bakry, A. S.; Mulkay, E. L. A novel branch and bound algorithm for optimal development of gas fields under uncertainty in reserves. Comput. Chem. Eng. 2006, 30 (67), 1076–1092. (14) Solak, S. Efficient Solution Procedures for Multistage Stochastic Formulations of Two Problem Classes. Thesis, Georgia Institute of Technology, Atlanta, 2007. (15) Balasubramanian, J.; Grossmann, I. E. Approximation to Multistage Stochastic Optimization in Multiperiod Batch Plant Scheduling under Demand Uncertainty. Ind. Eng. Chem. Res. 2004, 43 (14), 3695–3713. (16) Tarhan, B.; Grossmann, I. E. A multistage stochastic programming approach with strategies for uncertainty reduction in the synthesis of process networks with uncertain yields. Comput. Chem. Eng. 2008, 32 (45), 766–788. (17) Gupta, V.; Grossmann, I. E. Solution strategies for multistage stochastic programming with endogenous uncertainties. Comput. Chem. Eng. 2011, 35 (11), 2235–2247. (18) Gupta, A.; Maranas, C. D. Managing demand uncertainty in supply chain planning. Comput. Chem. Eng. 2003, 27 (89), 1219–1227. (19) Ben-Tal, A.; El Ghaoui, L.; Nemirovski, A. Robust Optimization; Princeton University Press: Princeton, NJ, 2009. (20) 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 (4), 931–946. (21) Mulvey, J. M.; Vanderbei, R. J.; Zenios, S. A. Robust Optimization of Large-Scale Systems. Oper. Res. 1995, 43 (2), 264–281. (22) Ahmed, S.; Sahinidis, N. V. Robust Process Planning under Uncertainty. Ind. Eng. Chem. Res. 1998, 37 (5), 1883–1892. 834


Planning under Demand and Yield Uncertainties in an Oil Supply Chain

Recommend Documents