Supply Chain Optimization for Refinery with Considerations of

276

Ind. Eng. Chem. Res. 2010, 49, 276–287

Supply Chain Optimization for Refinery with Considerations of Operation Mode Changeover and Yield Fluctuations Jiali Yang, Haijie Gu,* and Gang Rong* State Key Laboratory of Industrial Control Technology, Zhejiang UniVersity, Hangzhou, 310027, People’s Republic of China

Stochastic programming is employed to achieve optimization for the multiperiod supply chain problem in a refinery with multiple operation modes under the uncertainty of product yields. With dramatic fluctuations of product yields at the beginning of operation mode changeover, the product yields tends to stabilize after the changeover is finished. Markov chain is utilized here to describe the dynamic characteristic of product yield fluctuations. The distribution of yield fluctuation in each period is usually unknown since it depends on the decision variable of operation mode changeover. Therefore, the resulting chance constrained programming is more complicated than general situations where the distribution characteristic of stochastic variable is known in each period. This problem can be solved by the big-M method and by transforming chance constrained inequalities into a group of equivalent deterministic inequalities. This method provides a universal approach for similar chance constrained programming in which the distribution of stochastic variable depends on binary decision variables. Case studies show that the proposed modeling and solving approach can provide an effective decision-making guidance that balances confidence level and economic interests for supply chain optimization problems with multiple operation modes under yield uncertainty. 1. Introduction Supply chain risks include operational risks and disruption risks.1 Operational risks refer to inherent uncertainties such as customer demands, price of raw material or products, delays in delivery, fluctuation of product yields, process failure or defective products, etc.2,3 Among these, demand uncertainty has been well developed and widely discussed.2-8 Studies on other operation risks include balancing demand uncertainty and price fluctuation,9-12 crude oil operations under ship arrival uncertainty13 or freight rate uncertainty,14 and so on. Research about the uncertainty of product yields is relatively rare. Lababidi and Ahmed applied a scenario-based approach to discuss uncertainties in key parameters. They used three scenarios, namely “above average”, “average”, and “below average”, to present uncertainties such as product yields, with each scenario having an associated probability of occurrence. In this way, the expected cost under different scenarios and differing decision-making can be evaluated.15 This paper focuses on the influence of product yield fluctuation on decision-making in the supply chain operational level. Refinery product yields mainly depend on crude oil properties, cut points of crude oil distillation, and product blending.16-18 Product yields fluctuate dramatically during the period when there is a change of operation mode. This is largely due to the adjustment of crude oil and operation condition. The product yields gradually tend to be stable as production proceeds. Based on this fact, we assume that fluctuation of product yields is independent of the previous production condition during the period when there is a change of operation mode. Accordingly, the fluctuation of product yields within the same operation mode can be assumed to be Markov chain due to production continuity. Markov chain is a discrete process in terms of time and state and is characterized by the Markov property,19 which denotes that the future state of the process only depends on the current state. Markov chain has been widely applied to describe * To whom correspondence should be addressed. E-mail: grong@ iipc.zju.edu.cn (G.R.); [email protected] (H.G.).

production processes, such as the inventory dependence relationship in product recycle and reproduction process,20,21 production-distribution system with three nodes,22 state transition relationship of manufacturing systems,23 assessment of equipment failure and availability,24 etc. In this paper, Markov chain is applied to describe the dynamic evolution process of product yield fluctuation. However, different from the static treatment of product yield fluctuation in the literature,15 where the characteristic of product yield fluctuation is assumed to be fixed over a time period and it has nothing to do with operation condition, we consider the fluctuation of product yields from the dynamic point of view; i.e., the fluctuation characteristic is closely related to the changeover of operation modes. Approaches to dealing with optimization problems under uncertainty include scenario-based stochastic programming,25 chance constrained stochastic programming,26 simulation-based optimization,27 and fuzzy programming.28 Chance constrained programming is a very competitive tool for solving optimization problems under uncertainty. It provides a systematic analysis instrument for the trade-off between profits and reliability in decision-making.26 This method was first introduced decades ago.29,30 Its main feature is that the resulting decision ensures the probability of complying with constraints, i.e., the confidence level of being feasible. In recent years, applications of chance constrained programming have been seen in the field of process system engineering,31-35 as well as in the production planning optimization problem under supply or demand uncertainty.26,36,37 These papers mainly discuss supply or demand uncertainties without considering production fluctuation. They also assume the stochastic distribution in each period is given, including the distribution pattern and parameters and the covariance matrix that represents interactions of stochastic distribution among all periods. Therefore, these models can be easily converted into an equivalent deterministic model and solved by standard method. However, when we adopt chance constrained programming to deal with product yield uncertainty in supply chain operation, the distribution of stochastic variables in each period depends on the decision variable of operation mode changeover

10.1021/ie900968x CCC: $40.75  2010 American Chemical Society Published on Web 11/17/2009

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

Figure 1. Scheme of a refinery.

since we consider the dynamic characteristic of product yield fluctuation under multiple operation modes. Thus, converting our chance constrained programming model into an equivalent deterministic model is more complicated than normal conditions. A method is proposed in this paper to solve such a complicated problem, where a chance constrained inequality will be transformed into a group of deterministic inequalities with binary decision variables. Overall, the dynamic characteristic of product yields is considered in this paper. Markov chain is applied to represent the relationship between the dynamic characteristic of yield fluctuation and operation mode changeover. Furthermore, chance constrained programming is employed to solve supply chain operation problems under yield uncertainty. The proposed approach provides a reference for similar issues where the probability distribution of uncertainty over time depends on binary decision variables. The rest of the contents are organized as follows: the research object and concerned issue are presented in section 2 and the modeling methodology and solving approach are illustrated in section 3, followed by case studies where the proposed model is tested and experimental findings are discussed (section 4). A brief conclusion is made in the final section. 2. Problem Statements We consider a supply chain optimization problem for the process as shown in Figure 1, which is the description of a simplified refinery. The refinery process contains five main processing units: crude distillation unit (CDU), reformer, fluid catalytic cracker (FCC), gasoline blending pool (GB), and diesel blending pool (DB). Mixed crude oils are processed by these processing units and produce four final products, namely, no. 93 gasoline (G93), no. 90 gasoline (G90), diesel (D), and fuel oil (FO). Many intermediate products are generated and used in this process, including gross overhead (GO), heavy naphtha (HN), atmospheric gas oil (AGO), vacuum gas oil (VGO), bottom residua (BR), reformer gasoline (RG), crack gasoline (CG), and crack gas oil (CGO). Methyl tert-butyl ether (MTBE) is an important raw material which is used for blending. Crude oils are purchased from different oil fields, each with a unique true boiling point (TBP) curve and quality attributes (such as sulfur content, acid value, etc.). The properties of a certain crude oil usually remain stable. Refineries usually have high requirements for production stability and safety, including the quality attributes of crude oils and products, chemical and physical attributes of reactants in devices, as well as operating conditions of equipment. Several operation modes have been formed after years of process optimization and practice.38 Each operation mode tends to produce more volume of one or two products. Moreover, the mixing ratio of different crude oils, process

277

operating conditions, and product yields are relatively fixed under each operation mode. Thus, we refer to the expected product yields under each operation mode as standard yields. However, operation modes should adjust to the variation of market demands. At the time when there is a change of operation modes, process operating conditions cannot adjust immediately due to production continuity. Therefore, product yields may deviate from standard yields in the beginning, but product yield gradually approaches standard yield as production continues. Since the properties of crude oil may have certain levels of fluctuation, operators also have some autonomy to make temporary adjustments to control systems to meet the process operating conditions. Therefore, even if the process is stabilized, product yields may still fluctuate in the vicinity of standard yields. We assume that the fluctuation of product yields in a current period is only relevant to the previous scheduling period during the process that product yields gradually stabilize. Based on this description, the production uncertainty can be modeled as a Markov process. More details are addressed in the next section. Under given market demands, production uncertainty will cause uncertainty in product stock. In bad situations, the safety stock level may be exceeded and the refinery may fail to delivery products to customers on time. Therefore, a robust optimization approach is proposed when making short-term supply chain operation decisions; it allows operation mode changeover and considers production uncertainty. Chance constrained programming is introduced to achieve robust decision-making. Its objective is to ensure that the probability of inventory within the safety stock level reaches a predetermined confidence level. Usually, chance constrained programming considers a continuous-time Markov chain. However, this paper aims at a discretetime Markov chain. Besides, the distribution of stochastic variable is dependent on the decision variable due to the changeover of operation modes. Therefore, probability constraints in this paper cannot be simply converted into equivalent inequalities as a general solution algorithm. Note that, for simplicity, supply chain refers to short-term supply chain in the latter sections unless specified otherwise. 3. Chance Constrained Programming under Product Yield Uncertainty Two time scales are defined since the supply chain period is larger than the production scheduling period. The supply chain period is denoted by t, with total consideration of Nt periods. The production scheduling period is denoted by k, with Nk scheduling periods in each supply chain period. For simplicity, the compound scheduling period 〈t,k〉 is defined to represent the number k scheduling period of the number t supply chain period. To facilitate computation, a cumulative scheduling period j is introduced to indicate the continuous count of scheduling periods in the whole supply chain period. The formulation between cumulative scheduling period j and compound scheduling period 〈t,k〉 is j ) (t - 1)Nk + k, where the range of j is 1, ..., Nt Nk. 3.1. Markov Chain of Product Yield Fluctuation. Yi,m is the standard yield of product i under operation mode m, while Yî,m,t,k is the actual yield of product i under operation mode m in compound scheduling period 〈t,k〉. ξi,t,k represents the actual yield fluctuation of product i compared with the standard yield in compound scheduling period 〈t,k〉. Its range is a series of discrete values Λi ) {λi,1, λi,2, ..., λi,Nsi}, which is the Markov state space. The total state number of this space is therefore defined as Nsi. Thus

278


Yî,m,t,k ) Yi,m + ξi,t,k

(1)

We assume that the production loss in a refinery remains a constant. Then, the yield fluctuation of all Np types of products should meet the following constraint:

∑ξ

i,t,k

) 0,

∀t, k

(2)

i

It is assumed that the yield fluctuations of (Np - 1) types of products are independent of each other, while the remaining product yield entirely depends on them. Without the loss of generality, we suppose that the yield fluctuation of the first (Np - 1) types of products are independent of each other, while the yield fluctuation of product Np can be derived from eq 2. 3.1.1. Internal Transition within the Same Operation Mode. According to the assumption in section 2, if there is no changeover of operation modes and the current product yield fluctuation is known, then the product yield fluctuation of the next scheduling period is independent of historical states. Thus, the Markov chain can be applied to describe the state transition of yield fluctuation ξi,t,k within the same operation mode, where i is among the first (Np - 1) types of products. Based on the property of Markov,19 we can derive Pr{ξi,t,k+1 ) λi,b |ξi,t,k ) λi,a, ...,ξi,1,1 ) λi,x} ) Pr{ξi,t,k+1 ) λi,b |ξi,t,k ) λi,a}

(3)

where λi,a, λi,b, λi,x ∈ Λi. In eq 3, ξi,t,k+1 equals ξi,t+1,1 when k ) Nk. The same treatment will be done without repeating the details if a similar situation of period overflow happens in the following sections. The Markov state transition is assumed to be time-homogeneous, that is Pr{ξi,t,k+1 ) λi,b |ξi,t,k ) λi,a} ) Pr{ξi,t,k ) λi,b |ξi,t,k-1 ) λi,a},

∀t, k

(4)

For simplicity, the above one-step transition probability can be presented as pi,a,b ) Pr(ξi,t,k+1 ) λi,b |ξi,t,k ) λi,a)

(5)

The state transition matrix of product i within the same operation mode can be indicated by

(

λi,1 · · · λi,Nsi λi,1 pi,1,1 ... pi,1,Nsi Ρi ) l · ·· l l λi,Nsi pi,Nsi,1 · · · pi,Nsi,Nsi

)

(6)

where the abth entry of Ρi stands for the probability of product yield fluctuation transitioning from state λi,a to state λi,b in a single step. Based on theorem 11.1 specified by Grinstead and Snell,19 (n) the n-step transition probability pi,a,b equals the abth entry of n n matrix Ρi , where Ρi is the nth power matrix of state transition matrix Ρi. In addition, Ρ0i is a unit matrix with Nsi × Nsi dimensions. 3.1.2. State Transition Considering Operation Mode Changeover. It is assumed that only one operation mode can be executed in each scheduling period, and changeover of operation modes can only take place at the junction of two scheduling periods. For simplicity, we suppose that product yield fluctuations under different operation modes have the same state space and state transition matrix. In addition, the initial distributions are supposed to be the same when changing into

a new operation mode. Pbi(1) is the initial probability distribution vector of product i’s yield fluctuation when the operation mode has been changed. The corresponding probability locating at state λa is denoted as pbi,a(1). Therefore Pbi(1) ) [pbi,1(1), pbi,2(1), ...,pbi,Nsi(1)]

(7)

is a vector with 1 × Nsi dimensions. After a new operation mode has been executed continuously to the τth scheduling period, these τ states of product yield fluctuation will form a path. There may be altogether (Nsi)τ possible paths. The state transition path of λi,x1 f λi,x2 ... f λi,xτ (λi,x1, λi,x2 ... λi,xτ ∈ Λi) is defined as e. The probability of path e is therefore τ-1 pi,xl,xl+1. pbi,x1(1)∏l)1 To facilitate discussion, the operation mode in the first scheduling period (t ) 1, k ) 1) of the planning horizon is always treated as a new one. In this way, the probability distribution of product i’s yield fluctuation is fixed at Pbi(1) in the first scheduling period. This assumption is made only for the convenience of algorithm illustration. In fact, subsequent computations can be conducted if the initial operation mode of the entire horizon and the number of execution periods are given. So far, as long as the list of operation mode changeover is given throughout optimization horizon, we can readily calculate the probability distribution of each product’s yield fluctuation in each scheduling period. In addition, it is worth noting that the probability distribution of Np products’ yield fluctuation can be solved from eq 2 according to the distribution of the other (Np - 1) products’ yield fluctuations. 3.1.3. Parameter Evaluation and Hypothesis Testing. In the previous sections, Markov chain is adopted to describe the dynamics of product yield fluctuations with changeover of operation modes. In this subsection, a method is proposed to evaluate the parameters about Ρi and Pbi(1) and test the hypotheses about Markov chain based on the historical data of the refinery. First, the mean and standard deviation values of each product yield under every operation mode should be calculated to determine the standard yields and the discrete states of yield fluctuation. Then, the historical state value of yield fluctuation in each scheduling period can be obtained. The data should be split into two parts. One is used for calculating parameters, and the other is used for testing hypotheses. The elements of Ρi and Pbi(1) can be evaluated by the quotient of occurrence rate of the corresponding different items. For example, element pi,a,b of Ρi can be computed by fi,a,b/fi,a, where fi,a represents the frequency of state λi,a except the ones just before the changeover happens and fi,a,b represents the frequency of state λi,b in the condition that the state locates at λi,a in the previous scheduling period, and there is no changeover during these two periods. After the parameters of Ρi and Pbi(1) are obtained, the hypothesis that the product yield fluctuation based on Markov chain with the evaluation parameters should be tested from the testing data. To test the goodness of fit, the chi-square distribution (see chapter 9 of ref 39) will be adopted. In such a chi-square goodness-of-fit test, the expected frequency inferred from Ρi or Pbi(1) will be compared with the observed frequency from the testing data. If the two numbers are close enough, the value of approximated chi-square obtained from their differences tends to be smaller than the standard chi-square value with a specified confidence level. In this way, the hypothesis is proven to be true. 3.2. Deterministic Supply Chain Optimization. Before discussing stochastic programming, the supply chain optimization model without considering product yield uncertainty is


presented in this section. Since this paper focuses on production uncertainty, external networks such as transportation, procurement, and distribution are not taken into account. The adding or omitting of these factors will not impact the core issues we are dealing with here. Our deterministic optimization model is a common approach and is based on Li and Karimi’s production planning model.36,40,41 The description of the supply chain operation problem is specified as follows: (1) Known conditions or constraints: the number of supply chain periods throughout the optimization horizon and the number of scheduling periods in each supply chain period; the price of raw materials and products; changeover costs of operation modes; allowable purchasing amount of raw materials within the whole optimization horizon; mixing ratio of crude oils and product yields under each operation mode; refinery capacity constraints; initial inventory of each product as well as safety stock levels; demand for all products in each supply chain period (2) Objective function: pursuit of the maximization of profit under the premise that customer demands are satisfied as much as possible, since customer satisfaction level has a great impact on market share in the long-term perspective (3) Decision variable: purchasing amount of raw materials; refinery’s processing volume per period; list of operation mode changeover; sales volume of products. It is assumed that the processing volume Q is stable in each scheduling period throughout the optimization horizon since a refinery has high demands for maintaining safety and smooth production. Therefore, frequent changeover of operation modes is also forbidden. Each operation mode should last for at least N1 scheduling periods. We also assume that customer demands should be met during the required period without delay. In other words, back orders are not allowed. The mathematical formulation of the deterministic supply chain operation optimization model is given as follows: objective function: max

∑ γ (SL i

i,t

i,t)

-

∑γR

i i,t

- δQ(Nt)(Nk) -

t

∑ ηH

(t-1)Nk+k

t,k

(8) subject to: Ri,t )

∑Q

m,t,kXi,m

(9)

m,k

∑Q

(10*)

QL e Q e QU

(11)

0 e Qm,t,k e QU

(12)

Q e Qm,t,k + (1 - Zm,t,k)M

(13)

Qm,t,k e Q + (1 - Zm,t,k)M

(14)

Qm,t,k e Zm,t,kM

(15)

Gi,t )

m,t,kYi,m

m,k

∑Z

m,t,k

)1

(16)

m

CHm,t,k g Zm,t,k - Zm,t,k-1,

∀k g 2

(17)

CHm,t,k g Zm,t,k - Zm,t-1,K,

∀k ) 1

(18)

H(t-1)Nk+k )

∑ CH

m,t,k

m

(19)

279

j+N1-1

∑

Hj' e 1,

∀j ) 1, 2, ...,(Nt)Nk - N1 + 1

(20)

Hj ∈ {0, 1}

(21)

j')j

CHm,t,k,

Zm,t,k,

0 e SLi,t e Di,t

(22)

Vi,t ) Vi,t-1 + Gi,t - SLi,t

(23*)

VLi e Vi,t e VUi

(24*)

In the above formulas, the various symbols are specified as follows: SLi,t, sales volume of product i in supply chain period t; Ri,t, purchasing amount of raw material i in supply chain period t; γi, price of material i; Q, processing amount of total raw materials in each scheduling period; QU and QL, upper and lower bounds of processing amount respectively; δ, operation cost for per processing amount of raw materials; η, operation mode changeover cost; Xi,m, mixing ratio of raw material i under operation mode m; Qm,t,k, processing amount of total raw materials under operation mode m in compound scheduling period 〈t,k〉 (its value should be Q or 0); Gi,t, production volume of product i in supply chain period t; Zm,t,k, binary variable that denotes whether operation mode m is executed in compound scheduling period 〈t,k〉; CHm,t,k, binary variable that denotes whether operation mode is changed to m in compound scheduling period 〈t,k〉; M, a big constant; Di,t, demand for product i in supply chain period t; Vi,t, inventory level of product i at the end of supply chain period t; ViU and ViL, inventory upper and lower bounds of product i. The objective function in eq 8 is to maximize profit, which equals the product sales revenue subtracting the purchasing cost of raw materials, processing cost, and operation mode changeover cost. Equation 9 indicates that the purchasing amount of raw material i equals the total processing amount times the corresponding mixing ratio of this material. Equation 10* denotes that the production volume of product i equals the total processing amount times the standard yield of this product. Equations 11-15 represent that the processing amount remains stable at Q in each scheduling period, and is subject to the processing capacity constraint. If operation mode m is currently executed, then Qm,t,k ) Q; otherwise, Qm,t,k ) 0. Equations 13-15 realize this logic by implementing the big-M method, where M is a big constant and should at least meet the constraint of M > QU. The above logic works in this way: When Zm,t,k ) 1, eqs 13 and 14 jointly realize Qm,t,k ) Q, combined with M > QU; eq 15 becomes a redundant expression for the existence of eq 12. When Zm,t,k ) 0, eq 13 becomes a redundant expression for the upper bounds of eq 11 and M > QU. Similarly, eq 14 will lose its impact since its constraint is more slack than the upper bound of eq 12 while eq 15 ensures Qm,t,k ) 0. Equation 16 denotes that only one operation mode can be adopted in each scheduling period. Hj is a binary variable; it is used to describe whether there is a change of operation mode in compound scheduling period 〈t,k〉 or in its corresponding cumulative scheduling period j ) (t - 1)Nk + k. Equations 17-19 realize the calculation of Hj, where eqs 17 and 18 calculate whether operation modes have been changed to m in compound scheduling period 〈t,k〉. Moreover, eq 19 calculates whether there is any changeover of operation mode during this period. Equation 20 indicates the frequency constraint of operation mode changeover. As it shows, only one changeover is allowed in any continuous N1 scheduling periods to guarantee steady production. Equation 22 represents products must be supplied to the customer within the

280


specified period in the contract. Equation 23* shows that the inventory held at the end of period t is equal to that held at the end of the previous period (t - 1) plus production volume during period t and minus the sales volume during period t. Equation 24* denotes that the inventory level should meet the upper and lower bounds. 3.3. Chance Constrained Programming Model with Product Yield Uncertainty. In the previous subsection, the deterministic mathematical model is presented. Here the uncertainty of product yields is taken into account and the deterministic model should be revised properly. In order to cope with temporary changes in product demands (additional order or withdraw order), an extreme inventory level such as complete fullness or emptiness is not allowed when making a supply chain plan. Instead, upper and lower safety bounds of product inventory are set to increase flexibility. In the deterministic model, this requirement is expressed by constraint eq 24*. However, safety bounds of product inventory can no longer be fully guaranteed when considering the uncertainty of product yields, unless the worst scenario method is employed.42 However, the worst scenario method can not only severely impair profits, but also is hard to define when in the situation of multiple products, multiple operation modes, and dual bounds of safety stock levels. Under the case that safety stock levels cannot be fully ensured, it is hoped that the probability of safety stock levels being maintained reaches a certain extent. Therefore, chance constrained programming is employed to deal with the uncertainty of product yields. We assume that the probability of product i’s safety stock levels being maintained in each supply chain period should not be lower than confidence level Ri (0 e Ri e 1). This constraint can be expressed as follows: Pr{VLi e Vî,t e VUi } g Ri

(25)

where Vî,t denotes product i’s inventory with consideration of yield uncertainty. There is a trade-off between decision reliability and profitability. If the decision is robust at resisting risks, in other words, more conservative, then Ri should be relatively large while the profit margin would be small. Conversely, if the decision maker encourages risk taking, or is more aggressive, then Ri should be relatively small and the profit margin is also larger. However, once the realization of stochastic variable is worse than expected, the negative impact on customer satisfaction could be serious and enduring. Hence, selection of Ri depends on the decision maker’s risk preference, as well as the trade-off curve between decision reliability and profit. In the deterministic model, production volume is calculated by the standard yield Yi,m in eq 10*. When yield uncertainty is considered, the ˆ i,t should be obtained from the uncertain production volume G product yield Yî,m,t,k. This uncertainty is transformed further to Vî,t through eq 23*. In brief, the uncertainty of the model is caused by the product yield fluctuation ξi,t.k, which is also the fluctuation ofYî,m,t,k. Since the constraint of probability is bound on the safety inventory level, only the equations with asterisks (*) in subsection 3.2 need to be revised when transforming the deterministic model into the stochastic model. The relation between Vî,t and ξi,t,k can be formulated as follows based on equations in subsection 3.2 and eq 1:

ˆ i,t - SLi,t Vî,t ) Vî,t-1 + G Qm,t,k(Yi,m + ξi,t,k) - SLi,t ) Vî,t-1 +

∑ ∑ ∑Q ∑ ∑Q m,k

) Vi,0 +

+ ξi,t',k) -

m,t',k(Yi,m

t'et m,k

) Vi,0 +

∑ SL

i,t'

t'et

m,t',kYi,m

∑ ∑ Qξ

+

i,t',k

t'et m,k

t'et

-

k

∑ SL

i,t'

t'et

(26) where Vi,0 is the initial inventory of the whole optimization horizon. The relation between Qm,t,k and Q is employed in the derivation process. The corresponding stochastic chance constrained model can be derived by replacing eqs 10*, 23*, and 24* in the deterministic model with the following equation (27). Pr{VLi e Vi,0 +

∑ ∑Q

∑ ∑ Qξ ∑ SL e V } g R

+

m,t',kYi,m

t'et m,k

i,t',k

t'et

k

U i

i,t'

i

(27)

t'et

3.4. Solving the Chance Constrained Programming Model. There are two inequalities of upper and lower bounds in eq 27. Since the solving processes of these two inequalities are very similar, we first illustrate the solving process of the lower bound, which is the following inequality: Pr{VLi e Vi,0 +

∑ ∑Q

∑ ∑ Qξ ∑ SL

+

m,t',kYi,m

i,t',k

t'et m,k

t'et

-

k

i,t'}

g Ri (28)

t'et

The following equivalent equation can be obtained by transformation: Pr{

∑ ∑ξ

i,t',k

t'et

∑ SL ∑ ∑Q

g (VLi - Vi,0 +

k

i,t'

-

t'et

m,t',kYi,m)/Q}

g Ri

(29)

t'et m,k

Equation 29 does not have the standard form of chance constrained inequality as discussed in refs 26, 36, and 37 because the distribution of ξi,t′,k is unknown and it depends on the changeover of operation mode. As a result, eq 29 cannot be simply transformed into the corresponding deterministic inequality. According to the Markov chain model of product yield fluctuation in subsection 3.1, the following conclusion holds true: the probability distribution of the product yield fluctuation throughout the whole optimization horizon can be acquired as long as the decision variable of operation mode changeover list is given. This conclusion can also be applied to obtain the solution algorithm of eq 29. The main approach is to enumerate all possible combinations of the operation mode changeover, list the corresponding equivalent deterministic inequality of each changeover sequence, and then synthesize. If one possible sequence of the operation mode changeover is known, as well as the initial probability distribution Pbi(1) with operation mode changeover and the state transition matrix Pi, then the probability of ξi,t.k at any state λi,a (λi,a ∈ Λi) in any compound scheduling period 〈t,k〉 can be obtained. Two abbreviations are made as follows for the convenience of discussion: wi,t ) VLi - Vi,0 +

∑ SL

i,t'

t'et

-

∑ ∑Q t'et m,k

m,t',kYi,m

(30)


zi,t ) wi,t /Q

(31)

Meanwhile, F(zi,t) is employed to represent the probability distribution function: F(zi,t) ) Pr{

∑ ∑ξ

i,t',k

t'et

e zi,t}

(32)

k

Then, solving eq 29 is equivalent to solving the following equation (33): F(zi,t) e 1 - Ri

(33)

The next step is to calculate F(zi,t). Since the distribution of stochastic variable ξi,t,k is discrete, F(zi,t) can be obtained by the accumulation method. From the first compound scheduling period 〈1,1〉 to period 〈t,K〉, a state tree of stochastic variable ξi,t,k has formed with (Nsi)j branches within j ) (t - 1)Nk + K cumulative scheduling periods. Any branch Tri,t,e(e ∈ {1, 2, ..., (Nsi)j}) of the state tree can be denoted by sequence λi,x1 f λi,x2 ... f λi,xj (λi,x1, λi,x2 ... λi,xj ∈ Λi). The probability of this branch can be expressed by PTri,t,e j-1

∏ pb

(34)

pbi,xl(1), l ∈ {l|Hl ) 1} pi,xl-1,xl, l ∈ {l|Hl ) 0}

(35)

PTri,t,e ) pbi,x1(1)

i,xl,xl+1

l)1

where pbi,xl-1,xl )

{

Take Figure 2; for example, the sequence of states with circle symbol represents one branch of the state tree, while the one with triangle symbol represents another branch. H1 and H4 of eq 35 both equal 1 in Figure 2 since there is a change of operation mode at scheduling periods 1 and 4. The cumulative yield fluctuations ∑t′et ∑k ξi,t′,k of branch Tri,t,e can be denoted as the cumulative state value STri,t,e (overall deviations from the standard yields); then j

STri,t,e )

∑λ l)1

i,xl

(36)

It should be noted that the cumulative state values of many branches are the same. For example, the cumulative state values

Figure 2. Relation between state transition and operation mode changeover.

281

of both branch e1 and e2 equal 2 in Figure 2. It is assumed that set Ωi,t is made up of all unequal cumulative state values. The maximum number of elements in set Ωi,t is computed first. The approach is to randomly select j elements from Nsi numbers of states; each state can be selected repetitively. This combination can be noted as combinations with repetition by selecting j j+Ns elements from Nsi states. The number of combinations is CNsi-1i-1 ) (j + Nsi - 1)!/[j!(Nsi - 1)!], where “!” is the factorial operator. This conclusion is addressed in chapter 6 of Nievergelt’s book.43 The actual number of elements in set Ωi,t may be even smaller since the sum of different combinations of j states may be the same. In brief, the number of elements that set Ωi,t contains is far smaller than the number of branches in the state tree. The branches with the same cumulative state values are merged by adding their probabilities together. In this way, the probability of each element in set Ωi,t is obtained. The elements in set Ωi,t is then sorted ascendingly. Combined with the probability of each element, the solution of eq 33 can be determined: zi,t e F-1(1 - Ri)

(37)

where F-1(1 - Ri) must be an element of set Ωi,t. It is worth pointing out that, for each t, eq 33 may have multiple solutions since there may be several possible sequences of operation mode changeover. Each solution is necessarily corresponding to one or more changeover sequences. Therefore, eq 28 is different from the situation with a given distribution of stochastic variables. Instead, it is transformed into a series of deterministic inequalities. Each inequality contains the information of decision variable Hj which is the changeover sequence of operation modes. Of course, the ultimate decision making can only be one of these changeover sequences. In other words, although eq 28 generates a series of inequalities, only the one that matches the selected sequence of operation mode changeover will take effect; the rest of the inequalities will be invalid. The logical relation of “K out of N constraints must be satisfied” can be resolved by the big-M method. The approach of transforming eq 28 into several deterministic inequalities by big-M method is specified in Case 1. So far, the situation of the first (Np - 1) types of products has been addressed. The yield fluctuations of these products are assumed to be independent of each other. Based on assumptions in section 3.1, the yield fluctuation of product Np is related to other products and it can be calculated from eq 2. Thus, the computation of product Np’s probability is different from that of others. First, for the first (Np - 1) types of products, the cumulative yield fluctuations ∑t′et ∑k ξi,t′,k in set Ωi,t and the probability of each element can be obtained under possible sequences of operation mode changeover from compound scheduling period 〈1,1〉 to 〈t,K〉. We assume that element ai is randomly selected from set Ωi,t. Its probability of occurrence is pai. According to eq 2, the state value of product Np should be (-∑ieNp-1 ai). Its probability is ∏ieNp-1 pai. By traversing each element in set Ωi,t and merging the probability of product Np’s branches which have the same cumulative state value, set ΩNp,t of product Np and the probability of each element can be obtained. The pseudocode for traversing the possible product mode changeover sequence, solving the state tree, and calculating the inverse function of the probability distribution F-1(1 - Ri) is given elsewhere. Detailed case studies of the above processes are shown in the next section.

282


CHm,t,k e 2 - Zm,t,k - Zm,t,k-1,

Table 1. Parameters of Product 1’s Yield Fluctuation in Case 1 parameter

value

Markov state space Λ1 initial probability distribution when operation mode is changed Pb1(1) state transition matrix when operation modes executed continuously Ρ1

[-0.02 0 0.02] [0.35 0.45 0.2]

-0.02 0 0.02

(

-0.02 0 0.02 0.4 0.3 0.3 0.2 0.55 0.25 0.1 0.4 0.5

)

∑ ∑Q

m,t',kYi,m

+

t'et m,k

∑ ∑ Qξ ∑ SL t'et

i,t',k

e

VUi }

g Ri

(38)

t'et

Similarly, we assume that z'i,t ) (VUi - Vi,0 +

∑ SL

i,t'

∑ ∑Q

-

m,t',kYi,m)/Q

t'et

t'et m,k

(39) F(z'i,t) ) Pr{

∑ ∑ξ

i,t',k

t'et

e z'i,t}

(40)

k

However, the solution of eq 40 is quite different from that of eq 33. It should be F(z'i,t) g Ri

(41)

The final format of the expression is shown as eq 42: z'i,t g F-1(Ri)

(42)

The details for solving eq 42 are not specified here since it is very similar to the solving process of eq 37. It is interesting to note that, in the deterministic model, the variable of operation mode changeover CHm,t,k can be readily calculated by the variable of operation mode execution Zm,t,k from eq 17, eq 18, and the statements about binary variables in eq 21. For instance, if both two items on the right side of eq 17 equal 1, resulting in CHm,t,k g 0, CHm,t,k will automatically take the value of 0 to reduce the cost of operation mode changeover. This fits well with our logic. However, in the stochastic model, when transforming the probability constraint equation into an equivalent deterministic model, the value of CHm,t,k does not only affect the realization of uncertainty, but also affects product yields which further influence the production volume besides the changeover cost. Therefore, depending on the trade-off between the earnings from product yields and changeover cost, CHm,t,k may take a value of either 0 or 1 in the stochastic model if both items on the right side of eq 17 equal 1, instead of definitely taking the value of 0 as in the deterministic model. In order to guarantee an accurate presentation of logical relations under the stochastic situation, a more rigid method is adopted to implement an equivalent XOR (exclusive-or) relation by replacing eqs 17-19 with the following equations (43-51): CHm,t,k e Zm,t,k + Zm,t,k-1,

∀k g 2

(45)

CHm,t,k g Zm,t,k-1 - Zm,t,k,

∀k g 2

(46)

CHm,t,k e Zm,t,k + Zm,t-1,K,

∀k ) 1

(47)

∀k ) 1

(48)

CHm,t,k g Zm,t,k - Zm,t-1,K,

∀k ) 1

(49)

CHm,t,k g Zm,t-1,K - Zm,t,k,

∀k ) 1

(50)

H(t-1)Nk+k )

1 2

∑ CH

m,t,k

(51)

m

4. Case Studies

-

k

i,t'

(44)

CHm,t,k g Zm,t,k - Zm,t,k-1,

CHm,t,k e 2 - Zm,t,k - Zm,t-1,K,

In addition, the upper bound in eq 27 is quite similar to the lower bound we mentioned above. The difference lies in the direction of inequality and the solving of the inverse function from the probability distribution function. The upper bound is given by the following inequality: Pr{Vi,0 +

∀k g 2

∀k g 2

(43)

4.1. Case 1. A simple example is designed to illustrate the computation of product yield probability based on the discrete Markov chain under multiple operation modes, as well as transformation of stochastic constraints which relate to the decision variable into equivalent equations which can be solved by mixed integer linear programming (MILP). The scale of the model is as follows: Nt ) 2, Nk ) 2; one raw material and two products; two operation modes; minimum changeover period N1 ) 3; each product’s yield fluctuation has three possible states. We assume that the distribution of product 1’s yield fluctuation is independent, while the distribution of product 2’s yield is related to product 1’s according to eq 2. The yield fluctuation parameters of product 1 are shown as Table 1. It should be noted that this case can be solved manually instead of by use of software due to the small scale and limited number of variables. Based on the constraint of minimum changeover period, it can be seen that within four compound periods from 〈1,1〉 to 〈2,2〉 the sequence of operation mode should either be [1 0 1 0] or [1 0 0 1], which can be denoted by “seq1” and “seq2”, respectively. “1” signifies there is a change of operation modes. For simplicity, this value is fixed at “1” in the first compound period according to the assumption in subsection 3.1.2. In the case of seq1, the state tree of the stochastic variable’s distribution is shown in Table 2. State transition path 2s3s1s2 represents product 1’s yield fluctuation locating at state 2(λ1,x1 ) 0), 3(λ1,x2 ) 0.02), 1(λ1,x3 ) -0.02), and 2(λ1,x4 ) 0), respectively, in the four compound periods. The probability of each state transition path is calculated according to eqs 34 and 35. The cumulative state value of each path is calculated according to eq 36. For instance, for state transition path 2s3s1s2, the initial probability of state 2 is pb1,x1(1) ) 0.45; from state 2 to state 3, the one-step transition probability within the same operation mode should be pb1,x1,x2 ) p1,x1,x2 ) 0.25; from state 3 to state 1, the one-step transition probability with operation mode changeover should be pb1,x2,x3 ) pb1,x3(1) ) 0.35; from state 1 to state 2, the one-step transition probability within the same operation mode should be pb1,x3,x4 ) p1,x3,x4 ) 0.3. Therefore, PTr1,2,e ) Table 2. State Tree under the Changeover Sequence seq1 in Case 1 state transition path

probability of this path

corresponding overall state value of this path

1s1s1s1 1s1s1s2 ... 2s3s1s2 ... 3s3s3s3

0.35 · 0.4 · 0.35 · 0.4 0.35 · 0.4 · 0.35 · 0.3 ... 0.45 · 0.25 · 0.35 · 0.3 ... 0.2 · 0.5 · 0.2 · 0.5

-0.02 - 0.02 - 0.02 - 0.02 -0.02 - 0.02 - 0.02 + 0 ... 0 + 0.02 - 0.02 + 0 ... 0.02 + 0.02 + 0.02 + 0.02


z1,1 e (seq1)M - 0.02

Table 3. Cumulative State Distribution Vector and the Corresponding Probability until period t cumulative state vector Ω1,t 1 2

283

corresponding probability vector of Ω1,t

[-0.04, -0.02, 0, 0.02, 0.04] [0.14, 0.195, 0.3725, 0.1925, 0.1] [-0.08, -0.06, -0.04, -0.02, [0.0196, 0.0546, 0.14233, 0, 0.02, 0.04, 0.06, 0.08] 0.19918, 0.24183, 0.18241, 0.11156, 0.0385, 0.01]

pb1,x1(1)∏3l)1 pbi,xl,xl+1 ) 0.45 · 0.25 · 0.35 · 0.3. STr1,2,e ) ∑4l)1 λ1,xl ) 0 + 0.02 - 0.02 + 0 ) 0. Moreover, the cumulative state distribution vector Ω1,t in each supply chain period t and the corresponding probability of each element are given in Table 3. Ω1,t contains all the elements that have unequal cumulative state values. The probability of each element in Ω1,t is obtained by adding the probabilities of the state transition paths which have the same cumulative state values. For instance, element -0.06 at period t ) 2 in Ω1,t contains four state transition paths, namely, 1s1s1s2, 1s1s2s1, 1s2s1s1, and 2s1s1s1, respectively. Therefore, the probability of this cumulative state value equals the sum of these four state transition paths’ probability, which is 0.35 · 0.4 · 0.35·0.3+0.35·0.4·0.45·0.2+0.35·0.3·0.35·0.4+0.45·0.2·0.35·0.4 ) 0.0546. According to Table 3, the reverse function F-1(1 - R1) of the probability distribution with any given confidence level R1 can be resolved. For instance, if R1 ) 0.79, z1,1 in eq 37 should be -0.02 in the cumulative state vector Ω1,1 since 0.14 < 1 0.79 < 0.14 + 0.195. z1,2 should be -0.04 in Ω1,2 since 0.0196 + 0.0546 < 1 - 0.79 < 0.0196 + 0.0546 + 0.14233. Similarly, if given R1 ) 0.55, it can also be obtained from Table 3 that z1,1 ) 0 and z1,2 ) 0. In the case of seq2, product 1’s state tree, the cumulative state distribution vector in each supply chain period and the corresponding probability can be obtained following the above approach. Similarly, if R1 ) 0.79, z1,1 and z1,2 can also be computed, which is another set of solution. To distinguish the two cases, the solutions of z1,1 and z1,2 under seq1 are expressed as z1,1(1) and z1,2(1) and the solutions under seq2 are expressed as z1,1(2) and z1,2(2). Then, z1,1(2) ) -0.02 and z1,2(2) ) -0.02. This is the coupling situation of the stochastic variable’s state distribution and the decision variable where the distribution of stochastic variable over time cannot be determined before decision-making. Therefore, the relation between the solution of z1,1 and z1,2 and the decision variable of operation mode changeover needs to be formulated. All the possible changeover sequences can be expressed by variable Hj, which denotes whether there is a change of operation mode: seq1 ) H1 + H2 + H3 + H4

(52)

seq2 ) H1 + H2 + H3 + H4

(53)

where Hj is the complementary operation of Hj and it is equal to (1 - Hj). Based on this notation of changeover sequence, only the value of the selected changeover sequence equals 0 at any time, while the rest are not equal to 0. The reason is that each item on the right side of these formulations is either 0 or 1. The sum of these items will be 0 only if each item is 0. Since all changeover sequences are unique, there must be at least one item equal to 1 on the right side of each unselected sequence. Thus, the value of the unselected sequences would always be larger than 0. For example, if sequence [1 0 1 0] is selected, eq 52 equals 0 while eq 53 equals 2. The relation between zi,t and the changeover sequence can by expressed by the big-M method as follows:

z1,1 e (seq2)M - 0.02 z1,2 e (seq1)M - 0.04 z1,2 e (seq2)M - 0.02

(54)

where M is a big constant. Based on eq 31, eq 54 can be transformed into four equivalent linear inequalities as shown in eq 55, which can be solved directly by an MILP solver. In the final solution, either seq1 or seq2 will be chosen as the decision variable of operation modes changeover. Thus, only two of the four formulas in eq 55 are valid. w1,1 e (seq1)M - 0.02Q w1,1 e (seq2)M - 0.02Q w1,2 e (seq1)M - 0.04Q w1,2 e (seq2)M - 0.02Q

(55)

So far, the inequalities of product 1’s probability constraints have been successfully transformed into equivalent deterministic linear inequalities. The coupling relation between the probability distribution function and the decision variable of operation mode changeover has been linearly separated. The next step is to transform the inequalities of product 2’s probability constraints. Instead of conducting an independent computation based on the Markov chain and sequences of operation mode changeover, the probability distribution of product 2 can be derived from product 1’s distribution function directly. According to eq 2, the distribution of the last product is a joint probability distribution if the number of all products is no less than three. This is also the situation in case 2. The joint probability is relatively easy to solve in the case of discrete probability distribution. The cumulative state vector of product 2 in each supply chain period and its corresponding probability are shown in Table 4. These results can be derived directly from Table 3. For instance, the total production loss of the two products should meet constraint eq 2, while the probability of product 1’s cumulative state 0.04 at period ) 1 is 0.1 according to Table 3. Therefore, the probability of product 2’s cumulative state -0.04 at period t ) 1 should also be 0.1. The following computing process of product 2 is very similar to product 1’s, including solving the reverse function of the probability distribution function and transforming probability constraint inequalities into equivalent linear inequalities. Details are not specified here. 4.2. Case 2. Case 1 is conducted to demonstrate the detailed calculation process of the proposed method, while in case 2 the effect of different confidence levels in decision-making under the probability constrained method is analyzed. In this case, ILOG Cplex is used as the optimization solver. The proposed optimization model was complied by ILOG OPL language with software versions ILOG Cplex 11.0 and ILOG OPL Development Studio 5.5, respectively. The resulting MILP model has 718 constraints and is solved using Cplex in 0.75 s on a computer with 2.66 GHz Intel Core 2 Duo and 2 GB RAM. The following processes are realized by matlab programming with a total computation time of 55 s, including traversing sequences of operation mode changeover, solving the probability of Markov state transition and the inverse function of probability distribution function, etc. In the following case study, only the inventory lower bound is taken into account since the approaches for solving upper bound and lower bounds in eq 27 are quite

284


Table 4. Cumulative State Distribution Vector of Product 2 and Its Corresponding Probability until Each SC Period until period t 1 2

cumulative state vector Ω1,t

corresponding probability vector of Ω1,t

[-0.04, -0.02, 0, 0.02, 0.04] [0.1, 0.1925, 0.3725, 0.195, 0.14] [-0.08, -0.06, -0.04, -0.02, [0.01, 0.0385, 0.11156, 0.18241, 0, 0.02, 0.04, 0.06, 0.08] 0.24183, 0.19918, 0.14233, 0.0546, 0.0196]

Table 5. Some Basic Configurations of Material

material i

material type and numbering

price γi ($/bbl)

initial inventory IVi,t-1 (kbbl)

R1 R2 R3 MTBE G93 D FO G90

raw material 1 raw material 2 raw material 3 raw material 4 product 1 product 2 product 3 product 4

70 62 65 92 105 100 70 97

s s s s 150 70 100 80

lower bound of safety stock IViL (kbbl) s s s s 50 50 60 50

similar. We assume that the inventory upper bound is large enough and the fluctuation of production will not pose a danger to upper safety stock levels. The short-term supply chain optimization problem of a medium-sized refinery is considered. The refinery purchases three types of crude oil and another important raw material, methyl tert-butyl ether (MTBE), and produces four types of products. The production scheme is shown as Figure 1. Both G90 and G93 belong to the gasoline blending pool GB, while the gasoline blending pool GB competes with the diesel blending pool DB and fuel oil FO in terms of product yields; thus, the product yields of G93, D, and FO are chosen as independent variables while the distribution of G90 can be derived from the joint probability distribution of the other three products. We consider a problem with the scale of Nt ) 5, Nk ) 3, and N1 ) 6. The name, type, identifier, and price of the materials and the initial inventory of products as well as the lower safety stock levels are shown in Table 5. Table 6 shows the mixing ratio of raw materials and production standard yields under each operation mode, the demand of products in each supply chain period, and a few other parameters. The vectors of value are in accordance with the order of raw material or products. Experiments under different confidence levels are conducted. The confidence level of each product is assumed to be the same which is Ri ) R, ∀i ∈ IFP. Figure 3 shows the corresponding optimal solution of objective function under different confidence levels. The second last dot in both Figure 3 and Figure 4 represents the case under a confidence level of R ) 0.99. According to Figure 3, profit will decrease with the increase of the confidence level of probability constrained inequality. This exactly explains that, under an aggressive strategy, with the increasing probability of constraints being violated, it is more likely to achieve higher profit. In contrast, under a conservative strategy, the high consistency of constraints is maintained at the expense of relatively lower profit. In Figure 3, it can be seen that before R ) 0.95 the trade-off between maximal profit and decision reliability is relatively smooth. However, it accelerates drastically after this point. Apparently, if the constraints of a stochastic process are completely satisfied, then the economic loss will be significant. From the comparison of R ) 0.99 and R ) 1, it can be seen that the probabilities of many bad scenarios are actually very low. However, they have a great effect in decreasing the expected profit. Therefore, decision makers can choose the value of R no more than 0.99

Table 6. Basic Configuration and Product Demands under Each Operation Mode parameter name raw material mixing ratio under operation mode 1 Xi,1 product yield ratio under operation mode 1 Yi,1 raw material mixing ratio under operation mode 2 Xi,2 product yield ratio under operation mode 2 Yi,2 raw material mixing ratio under operation mode 3 Xi,3 product yield ratio under operation mode 3 Yi,3 demand of product 1 in each supply chain period D1,t (kbbl) demand of product 2 in each supply chain period D2,t (kbbl) demand of product 3 in each supply chain period D3,t (kbbl) demand of product 4 in each supply chain period D4,t (kbbl) production capacity upper/lower bound in each scheduling period QU/QL (kbbl) operation mode changeover cost η (k$/times) operation cost for per processing amount of raw materials δ (k$/kbbl)

value [0.2, 0.3, 0.5, 0] [0.19, 0.38, 0.19, 0.18] [0.4, 0.4, 0.2, 0] [0.22, 0.34, 0.21, 0.17] [0.45, 0.2, 0.3, 0.05] [0.17, 0.36, 0.2, 0.21] [0, 1950, 0, 2250, 1200] [1900, 0, 3650, 0, 3500] [950, 0, 0, 2100, 2150] [0, 0, 2650, 0, 1750] 2000/1500 100 3

according to their own risk preference based on the experiment result of Figure 3. Figure 4 gives the processing amount Q under different confidence levels. With the increase of the confidence level, the processing amount also increases except the point of R ) 0.95. This is because when the processing amount increases, even if product yields decrease a certain range, the production volume can still be guaranteed, as well as the lower bound of safety stock levels. The operation modes within 15 compound periods are given in the vector after the colon as shown in eq 56. It should be noticed that the operation mode changeover strategy at the point of R ) 0.95 is different from others. Therefore, the exception of point R ) 0.95 is due to the adjustment of operation mode changeover.

Figure 3. Maximal profit under different confidence levels.

Figure 4. Processed crude in every scheduling period under different confidence levels.


285

R ) 0.95: [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1] R ) 1: [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2] R ) others: [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2] (56) In order to make sure that the lower bound of safety stock level is not violated, with the increase of confidence level, the optimization results not only tend to increase the production volume, but also tend to make more a conservative sales plan. This tendency is shown in Figure 5, where the vertical axis is the ratio of total product sale/demand within all periods. Except for the sharp decline at R ) 1, the tendency for R between 0.5 and 0.95 is relatively smooth, and the order fulfillment rate is also high. It is worth noting that this is just the sales plan in Figure 5. In fact, in later simulations, there are situations when the sales plan cannot be fulfilled as planned under low confidence level; even the safety stock level is designed to avoid the risk of production uncertainty. The inventories under different confidence levels are analyzed by simulation which is also the situation of constraints eq 28 being violated. For simplicity, only the inventory levels of product 1 (G93) under confidence levels R ) 0.5 and R ) 0.95 are presented. Figures 6 and 7 show the inventory evolution process over time through 100 times of simulations. Each simulation is a realization of one branch in the stochastic variable’s state tree, while decision variables in each simulation are the same due to the predetermined confidence level. It should be noted that the lines between the points do not show how inventory is varying between them. We keep the lines only because it is easier to see the points that belong to the same path as well as the inventory trend in each simulation (the same treatment is conducted in Figure 8 for this reason). The dashed-dotted horizontal line in Figures 6 and 7 represents the lower bound of the safety stock level, which should be 50 according to the information in Table 5. In periods t ) 1 and t ) 3, since there is no demand and thus no sales, products are stocked, resulting in a relatively large inventory level. In this case, it will not pose a danger to the lower bound of safety stock. In periods t ) 2, t ) 4, and t ) 5, products are sold and

Figure 5. Ratio of planned overall sales/demand under different confidence levels.

Figure 6. Evolution process of product 1’s inventory at R ) 0.5 through 100 times of simulation.

Figure 7. Evolution process of product 1’s inventory at R ) 0.95 through 100 times of simulation.

Figure 8. Sales shortage of each period through 100 simulations at R ) 0.5. Table 7. Times of Product 1’s Safety Stock Level Being Violated and Sales Not Fulfilled As Planned through 100 Times of Simulation t)2 times of inventory lower bound being violated at R ) 0.5 times of inventory lower bound being violated at R ) 0.95 total amount of sales plan not fulfilled at R ) 0.5 total amount of sales plan not fulfilled at R ) 0.95

t)4

t)5

17

29

43

0

0

5

132.472

908.534

1371.637

0

0

3.072

inventory is close to the lower bound. It can be seen from the two figures that the lower bound of the safety stock level is violated in a large number of simulations at R ) 0.5, while violation rarely happens at R ) 0.95. The times of inventory lower bound being violated at periods 2, 4, and 5 in 100 simulations is shown in Table 7. It can be seen that the number of violations at R ) 0.5 is apparently larger than that at R ) 0.95. Moreover, it can be seen that violation becomes more frequent as t increases. This is because inventory is a cumulative variable and its fluctuation will be amplified as production proceeds. The results in Figures 5 and 6 show that, with the increase of confidence level, the optimization result tends to increase production volume and make a conservative sales plan in order to comply with the constraints of lower inventory bound. This leads to another question: for the decision-making of low confidence level, is it possible to fulfill the sales plan at high standard by making use of relatively low production volume? An in-depth analysis is conducted for the above experiments. In this simulation, if the inventory level at the end of one period is smaller than 0, then it means the sales plan cannot be completed on time. In addition, the absolute value of the negative part represents the amount that cannot be fulfilled on time, which is also the sales shortage. Sales shortage of product 1 at R ) 0.5 in each period through 100 times of simulations is shown as Figure 8. It is noted that there are many times that the sales plan has not been correctly fulfilled. Besides, the sales shortage of a single period in several simulations is over 100 kbbl. In contrast, the cumulative sales shortage at R ) 0.95 is

286


only 3.072 kbbl through 100 times of simulations. To make it clearer, the amount of sales shortage in each period throughout 100 times of simulations is shown in Table 7. It can be seen that to achieve a high sales volume with low production volume under low confidence level is simply an illusion. Often in this situation, the sales plan cannot be fulfilled as planned. 5. Conclusions This paper focuses on the decision-making of a refinery supply chain operation with consideration of production uncertainty in the situation of multiple products and multiple operation modes. Markov chain is applied to describe the dynamic characteristic of product yields with the changeover of operation modes. Chance constrained programming is employed to reveal the trade-off between economic interest and reliability of decision results. It can be seen from the case studies that the decision maker can obtain abundant information about profitability, the extent of constraints violation, and sales shortage under different confidence levels. Decision makers can make a more thoughtful decision based on their own risk preferences. It should be emphasized that the method proposed in this paper to solve complex chance constrained problems can be extended to other similar situations, where the probability distribution of stochastic variables in each period depends on some binary decision variables. Acknowledgment The authors acknowledge financial support from the National High Technology Research and Development Program of China (863 Program) (Nos. 2007AA040702 and 2007AA04Z191). Nomenclature CHm,t,k ) binary variable that denotes whether operation mode is changed to m in compound scheduling period 〈t,k〉 Di,t ) demand for product i in supply chain period t e ) brief representation that indicates a state transition path of state tree f ) frequency of a state occurrence during a time interval F ) probability distribution function Gi,t ) production volume of product i in supply chain period t H ) binary variable that denotes whether there is any operation mode changeover M ) big constant that is used for big-M method Nk ) number of scheduling periods that each supply chain period contains Np ) number of products Nsi ) state number of product i’s yield fluctuation Nt ) number of supply chain periods in the whole planning horizon pi,a,b ) one-step transition probability of product i’s yield fluctuation from state a to state b under the same operation mode pbi,a(1) ) probability that product i’s yield fluctuation locates at state a in the period when there is a change of operation modes Pi ) state transition matrix of product i’s yield fluctuation Pbi(1) ) initial state vector that denotes the probability distribution of product i’s yield fluctuation in the period when there is a change of operation modes Pr ) probability of occurrence for an event PTri,t,e ) probability of path e in product i’s state tree in supply chain period t Q ) processing amount of total raw materials in each scheduling period Qm,t,k ) processing amount of total raw materials under operation mode m in compound scheduling period 〈t,k〉

QL, QU ) lower and upper bounds of processing amount, respectively Ri,t ) purchasing amount of raw materials i in supply chain period t SLi,t ) sales volume of product i in supply chain period t STri,t,e ) overall amount that actual yield deviates from standard yield along path e Tri,t,e ) detailed representation of a path that enumerates all possible evolution processes of product i’s yield fluctuation states in supply chain period t ViL, ViU ) lower and upper bounds of product i’s inventory level, respectively Vi,t ) inventory level of product i during supply chain period t wi,t, zi,t ) brief representations of the corresponding formulations Xi,m ) mixing ratio of raw material i under operation mode m Yi,m ) standard yield of product i under operation mode m Yî,m,t,k ) actual yield of product i under operation mode m in compound scheduling period 〈t,k〉 Zm,t,k ) binary variable that denotes whether operation mode m is executed in compound scheduling period 〈t,k〉 Greek Symbols R ) confidence level γi ) price of material i δ ) operation cost for per processing amount of raw materials η ) changeover cost of operation modes λ ) one possible state of product i’s yield fluctuations which is also an element of Λi Λi ) set that contains all possible states of product i’s yield fluctuations ξi,t,k ) product i’s yield fluctuation in compound scheduling period 〈t,k〉 τ ) number of periods that an operation mode is executed continuously Ωi,t ) set that contains all unique cumulative state values of product i’s yield fluctuations in supply chain period t Subscripts and Superscripts i ) material j ) cumulative period k ) scheduling period L ) lower bound m ) operation mode n ) power index t ) supply chain period U ) upper bound

Literature Cited (1) Tang, C. S. Perspectives in supply chain risk management. Int. J. Prod. Econ. 2006, 103 (2), 451–488. (2) Subrahmanyam, S.; Peknyt, J. F.; Reklaitis, G. V. Design of Batch Chemical-Plants under Market Uncertainty. Ind. Eng. Chem. Res. 1994, 33 (11), 2688–2701. (3) Gupta, A.; Maranas, C. D. A two-stage modeling arid solution framework for multisite midterm planning under demand uncertainty. Ind. Eng. Chem. Res. 2000, 39 (10), 3799–3813. (4) Grossmann, I. E.; Sargent, R. W. H. Optimum Design of ChemicalPlants with Uncertain Parameters. AIChE J. 1978, 24 (6), 1021–1028. (5) Gupta, A.; Maranas, C. D.; McDonald, C. M. Mid-term supply chain planning under demand uncertainty: customer demand satisfaction and inventory management. Comput. Chem. Eng. 2000, 24 (12), 2613– 2621. (6) Tsiakis, P.; Shah, N.; Pantelides, C. C. Design of multi-echelon supply chain networks under demand uncertainty. Ind. Eng. Chem. Res. 2001, 40 (16), 3585–3604.

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010 (7) 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. (8) Jung, J. Y.; Blau, G.; Pekny, J. F.; Reklaitis, G. V.; Eversdyk, D. A simulation based optimization approach to supply chain management under demand uncertainty. Comput. Chem. Eng. 2004, 28 (10), 2087–2106. (9) Escudero, L. F.; Quintana, F. J.; Salmeron, J. CORO, a modeling and an algorithmic framework for oil supply, transformation and distribution optimization under uncertainty. Eur. J. Oper. Res. 1999, 114 (3), 638–656. (10) Pongsakdi, A.; Rangsunvigit, P.; Siemanond, K.; Bagajewicz, M. J. Financial risk management in the planning of refinery operations. Int. J. Prod. Econ. 2006, 103 (1), 64–86. (11) Lakkhanawat, H.; Bagajewicz, M. J. Financial risk management with product pricing in the planning of refinery operations. Ind. Eng. Chem. Res. 2008, 47 (17), 6622–6639. (12) Hsieh, S.; Chiang, C. C. Manufacturing-to-sale planning model for fuel oil production. Int. J. AdV. Manuf. Technol. 2001, 18 (4), 303– 311. (13) Li, J.; Karimi, I. A.; Srinivasan, R., Robust Scheduling of Crude Oil Operations under Demand and Ship Arrival Uncertainty. Presented at the AIChE Annual Meeting, San Francisco, CA, 2006. (14) 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. (15) 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. (16) Brooks, R. W.; van Walsem, F. D.; Drury, J. Choosing cutpoints to optimize product yields. Hydrocarbon Process. 1999, 78 (11), 53–60. (17) Li, W. K.; Hui, C. W.; Li, A. X. Integrating CDU, FCC and product blending models into refinery planning. Comput. Chem. Eng. 2005, 29 (9), 2010–2028. (18) Mendez, C. A.; Grossmann, I. E.; Harjunkoski, I.; Kabore, P. A simultaneous optimization approach for off-line blending and scheduling of oil-refinery operations. Comput. Chem. Eng. 2006, 30 (4), 614–634. (19) Grinstead, C. M.; Snell, J. L. Introduction to Probability, 2nd ed.; AMS Bookstore: Providence, RI, 1997. (20) van der Laan, E.; Salomon, M. Production planning and inventory control with remanufacturing and disposal. Eur. J. Oper. Res. 1997, 102 (2), 264–278. (21) Kiesmuller, G. P.; van der Laan, E. A. An inventory model with dependent product demands and returns. Int. J. Prod. Econ. 2001, 72 (1), 73–87. (22) David, F. P.; Morris, A. C. Performance characteristics of stochastic integrated production-distribution systems. Eur. J. Oper. Res. 1993, 68, 23– 48. (23) Haurie, A. Time scale decomposition in production planning for unreliable flexible manufacturing systems. Eur. J. Oper. Res. 1995, 82, 339– 358. (24) Cochran, J. K.; Murugan, A.; Krishnamurthy, V. Generic Markov models for availability estimation and failure characterization in petroleum refineries. Comput. Oper. Res. 2001, 28 (1), 1–12.

287

(25) Kall, P.; Wallace, S. W. Stochastic programming; Wiley: New York, 1994. (26) Li, P.; Arellano-Garcia, H.; Wozny, G. Chance constrained programming approach to process optimization under uncertainty. Comput. Chem. Eng. 2008, 32 (1-2), 25–45. (27) Fu, M. C. Optimization for simulation: Theory vs. practice. INFORMS J. Comput. 2002, 14 (3), 192–215. (28) Liu, B. D.; Iwamura, K. Chance constrained programming with fuzzy parameters. Fuzzy Sets Syst. 1998, 94 (2), 227–237. (29) Charnes, A.; Cooper, W. W. Chance-constrained programming. Manage. Sci. 1959, 6 (1), 73–79. (30) Miller, B. L.; Wagner, H. M. Chance Constrained Programming with Joint Constraints. Oper. Res. 1965, 13 (6), 930–945. (31) Schwarm, A. T.; Nikolaou, M. Chance-constrained model predictive control. AIChE J. 1999, 45 (8), 1743–1752. (32) Petkov, S. B.; Maranas, C. D. Quantitative assessment of uncertainty in the optimization of metabolic pathways. Biotechnol. Bioeng. 1997, 56 (2), 145–161. (33) Henrion, R.; Moller, A. Optimization of a continuous distillation process under random inflow rate. Comput. Math. Appl. 2003, 45 (1-3), 247–262. (34) Li, P.; Wendt, M.; Arellano-Garcia, H.; Wozny, G. Optimal operation of distillation processes under uncertain inflows accumulated in a feed tank. AIChE J. 2002, 48 (6), 1198–1211. (35) Li, P.; Wendt, M.; Wozny, G. A probabilistically constrained model predictive controller. Automatica 2002, 38 (7), 1171–1176. (36) Li, P.; Wendt, M.; Wozny, G. Optimal production planning for chemical processes under uncertain market conditions. Chem. Eng. Technol. 2004, 27 (6), 641–651. (37) Li, W. K.; Hui, C. W.; Li, P.; Li, A. X. Refinery planning under uncertainty. Ind. Eng. Chem. Res. 2004, 43 (21), 6742–6755. (38) Pitty, S. S.; Li, W. K.; Adhitya, A.; Srinivasan, R.; Karimi, I. A. Decision support for integrated refinery supply chains. Part 1. Dynamic simulation. Comput. Chem. Eng. 2008, 32 (11), 2767–2786. (39) Montgomery, D. C.; Runger, G. C. Applied Statistics and Probability for Engineers, 3rd ed.; Wiley: New York, 2002. (40) 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. (41) Karimi, I. A.; McDonald, C. M. Planning and scheduling of parallel semicontinuous processes. 2. Short-term scheduling. Ind. Eng. Chem. Res. 1997, 36 (7), 2701–2714. (42) Hlavacek, I.; Chleboun, J.; Babuska, I. Uncertain input data problems and the worst scenario method, 1st ed.; Elsevier:: New York, 2004. (43) Nievergelt, Y. Foundations of logic and mathematics: applications to computer science and cryptography; Springer: New York, 2002.

ReceiVed for reView June 16, 2009 ReVised manuscript receiVed October 8, 2009 Accepted October 28, 2009 IE900968X

Supply Chain Optimization for Refinery with Considerations of

Recommend Documents