Refinery Planning under Uncertainty - American Chemical Society

6742

Ind. Eng. Chem. Res. 2004, 43, 6742-6755

Refinery Planning under Uncertainty Wenkai Li and Chi-Wai Hui* Chemical Engineering Department, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong

Pu Li Institut fu¨ r Prozess und Anlagentechnik, KWT 9, Technische Universita¨ t Berlin, 10623 Berlin, Germany

An-Xue Li Daqing Refining & Chemical Company, PetroChina Company Limited, Maanshan, Ranghulu District, Daqing City, Heilongjiang Province, China 163411

This paper presents a novel approach to refinery planning under uncertainty. To calculate the expectation of plant revenues, which is the main difficulty of the problem, loss functions are derived and applied. Different standard loss function approximation methods are compared and integrated into the planning model. The results show that the piecewise-linear approximation of the loss function can obtain good accuracy with improved solution speed. A general formulation for revenue and cost calculations is proposed by considering uncertainty in both the raw material and the product demand. To handle possible unmet customer demands, the hard-to-specify penalty functions of the two-stage programming are avoided and replaced by two of the decision maker’s service objectives, namely, confidence level and fill rate. Confidence level, which is the probability of satisfying customer demands, is commonly used in chance-constrained programming. However, fill rate, which is the proportion of demands that are met by a plant, is a greater concern of most managers. In this paper, fill rates are effectively calculated using the loss function. The maximum plant profit that satisfies certain fill rate objectives can then be obtained. Case studies show that a planning strategy that satisfies certain confidence level objectives might be overly lenient compared to a strategy that satisfies a fill rate objective. Because very few refinery planning models consider the influences of uncertainty, case studies including realworld large-scale refinery planning problems are used to illustrate the effectiveness of the proposed approach. 1. Introduction Most of the current plant planning/scheduling models are based on deterministic programming.1 However, because of volatile raw material prices, fluctuating products demands, and other changing market conditions, many parameters in a planning/scheduling model are uncertain. Upon realization of uncertainty, a schedule built on a deterministic approach will be nonrobust or in some cases even infeasible.2,3 For example, failure to account for significant demand fluctuations could lead to either unsatisfied customer demand, translating into loss of market share, or excessively high inventory holding costs.4 Even though almost 50 years have passed since the seminal works on uncertainty appeared,5,6 Dantzig still considers planning under uncertainty as one of the most important open problems in optimization.7 Thus, the consideration of uncertainty is critical in planning/scheduling models and is now attracting the attention of increasing numbers of researchers. To date, few researchers have studied refinery planning under uncertainty. Yet, refineries are vital components of national economies, and the fluctuations in the prices and demands of crude oil, gasoline, and diesel * To whom correspondence should be addressed. E-mail: [email protected]

oil are highly uncertain in reality, arising from uncertain global and national economic situations and indeterminate factors such as outbreaks of war, strikes, and diseases. Clay et al.2 used refinery planning in their case studies to demonstrate their solution algorithm. Neiro et al.8 performed supply chain optimization of refineries with the consideration of uncertainty using a scenariobased approach. 2. Solution Approaches to Problems with Uncertainties Stochastic programming deals with problems in which some parameters incorporated into the objective or constraints are uncertain. These uncertain parameters are usually described by probability distributions9,10 or by possible scenarios2,11,12 in stochastic programming. Stochastic programming mainly consists of recourse models2,11 and chance-constrained programming,9,10 distinguished by the methods used to describe the uncertain parameters and the algorithms used to solve the model. These modeling techniques are discussed in the next sections. 2.1. Recourse Models. Many research studies have employed recourse models for design/planning/scheduling under uncertainty. Recourse models use corrective actions (usually penalty functions) to compensate for the violation of constraints arising during the realization

10.1021/ie049737d CCC: $27.50 © 2004 American Chemical Society Published on Web 09/10/2004

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004 6743

Figure 1. Simple example.

of uncertainty. The two-stage model is one of the main paradigms of recourse models.13 The two-stage model divides the decision variables into two stages. The firststage variables are those that have to be decided right away, before any future realization of uncertain parameters. Then, the second-stage variables are those used as corrective measures or as a recourse against any infeasibilities arising during the realization of the uncertainty. In a design problem, the first-stage variables include capital investment, unit expansion, etc. After the unit capacities and the capital investments have been determined, further design or operational improvements can be achieved by selecting the values of the second-stage variables. This corrective action is known as recourse. In a planning problem, the firststage variables are the production rates required, the amounts of raw materials needed, etc. The second stage is to determine the amount of the product from other producers required to meet the actually realized product demand or the amount of raw material required from other suppliers to meet production requirements. Consider the simple example illustrated in Figure 1. Assume that the plant processes a single raw material to produce one product. We wish to determine the production rate, P, under uncertain customer demand, D. As shown in Figure 1, one unit of raw material can produce two units of product. The price of the raw material is $1/unit. The operating cost of each unit of product is $0.25. There are three demand scenarios: (1) D ) 10 units with probability 0.3, (2) D ) 20 units with probability 0.5, and (3) D ) 30 units with probability 0.2. Any unmet demand for the product is penalized at a cost of $3/unit. In terms of the two-stage model, in this example, the first-stage variable or the decision variable is the production rate, P, and the second-stage variable is the unmet demand. Two-stage programming divides the total cost into two categories: the first-stage cost, which includes the raw material cost and the operating cost, and the second-stage cost, which is the penalty cost for unmet demand. The objective of two-stage programming is then to find the optimal production rate, P, that will minimize the total cost. The two-stage model is formulated as follows 3

min cost ) 0.5P + 0.25P + s.t.

s

3PrsYs ∑ s)1

min cΤx s.t.

Ax e b xg0

(2)

where x ∈ Rn and c ∈ Rn are vectors and b ∈ Rm and A ∈ Rm×n. Assume that c and A are deterministic parameters and b is an uncertain parameter with the cumulative distribution function Φ. Then, assign a confidence level vector, r ∈ Rn, to reformulate constraint 2 as m

Aijxi e bi} g Ri, ∑ j)1

Pr{

i ) 1, 2, ..., n

(3)

In constraint 3, Pr is the operator of the probability computation. By applying the cumulative distribution function of b, constraint 3 can further be reformulated to m

Aijxi eΦi-1(1 - Ri), ∑ j)1

s

Y gD -P P g 0, Ys g 0

2.2. Chance-Constrained Programming. Some recourse models seek a solution to satisfy all scenarios; the solution is then too conservative. Some other recourse models allow the violation of soft second-stage constraints by including a penalty term (for example, the penalty of $3 for unmet demand in the simple example in Figure 1) in the objective function. That is, infeasibilities are allowed in the second stage provided that a certain penalty is imposed. However, the exact values of the penalty terms are difficult to determine because they include intangible components such as loss of goodwill and the costs of off-specification products or outsourcing of production requirements. Thus, in many cases of process operations, such a penalty term is not available. This difficulty is overcome in chance-constrained programming. Chance-constrained programming (or the probabilistic approach) seeks to satisfy the constraints at a predetermined confidence level using the known probability density/cumulative distribution of random variables.9 That is, rather than requiring that constraints containing the uncertain parameters always be satisfied or imposing penalties for infeasibilities, a probability of constraint satisfaction (usually called the confidence level) can be specified by the decision maker. This approach provides comprehensive information on economic achievements as a function of the desired confidence level. For complex plant operations under multiple uncertainties, the sources of risk that have the most significant impacts on the profitability can be identified. The basic procedure through which chance-constrained programming handles uncertainty is illustrated using the following simple linear model

(1)

where P is the production rate, Prs is the probability of demand Ds in scenario s, and Ys is the unmet demand in scenario s. The cost of the first stage, 0.5P + 0.25P, is the sum of the raw material cost and the operating 3 3PrsYs. cost. The expected second-stage cost is ∑s)1 When the above model is solved, the optimal production rate, P, turns out to be 20, with a corresponding total cost of $21.

i ) 1, 2, ..., n

(4)

In constraint 4, Ri and Φi-1 are known, and thus the right-hand side of constraint 4 is known and constraint 3 is reduced to deterministic linear constraints. Using constraint 4, the preceding probabilistic model is then converted to an ordinary deterministic linear programming model. In constraint 3, there is an assigned confidence level, Ri, for each bi. This implies that each bi can be satisfied at a different confidence level. For problems in which the satisfaction of all constraints as a whole is empha-

6744

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004

sized, all bi’s are assigned the same confidence level, λ. Constraint 3 is then replaced by a joint probabilistic constraint m

Pr{

Aijxi e bi} g λ ∑ j)1

(5)

When the uncertain bi’s are independent random variables, the joint probabilistic constraint can be decomposed into the multiplication of the single probabilities n

∏ i)1

m

[1 - Φi(

Aijxi)] g λ ∑ j)1

(6)

Nonlinear constraint 6 transforms the model into a nonlinear programming model. If the uncertain bi’s are correlated, simultaneous integration of multivariate probability distributions is required. Some methods, such as Gaussian quadrature and Monte Carlo sampling, can be used to approximate the multivariate integration. Again consider the example illustrated in Figure 1. Because the penalty for unmet demand is difficult to know in reality, we do not apply a penalty term for unmet demand in the objective function. Instead, we require that the demand be satisfied under a certain probability. Assuming that the demand conforms to a normal distribution with mean 19 and standard deviation 2, the example is reformulated by chance-constrained programming as follows

min cost ) 0.5P + 0.25P s.t.

Ρr{P g ζ} g R Pg0

(7)

where ζ represents the random demand and R is the confidence level. The deterministic equivalent of problem 7 is then

min cost ) 0.5P + 0.25P s.t.

P g Φ-1(R) Pg0

(8)

If we set the confidence level, R, to 0.95, then Φ-1(R) ) 22.29. After problem 8 is solved, the production rate, P, is 22.29, and the total cost is $16.7. 3. Applications of Modeling with Uncertainty Problems that include uncertainty mainly focus on plant design, plant planning/scheduling, and the supply chain. (a) Plant Design. For plant design problems, twostage programming is commonly applied in the literature. Wellons and Reklaitis14 investigated the design of multiproduct batch plants under uncertainty. They suggested a distinction between “hard” and “soft” constraints and introduced penalty terms for the latter. Analytical expressions for the expected profit objective were developed. Liu et al.11 studied the design of a chemical processing plant using two-stage (operating and design stages) programming. The first stage determines the investment decisions (e.g., equipment size). In the second stage, the investment decisions are made by considering all possible different future scenarios

comprehensively so that the total cost is minimized as a whole. Rooney and Biegler15 studied the optimal processing plant design that incorporated two types of uncertain parameters, unknown parameters and variable parameters. The entire domain of unknown parameters should always be feasible for a process, whereas control variables can be used to compensate for variable parameters. An extended two-stage programming model was proposed to address two types of uncertain parameters in their paper. (b) Plant Planning. Planning is essential for plants after the optimal design is determined. Such planning addresses applications such as feedstock selection and disposition, as well as overall material balance and conversion optimization.2 Clay et al.2 studied production planning using linear two-stage programming. Uncertain parameters were presented by finite discrete probability distribution functions. Refinery planning examples were used to demonstrate their partitioning algorithm. Lee et al.16 proposed a general strategy for treating open-shop batch process planning in which inventory was not considered. A discrete demand pattern was used, and a hybridization of the Monte Carlo simulation and simulated annealing techniques was utilized in their flexible planning algorithm. They claimed that it is possible to use only one-half the maximum capacity of a batch process because of uncertainty. They also concluded that the open-shop mode is preferred in a batch process when the inventory cost is large. As mentioned by some researchers,10 the penalty terms used in two-stage programming are difficult to specify. Moreover, two-stage programming provides no obvious information on the relation between reliability and profitability, which is crucial for decision makers. They proposed the use of chance-constrained programming in plant planning problems. Chance-constrained programming has been applied in many disciplines such as finance and management.17 However, very few applications have been made to chemical processes.18 In this paper, refinery planning is studied by applying chance-constrained programming. The tradeoff between the plant service level objectives, which take into account the confidence level and the fill rate, and overall plant profitability is captured. (c) Planning/Scheduling of the Supply Chain. The planning/scheduling of the supply chain under uncertainty is important in light of the ever-changing market conditions. Because there exist considerable lead times (the difference between the request and fulfilment times of an order) in the supply chain, the production variables have to be determined prior to the demand realization. The challenge is how to set optimal operating strategies to maximize the profit with satisfactory customer service levels. Gupta et al.19,20 studied supply chain planning using two-stage programming. Inventories were considered in the model, and penalty terms were used for stockout or inventory levels that are too low. In their approach, the penalty terms (such as a safety stock violation penalty) were difficult to estimate. Applequist et al.21 introduced a new metric for evaluating supply chain design and planning risk under uncertainty. A rational balance between the return and risk can thus be obtained. 4. Different Methods of Revenue Calculation Most commonly, the dominant uncertain parameter in chemical plant planning problems under uncertainty


is the product demand.21 The revenue of a plant is computed by the expected value of the operating variables whose values are adjusted according to the realization of uncertain parameters such as demand. Evaluating uncertain functions and their expectations generates the main difficulty in stochastic programming,22 and approaches in the literature differ primarily in how these expectations are computed.21 Clay et al.2 applied the certainty equivalent transformation (CET) to yield a deterministic equivalent problem. Ierapetritou et al.12 used the Gaussian quadrature formula to approximate the expected revenue. Liu et al.11 used Monte Carlo sampling to estimate the expectation of the objective function. Hene et al.23 used a specialized cubature integration formula that outperforms the Gaussian quadrature formula. They claimed that, when four uncertain parameters were used in the model, a fifth-degree Gaussian quadrature formula requires 625 integration points whereas their specialized cubature integration formula requires only 24 points to gain similar accuracies. However, all of these methods are somewhat complicated to use, and their solution speeds might not be satisfactory. In chance-constrained programming, when the uncertainty is represented by continuous distributions, complicated multiple integration techniques must be used.24,25 Thus, a linear approximation of these integrals would be useful in an optimization problem.21 Some works13 have employed approximation methods that used some functions to approximate the recourse function of a two-stage model. In this paper, different methods for calculating the expectations are discussed, and the expectations are approximated using simple piecewise-linear functions, which significantly reduce the computation time. In this paper, the objective of planning under uncertainty is to determine the current operating strategy, such as the unit operating load, the amount of raw material required, and the product production rate, by considering the uncertain product demand and raw material supply. In general, planning problems under uncertainty that have appeared in the literature can be modeled in following forms: (i) Minimize cost

min TC ) Crw + Cop + Cin + Civ + other costs

(9)

where TC is the total cost, Crw is the raw material cost, Cop is the operating cost, Cin is the inventory cost, and Civ is the investment cost. The objective is to minimize the total cost. Other costs, such as the transportation cost in the supply chain planning, should also be included. In eq 9, the plant revenue is not considered and is thus avoided in the model. (ii) Maximize profit I11

max revenue I - Crw - Cop - Cin - Civ - other costs (10) In eq 10, revenue I is calculated as

revenue I ) E[

∑j CjSj]

(11)

where Cj is the product price and Sj is the amount of product j produced by the plant and sold in the market. Underlying the above calculation is the assumption that

Table 1. Revenue Calculation of the Simple Example (Figure 1) demand

probability

product sold

revenue ($)

0 3 6 9 12 15 18 21 24 27 30 33

1.512 × 10-20 7.578 × 10-15 4.004 × 10-10 2.23 × 10-6 0.001309 0.0809864 0.528098 0.3629561 0.0262925 0.0002007 1.615 × 10-7 1.37 × 10-11

0 3 6 9 12 15 18 20 20 20 20 20

0.0 0.0 0.0 0.0 1.6 121.5 950.6 725.9 52.6 0.4 0.0 0.0 1852.5

the products produced by the plant are always absorbed by the market or that the amount of product j produced under each scenario is less than the market demand. However, this is not always the case in the real world. In many cases, the amount of product produced is larger than the market demand, and hence, the plant’s revenue will be reduced because of unsold products. (iii) Maximize profit II

max revenue II - Crw - Cop - Cin - Civ other costs (12) In eq 12, revenue II is calculated as

∑j cjdj]

revenue II ) E[

(13)

where dj is the demand for product j. The revenue is calculated from the expectation of the product demand. Underlying eq 13, it is assumed that the amount of product j produced under each scenario is greater than the market demand. Again consider the example illustrated in Figure 1. To explain how revenue is calculated, we first assume that the production rate, P, is 20 tons/day according to the production plans, and we assume that the demand conforms to a normal distribution with mean 19 ton/ day and standard deviation 2 ton/day, i.e., d ∼ N(19,2). We further assume that the price of the product is $100/ ton, which is a deterministic parameter. No inventory is considered in this example. Overproduced products are assumed to be valueless. This issue will be further discussed in section 8. The revenue results for this example are listed in Table 1. Even though the plant plans to produce 20 tons of the product today, the actual demand for the product tomorrow is still uncertain. The actual demand is listed in the first column, with a range of (0, 33). The probability corresponding to each demand is listed in the second column, which is lumped from the density function of the normal distribution, N(19,2). The third column lists the actual amount of product sold to the customers. For example, in the first scenario (the first row), the demand for the product is 0. The actual product sold to the market is also 0 no matter how much product is produced. On the contrary, in the eighth row, the demand for the product is 21 tons, but the production rate of the product is still 20 tons/day. Thus, the actual amount of the product sold to the market is 20 tons/day. The last column lists the revenue with consideration of its corresponding probability (probability × actual amount sold).

6746


From Table 1, we can see that the expected revenue is $1852.5, which is neither the value calculated using revenue I () 100 × 20 ) 2000) nor the value calculated using revenue II () 100 × 19 ) 1900). This difference is also noted by Petkov et al.4

Then, by using L(‚) to replace the loss function in eq L4, the revenue can be calculated as

5. Revenue Calculation Using the Loss Function

Equation L6 is the formulation for calculating the revenue of a single-product plant whose demand conforms to a normal distribution. From eq L6, the actual amount of product sold to customers, S, is calculated as

5.1. Derivation of the Loss Function. From Table 1, it can be seen that the actual amount of product sold to customers is the minimum value of the production rate and the product demand. Thus, the expected revenue can be represented as (assume that the plant produces only one product for the sake of clarity)

(P -σ θ)] ) Ch [θ - σL(z)]

[

revenue ) C h θ - σL

S ) θ - σL(z)

∑C × min(P, x)] x)0

revenue ) C h

∫0∞min(P, x) F(x) dx

(L2)

That is, if x e P, the revenue is C h ∫P0 xF(x) dx; if x > P, ∞ the revenue is C h ∫P PF(x) dx. Thus eq L2 becomes

∫0PxF(x) dx + ∫P∞PF(x) dx]

revenue ) C h[

(L3)

Assuming that the mean of the demand is θ ) ∫∞0 xF(x) dx, then eq L3 becomes

revenue ) C h [θ )C h [θ -

∫P∞xF(x) dx + ∫P∞PF(x) dx] ∫P (x - P)F(x) dx] ∞

(L4)

The term ∫∞P (x - P)F(x) dx is usually called the loss function. In eq L4, F(x) can be the density function of any distribution such as the normal distribution, the Poisson distribution, etc. For discrete distributions, the ∞ (x - P)F(x). The loss function can be written as ∑x)P loss function can be expressed in various forms for different distributions.26 For normal distributions with mean θ and standard deviation σ, the loss function becomes

∫P∞(x - P)F(x) dx ) ∫P∞(x - P)x 1

2πσ

e-(x-θ) /2σ dx 2

2

Let t ) (x - θ)/σ. Then the loss function becomes ∞ 1 (σt + θ - P) ∫P∞(x - P)F(x) dx ) ∫(P-θ)/σ x ∞ ∫(P-θ)/σ (t - P -σ θ)x1

2π

) σL

e-t /2σ dt 2

2πσ

)σ

e-t /2 dt 2

(P -σ θ)

where L(‚) is the standard loss function. We define

z)

P-θ σ

J

(L1)

where C is the product price, P is the rate of production of the product, and x is the demand. Assuming that the density function of the uncertain demand is F(x) and replacing the price C with its expected value, C h , eq L1 can be reformulated as

(L5)

(L6a)

For a multiproduct plant, eq L6 can be written as

∞

revenue ) E[

(L6)

revenue )

Cj[θj - σjL(zj)] ∑ j)1

(L7)

where J is the number of products and C h j is the expected price of product j. The demand for product j conforms to a normal distribution with mean θj and standard deviation σj. Pj is the production rate of product j, and we have zj ) (Pj - θj)/σj. For other probability distributions, the loss function can be rewritten in different forms, and corresponding formulations can be derived from eq L4. In this paper, we focus on the normal distribution because it captures the essential features of demand uncertainty. Petkov et al.4 developed a formula for the calculation of the revenue (index t is eliminated in a single-period planning problem)

revenue ) J

Cj{θj + σj[-φ(Kj) + (1 - Φ(Kj))Kj]} ∑ j)1

(L7a)

where Kj ) (Pj - θj)/σj and φ(‚) is the standard normal density function. Because functions φ(‚) and Φ(‚) are included, eq L7a is a highly nonlinear equation. In fact, the term -φ(Kj) + [1 - Φ(Kj)]Kj in eq L7a equals the negative of L(Kj). The loss function, L(‚), is much easier to handle in the planning problem shown later in this paper. If multiperiod and multiplant scheduling is necessary, eq L6 can be rewritten as J

revenue )

T

[

∑ ∑ Cj,t θj,t - σj,tL j)1 t)1

(

)]

Pj,t - θj,t σj,t

(L7b)

where t ) 1, ..., T, represents the time periods; C h j,t is the expected price of product j in time period t; and Pj,t is the production rate of product j in time period t. The demand for product j in time period t conforms to a normal distribution with mean θj,t and standard deviation σj,t. In eq L7, C h j, θj, and σj are known parameters. To apply eq L7 in a planning model, the only difficulty lies in how to express the standard loss function, L(‚). Standard tables are available for the values of L(‚).27 However, these tables are difficult to use in planning models. In this paper, two methods (described in sections 5.4 and 5.5) are used to approximate L(‚) after some discussions of the loss function in sections 5.2 and 5.3. Comparisons of different approximation methods are also given. 5.2. Properties of the Loss Function. The physical meaning of a loss function can be considered as the


Figure 4. Loss function of the exponential distribution.

Figure 2. Standard loss function.

Figure 5. Loss function of the Poisson distribution (mean ) 10). Figure 3. Exponential distribution with λ ) 0.01 and θ ) 300.

amount of unmet demand (the backorder level) of a plant facing uncertain demand.26 The loss function, ∫∞P (x - P)F(x) dx, is a convex function of P. This can be seen from the nondecreasing first-order derivative [) -∫∞P F(x) dx] of the loss function of P and the nonnegative second-order derivative [) F(P)] of the loss function of P. If F(x) is the density function of a normal distribution, N(θ,σ), the value of the loss function can be calculated by the standard loss function, i.e.

∫P∞(x - P)F(x) dx ) σL(P -σ θ) ) σL(z) This is an important property because we can confine our study to the characteristics of the standard loss function, L(‚). The standard loss function, L(‚), is shown in Figure 2. As z goes to infinity, the value of L(z) goes to zero

lim L(z) ) 0

zf+∞

As z goes to minus infinity, the value of L(z) goes to -z

lim L(z) ) -z

zf-∞

5.3. Loss Functions of Other Distributions. In eq L4, F(x) can be the density function of other distributions. For an exponential distribution (Figure 3), the loss function can be represented as (see Appendix I for the derivation) λ(θ-P)

∫P∞(x - P)F(x) dx ) e

λ

where λ is the lambda of the exponential distribution, θ is the location parameter, and P is the production rate. Figure 4 shows the loss function of the exponential distribution with λ ) 0.01 and θ ) 300. For discrete distributions such as the Poisson distribution, the loss function can be written as26 ∞

∑ (x - P)F(x)

or θF(P) + {(θ - P)[1 - G(P)]}

x)P

where θ is the mean of the Poisson distribution and F(P) and G(P) are the probability mass function and cumulative distribution function, respectively. P is the production rate. Figure 5 shows the loss function of the Poisson distribution with mean 10. 5.4. Nonlinear Approximation of the Standard Loss Function. A straightforward method for approximating L(‚) is to use one polynomial function. However, in Figure 2, it can be seen that the discrepancy between the standard loss function and the polynomial approximation function (approximation ranges from -10 to 50) is significant. Thus, it is not appropriate to approximate L(‚) using a single polynomial function. Multipolynomial functions are used in this paper to approximate L(‚) accurately. The whole range was divided into four parts, i.e., (-∞, -3], (-3, 0], (0, 3], and (3, +∞). Four polynomial functions [LN1(z)-LN4(z)] are used for these parts. The coefficients of the polynomial functions obtained by curve fitting of the standard loss function using the LSM (least-squares method) are listed in Table 2. Sixth-order polynomial functions in the form of eq L8 are used because they have high accuracy with low order

LNi(z) ) ai + biz + ciz2 + diz3 + eiz4 + fiz5 + giz6, i ) 1-4 (L8) where ai-gi are the coefficients listed in Table 2 and z is calculated using eq L5. Because L(z) is very close to

6748

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004 Table 3. Coefficients of the Ten Linear Functions LL1 (-∞, -3] LL2 (-3, -2] LL3 (-2, -1] LL4 (-1, -0.5] LL5 (-0.5, 0] LL6 (0, 0.5] LL7 (0.5, 1] LL8 (1, 2] LL9 (2, 3] LL10 (3, +∞)

Figure 6. Nonlinear approximation functions of the standard loss function.

ai bi ci di ei fi gi sum of the squared error

0 -1 0 0 0 0 0 N/A

LN2 (-3, 0]

LN3 (0, 3)

0.39893 0.398923956 -0.5004 -0.49952273 0.19653 0.196346393 -0.0084 0.008483501 -0.028 -0.027947893 -0.0078 0.007732298 -0.0007 -0.000701398 2.71 × 10-9 5.01174 × 10-8

LN4 (3, +∞) 0 0 0 0 0 0 0 N/A

-z when z < -3, L(z) ) -z is used in this range. When z > 3, L(z) is a very small value (less than 4 × 10-4), L(z) ) 0 is used in this range. For ranges (-3, 0) and (0, 3), two-sixth-order polynomial functions are used. These two functions together with the standard loss function L(‚) are shown in Figure 6. For the revenue maximizing problem, the value of L(z) is to be minimized in most cases (further discussions are given in section 8). The value of L(z) can thus be modeled as

L(z) g LNi(z), i ) 1-4

(L9)

5.5. Linear Piecewise Approximation of a Standard Loss Function. Even though nonlinear functions approximate the standard loss function with high accuracy, they transform a linear planning model into a nonlinear one. As a result, the solution speed could be reduced. Several piecewise-linear functions are used here as an alternative to the approximation of the

bi

0 0.024708634 0.158140087 0.312277656 0.39894228 0.39894228 0.312277656 0.158140087 0.024708634 0

-1 -0.991891104 -0.925175378 -0.771037809 -0.59770856 -0.40229144 -0.228962191 -0.074824622 -0.008108896 0

standard loss function. The whole range of the independent variable, z, is divided into 10 parts. Linear functions (LL1-LL10) are used for these parts. The coefficients of the linear functions are listed in Table 3. Each linear function has the following form

Table 2. Coefficients of the Four Polynomial Functions LN1 (-∞, -3]

ai

LLi(z) ) ai + biz, i ) 1-10

(L10)

where ai and bi are coefficients obtained using the LSM. Figure 7 illustrates the four linear functions in the range (0, 3). Similarly, the value of the standard loss function can be calculated as

L(z) g LLi(z), i ) 1-10

(L11)

The approximation methods are suitable for demands with normal distributions. For other distributions, we can use methods similar to the above to approximate their loss functions such as the loss function of the exponential distribution shown in Figure 4 and the loss function of the Poisson distribution shown in Figure 5. The difference is that one needs to perform the approximation once the parameter of the distribution changes. 6. Raw Material and Utility Cost Calculations The inflows of a plant, such as raw materials and utilities, are from different suppliers. The supply of these inflows can also be uncertain. For example, the supply of crude oil to a refinery from oil fields might change as a result of the fluctuating demands of other refineries and the changing production rate of the oil fields. The calculations of raw material costs and utility costs are done by methods similar to that proposed above if the supply of the raw material or utility is uncertain. The following formulations are used to

Figure 7. Linear approximation functions of the standard loss function.


calculate the costs of uncertain raw materials and utilities

raw material cost ) M

∑ Ch m[θm - ∫R (x - Rm)Fm(x) dx] m)1 ∞

m

L

C h l[θl - ∫U (x - Ul)Fl(x) dx] ∑ l)1

utility cost )

∞

(L12) Figure 8. Plant with inventory.

(L13)

l

where M is the number of uncertain raw materials, C hm is the expected price of raw material m, and Rm is the consumption of raw material m. θm is the mean of the uncertain raw material m, and Fm(x) is its density distribution function. L is the number of uncertain utilities, C h l is the expected price of utility l, and Ul is the consumption of utility l. θl is the mean of the uncertain utility l, and Fl(x) is its density distribution function. If the uncertain raw materials and utilities conform to discrete distributions, then the integrals in eqs L12 and L13 have to be replaced by their corresponding summation forms. If the uncertain raw materials and utilities conform to normal distributions, eqs L12 and L13 can be simplified to eqs L14 and L15, respectively

raw material cost ) M

(

[

∑ Ch m θm - σmL m)1 L

utility cost )

∑ l)1

σm

[ ( )]

C h 1 θl - σlL

)]

Rm - θm

(L14)

Ul - θl

(L15)

σl

7. Calculation of the Total Revenue and Raw Material/Utility Cost Aside from the uncertain inflows, the supplies of some inflows to a plant are considered to be stable and unlimited. Similarly, the market demands for some products of a plant can also be considered as unlimited compared to the production rate of this product. Thus, if only these certain inflows and products are taken into account, the plant can obtain as much of any raw material or utility as it wants and sell as much of any product as it produces. The revenue of a plant is the sum of the revenues of uncertain products and certain products. Similarly, the raw material/utilities cost of a plant is the sum of uncertain raw material/utilities costs and certain raw material/utilities costs. The following objective function gives the sum of these revenues and costs (a normal distribution is assumed for all uncertain inflows and products) J

max

[ ( )] ∑ [ ( )] ∑ ∑ [ ( )] ∑

Cj θj - σjL ∑ j)1

Pj - θj

∑ Ch hPh h)1

C h m θm - σmL

Rm - θ m

L

-

C h l θl - σlL

l)1

- other costs

N

-

σm

m)1

Ul - θl σl

8.1. Loss Functions without Inventory. In a plant operation, two operating modes, the closed-shop mode and the open-shop mode, are generally used.16 In the closed-shop mode, all of the customers’ requests are fulfilled by inventory, and production tasks are generally the result of inventory replenishment decisions. In the open-shop mode, products are delivered to customers directly, and no inventory is stocked. Examples of the open-shop mode, such as the news vendor model26 in which the sales of holiday lights disappear after Christmas, for instance, can be found in the real world. The closed-shop mode has received extensive attention, whereas only a few researchers have yet studied the features of the open-shop operating mode.16 The major advantage of the open-shop mode is that it has a much lower inventory cost than the closed-shop mode has. For cases in which the inventory cost is large, the open-shop mode is preferred.16 The simple example illustrated in Figure 1 assumes that an open-shop mode is used at the plant. For an open-shop mode plant, the methods described in section 5 can be used for revenue calculations. 8.2. Loss Functions with Inventory in a Single Period. A plant with inventory is operated in the closed-shop mode. Figure 8 shows the configuration of a typical plant with inventory. Each product j is allocated a product tank with inventory Ij. Pj is the production rate of product j. The plant faces an uncertain customer demand dj. It can be seen in Figure 8 that, different from the open-shop mode, the maximum amount of product j that can be sold to customers (the deliverable amount, Dj) is the sum of the production rate, P, and the initial inventory, I0,j. In this case, the variable P in eqs L1-L6 should be replaced by the deliverable amount Dj () Pj + I0,j). Equation L5 should be replaced by

zj )

+ other revenues -

8. Loss Functions with and without Inventory

H

+

σj

M

where H is the number of certain products, C h h is the expected price of certain product h, and Ph is the production rate of product h; N is the number of certain raw materials, C h n is the expected price of certain raw material n, and Rn is the consumption of raw material n; K is the number of certain utilities, C h k is the expected price of certain utility k, and Uk is the consumption of utility k.

C h nR n

n)1

K

-

C h kUk

k)1

(L16)

Pj + I0,j - θj σj

(L17)

8.3. Feasible Region of Loss Value in the Multiperiod Case with Inventory. For a multiperiod case, a time interval index, t, should be added to eqs L1-L6. The deliverable amount of a product in a single-period case with inventory is the sum of the initial inventory and the production rate. Because the initial inventory is a fixed parameter, the deliverable amount, Dj, depends solely on the production rate, Pj. In the multiperiod case with inventory, the maximum deliverable

6750


9. Type 1 and Type 2 Service Levels

Figure 9. Feasible region of the loss value.

amount of a product in a time interval is the sum of the amount of production plus the amount of inventory available in that time interval. That is, the deliverable amount of product j at time t is

Dj,t ) Pj,t + Ij,t-1

(L18)

Equation L17 is then replaced by

zj,t )

Pj,t + Ij,t-1 - θj,t σj,t

(L19)

In the multiperiod case, the inventory at the end of the previous period, It-1, is variable. It is possible for a plant to adjust inventory levels or production rates among time intervals to maximize the profit. This can be realized by increasing the loss value obtained from eq L9 or L11. For normally distributed demand, the actual amount of product j sold to customers, Sj,t, during time interval t is calculated by

Sj,t ) θj,t - σj,tL(zj,t)

(L20)

where θj,t and σj,t are the mean and standard deviation, respectively, of the normally distributed demand in time interval t. If L(zj,t) increases, then Sj,t decreases. The inventory level, Ij,t, of product j is calculated as

Ij,t ) Dj,t - Sj,t ) Ij,t-1 + Pj,t - Sj,t

(L21)

Thus, decreasing Sj,t might subsequently increase Ij,t if Pj,t is unchanged. Because Pj,t is also a variable, another function of the increasing loss value is that, if the plant wants to produce less product j so that it can increase the production rates of other profitable products, the production rate of product j can be decreased by decreasing the amount of the product sold. From eq L21, if Ij,t is unchanged, then decreasing Sj,t can decrease Pj,t. This is reasonable in practice as a plant can choose not to sell its products during certain time intervals to increase the total profit even though the product demand is large. The fill rates or confidence levels during these time intervals can be set to lower values. Because the actual amount of product j sold is Sj,t g 0, then from eq L20, we have L(zj,t) e θj,t/σj,t. This represents the upper limit of the loss value of product j in time interval t. The striped area in Figure 9 shows the feasible region of the loss value, L(zj,t). The value calculated from the loss function is the lower bound of the loss value because it represents the maximum amount of product that can be sold. This is the primary reason that the linear approximation functions (see Figure 7) are located above L(‚).

9.1. Definitions of the Two Service Levels and the Difference between Them. There are two types of service levels (or customer satisfaction levels) that are commonly used in the processing industry. A simple example is given in Appendix II to show the difference between the two service levels. The type 1 service level is the probability of not stocking-out in all scenarios or horizons.27 Stocking-out occurs whenever a demand is not satisfied from a plant, regardless of the amount of unmet demand. The type 2 service level (also often called the fill rate) measures the proportion of demands that are met by a plant.27 A type 1 service level is suitable when a stocking-out occurrence has the same consequences regardless of its time or amount. However, the type 1 service level is not how service is interpreted in most applications.27 When a manager mentions that a certain percentage of service is provided, he usually means that a certain percentage of demand is filled by his plant. 9.2. Applications of the Two Service Levels. Up to now, the customer satisfaction level that has been used in chance-constrained programming is the type 1 service level. Many researchers10 have studied the application of chance-constrained programming using the type 1 service level (usually also called the confidence level). A constraint that represents the objective of a type 1 service level has the following form

Pr(Dj g dj) ) R where Dj is the production rate (no inventory) or the deliverable amount (with inventory) of product j, dj is the market demand for product j, and R is the type 1 service level or confidence level defined by the decision maker. The above constraint becomes the following when chance-constrained programming is applied10

Dj g Φ-1(R)

(L22)

In the above constraint, Φ-1 is the known reverse cumulative distribution function of the product demand. Thus, the right-hand side of the above constraint is a known value, and the above constraint is reduced to a deterministic constraint. Even though it is straightforward to apply a type 1 service objective in a planning model, the type 1 service level is not the one that most managers want. Furthermore, the type 1 service level cannot accurately approximate a type 2 service objective. Thus, calculating the type 2 service level remains an essential challenge in refinery planning. In this paper, the type 2 service level is calculated using a loss function. From its definition, the loss function measures the amount of unmet demand of a plant facing uncertain demand. Thus, for a single-period, single-product plant, the type 2 service level or fill rate can be calculated as

β)

S θ - σL(z) ) θ θ

where θ and σ are the mean and standard deviation, respectively, of the normally distributed demand. L(z) is the loss value that can be calculated using eq L9 or L11. S is the actual amount of product j sold to customers. z is calculated using eq L5. β is the type 2 service level or fill rate defined by the decision maker.


Figure 11. Configuration of the blending example. Table 4. Prices and Octane Numbers for the Illustrative Example raw materials

products

gasoline MTBE 90# gasoline 93# gasoline price (Yuan/ton) octane number

1400 70

3500 101

3215 90

3387 93

conforms to a normal distribution. A similar procedure with some modifications can also be used to deal with other types of distributions and the raw material/utility cost calculation.

Figure 10. Main flow diagram for using the loss function.

If a type 2 service objective has to be satisfied, then the following constraint is added

θ - σL(z) g βθ

(L23)

If linear functions are used to approximate the standard loss function from the production rate or the deliverable amount, the above constraint is linear and deterministic. The application of these two service levels is illustrated in detail in case study 2. 10. Main Flow Diagram for Using the Loss Function The main flow diagram for using the lose function in a refinery planning model is illustrated in Figure 10. The whole procedure consists of the following steps: (i) Obtain the input data. This includes the mean and standard deviation of the product demand. (ii) Determine the deliverable amount D. In a problem without inventory, D is the production rate of the product. In a single-period problem with inventory, D is the production rate plus the initial inventory. In a multiperiod problem with inventory, D is determined using eq L18. (iii) Calculate z. Equations L5, L17, and L19 are used for a problem without inventory, a single-period problem with inventory, and a multiperiod problem with inventory, respectively. (iv) Determine the loss value, L(z), using eq L9 or L11. (v) Determine the actual amount of product sold to customers, S. Equations L6a and L20 can be used for different problems. (vi) Add constraints for the objective of a type 1 service level using eq L22, or add constraints for the objective of a type 2 service level using eq L23. (vii) Calculate the revenues and costs of the plant using eq L16. Note that the above procedure deals with the uncertain product demand and assumes that the demand

11. Case Studies Three examples with different problem sizes and complexities demonstrate the capabilities and effectiveness of the new approaches and formulations proposed in this paper. The models were formulated with GAMS28 on a 933-MHz Pentium III PC. The solver OSL in GAMS 2.25 was used for MILP, and the solver DICOPT was used for MINLP. 11.1. Case 1: Blending Case. A simple example in Figure 11 is used to illustrate how the new approach is formulated. Figure 11 consists of one gasoline-blending (GB) unit, which are very common in refineries. MTBE and gasoline enter the GB to produce two products: 90# and 93# gasolines. The blending requirement is that the octane number of each product should be equal to or greater than the required octane number of that product. Table 4 lists the prices and octane numbers of the feed and product streams. The demands for 90# and 93# gasoline are assumed to conform to N(50,10) and N(70,10), respectively. The objective of this case study is to determine the optimal production plan prior to the realization of uncertain demand. No inventory is considered, and the overproduced products are assumed to be valueless. The mathematical model that applies the nonlinear approximation approach consists of constraints 1-12 and obj from the following list

101 × M90 + 70 × G90 g 90 × GASO90 101 × M93 + 70 × G93 g 93 × GASO93

(1) (2)

LossRt90 ) (GASO90 - DM_mean90)/DM_var90 LossRt93 ) (GASO93 - DM_mean93)/DM_var93 LossV90 g LNi(LossRt90), i ) 1-4

(3) (4) (5)

LossV93 g LNi(LossRt93), i ) 1-4

(6)

LossSL90 ) DM_mean90 - DM_var90 × LossV90 (7) LossSL93 ) DM_mean93 - DM_var93 × LossV93 (8)

MTBE ) M93 + M90 GASO ) G90 + G93

(9) (10)

GASO90 ) M90 + G90 GASO93 ) M93 + G93

(11) (12)

profit ) LossSL90 × 3215 + LossSL93 × 3387 GASO × 1400 - MTBE × 3500 (obj)

6752


Table 5. Comparison of the Three Methods in Case 1 method

profit (Yuan)

solution time (s)

scenario Petkova nonlinear approximation piecewise-linear approximation

38830.94 38733.69 38733.68 38642.09

0.6 0.05 0.04 0.03

a

Method used by Petkov et al.4

In the above model, MTBE and GASO signify the MTBE and crude gasoline used for blending the 90# and 93# gasoline products (GASO90 and GASO93, respectively). M90 and M93 represent the flow rates of MTBE that are needed for the blending of 90# and 93# gasolines, respectively. G90 and G93 represent the flow rates of crude gasoline that are needed for the blending for 90# and 93# gasolines, respectively. DM_mean90 and DM_mean93 are the means of the demands for GASO90 and GASO93, respectively. DM_var90 and DM_var93 are the standard deviations of the demands for GASO90 and GASO93, respectively. LossRt90 and LossRt93 are the standardized production rates of 90# and 93# gasolines, respectively. LossV90 and LossV93 are the values of the standard loss function for 90# and 93# gasolines, respectively. LossSL90 and LossSL93 are the actual amounts of GASO90 and GASO93, respectively, sold to customers. profit represents the total profit of the system. Constraints 1 and 2 ensure the octane numbers of the 90# and 93# gasolines, respectively, meet the product specifications. Constraints 5 and 6 obtain the values of the standard loss function using LossRt90 and LossRt93 calculated by constraints 3 and 4, respectively. The nonlinear approximation method is shown here. LNi(‚) in eqs 5 and 6 represents the polynomial approximation functions defined by eq L8. Constraints 7 and 8 are used to calculate the actual amounts of GASO90 and GASO93, respectively, sold to customers. Constraints 9-12 are for the mass balance of the blending unit. If the piecewise-linear approximation method is used, then constraints 5 and 6 should be replaced by eqs 5L and 6L, respectively

LossV90 g LLi(LossRt90), i ) 1-10

(5L)

LossV93 g LLi(LossRt93), i ) 1-10

(6L)

where LLi(‚) in eqs 5L and 6L represents the linear approximation functions defined by eq L10. For comparison, the scenario method is also used. The whole demand space is lumped into 36 scenarios. Their corresponding probabilities are also obtained from the normal distribution of the demands for GASO90 and GASO93. The results reported shown in Table 5. The Petkov method uses eq L7a to calculate the loss function; its results are thus accurate. Table 5 shows that the objective values obtained from the three approximation methods are very close (the largest difference among the three methods is less than 0.5%). Thus, the two approximation methods are suitable for solving a planning problem under uncertainty. The nonlinear approximation method obtains almost the same result as the Petkov method. The piecewise-linear approximation method underestimates the plant profit because the linear functions overestimate the value of the standard loss function (Figure 7). However, this inaccuracy can be overcome by using more linear functions to ap-

Figure 12. Overall configuration of the D refinery complex.

proximate the loss function. The solution time of the piecewise-linear approximation method is the best. In this very small example, except for the scenario method, the other three methods have similar solution results and solution times. However, this is not the case for large-scale problems, as shown in cases 2 and 3. 11.2. Case 2: Monthly Planning for a Real-World Refinery Complex. Case 2 is a real industrial example provided by a refinery complex located in P.R. China (referred to as the D refinery complex in this paper) with some modifications in data to maintain confidentiality. The D refinery complex can process up to 5.5 million tons of crude oil per year. It uses only one type of crude oil, which is bought from a large oil field nearby. It mainly consists of four plants, namely, the refinery plant, the lubricant oil plant, the polymerization plant, and the utility plant. Figure 12 shows the configuration of the D refinery complex. In Figure 12, crude oil is processed in the refinery plant to produce different intermediate or end products. The intermediate products are later sent to the blending units to produce end products such as gasoline and diesel oil. Propylene from the refinery plant is sent to the polymerization plant to produce different grades of polymer. Part of the heavy oil from the crude distillation unit (CDU) in the refinery plant is sent to the lubricant oil plant to produce different grades of lubricant oil. All three processing plants (the refinery plant, the polymerization plant, and the lubricant oil plant) are closely connected to the utility plant. The fuel oil and fuel gas from the refinery plant and the lubricant oil plant are sent to the utility plant as raw material. Steam and electricity generated in the utility plant are sent to the three processing plants as energy sources for the processing units. Electricity and steam can be imported from the market if needed. An integrated monthly planning model has been developed for all of these plants to consider the interactions among these plants. Models for the units in the four plants are included. For example, an effective linear model for determining CDU WTRs (weight transfer ratios) is applied.29 The CDU model has higher accuracy than traditional linear CDU models and higher solution speeds than rigorous simulation models.29 Octane number and pour point are used as the quality indexes of gasoline and diesel oil, respectively. The common linear blending rule is used for gasoline blending. Diesel oil blending is linearized by converting the pour points of intermediate diesel oils into pour point indexes.30,31 The on/off status of some processing units, such as the CDU and fluid catalytic cracking (FCC) unit in the refinery plant, boilers and turbines in the utility plant, etc., are determined by binary variables. In case 2, no inventory is considered. Products are sold to customers directly. The market demands for 90# gasoline and 93#

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004 6753 Table 6. Monthly Planning Results of Case 2 method Petkov nonlinear approximation piecewise-linear approximation Petkov nonlinear approximation piecewise-linear approximation Petkov nonlinear approximation piecewise-linear approximation

CLa FLc solution profit (106 Yuan) time (s) 90#b 93#b 90#b 93#b Case 2A (without CL) 154.0 1.9 0.57 0.59 0.971 0.971 154.0 2.4 0.57 0.59 0.971 0.971 153.8

0.4

0.50 0.69 0.964 0.980

Case 2B (with 0.95 CL) 134.7 1.7 0.95 0.95 0.998 0.998 134.7 2.1 0.95 0.95 0.998 0.998 134.4

0.9

0.95 0.95 0.997 0.996

Case 2C (with 0.980 FL) 152.0 2.2 0.72 0.74 0.985 0.984 153.4 1.9 0.67 0.69 0.980 0.980 153.3

0.4

0.67 0.69 0.980 0.980

a CL represents the confidence level or the type 1 service level. 90# represents 90# gasoline; 93# represents 93# gasoline. c FL represents the fill rate or the type 2 service level.

b

gasoline are assumed to conform to normal distributions N(55000,5000) and N(30000,3000), respectively. For the sake of simplicity, the market demands for all other products are assumed to be deterministic without loss of generality of the model. The objective of this case study is to determine the optimal monthly planning strategy for the D refinery complex under uncertain market conditions. In case 1, there is no information on how the customer demands are met. A high proportion of customer demands might not be satisfied, and a very low customer satisfaction level could be obtained from the planning results of case 1. Thus, the plant might suffer by losing future customers because of the low customer satisfaction level. This problem is handled in case 2. The two approximation methods and the Petkov method are applied for the monthly planning of case 2. Table 6 lists the results. The model involves 17 discrete variables, 399 single continuous variables, and 333 single constraints for the piecewise-linear approximation method. (The Petkov method involves a slightly different number of single constraints). In Table 6, case 2A shows the results for monthly planning without consideration of the customer satisfaction level. Cases 2B and 2C are the results with consideration of type 1 and type 2 customer satisfaction levels, respectively. The last four columns list the values of type 1 and type 2 customer satisfaction levels for 90# and 93# gasoline. It can be seen that the planning results without consideration of the customer satisfaction level obtained the highest profit. In fact, the production rates of 90# and 93# gasoline in case 2A (when the piecewise-linear approximation method is applied) are 55000.0 and 31500.0 ton/month, respectively. The corresponding type 1 customer satisfaction levels for 90# and 93# gasoline are 50% and 69%, respectively. These type 1 customer satisfaction levels are low. In case 2B, another constraint, which sets the type 1 customer satisfaction level for each product to 95%, was added to the planning model. The corresponding production rates of 90# and 93# gasoline (when the piecewise-linear approximation method is applied) are 63224.3 and 34934.6 ton/month, respectively. In this case, we increased the type 1 customer satisfaction level

to a high level (0.95). However, we also suffered from significantly decreased profit (12.4% lower). In case 2C, the constraint for the type 2 customer satisfaction level was added to replace the constraint for the type 1 customer satisfaction level. The type 2 customer satisfaction level was set to 0.98. We can see that the profit for case 2C is slightly decreased compared to that for case 2A. This is because, in case 2A, the type 2 customer satisfaction level is already near the required value (0.98). Because many managers use the fill rate or the type 2 customer satisfaction level as their criterion, it is important to calculate the actual value of the fill rate. In our case study, a very high type 2 customer satisfaction level (0.998 in case 2B) was achieved when the type 1 customer satisfaction level was set to 0.95. In general, the actual fill rate is higher than the actual type 1 customer satisfaction level. However, there is no direct method available that can be used to calculate the fill rate from the type 1 customer satisfaction level. If the managers do not require a very high fill rate and the planners set the type 1 customer satisfaction level to a certain value for which the final planning results achieve an unnecessarily high fill rate, the plant operating according to such an overly generous plan will suffer from profit losses. It can be seen from Table 6 that the piecewise-linear approximation method obtains final results in a much shorter solution time than does the Petkov method or the nonlinear approximation method. The final optimal results from all methods are generally very close. In cases 2A and 2B, the two nonlinear methods obtained slightly higher profits. However, this is not always the case. In case 2C, the Petkov method obtained the lowest profit because the MINLP solver stopped searching before finding the optimal results. If all other parts of the planning model (unit models, blending model, etc) are linear, the piecewise-linear approximation is a better choice because it has very high solution speed with good accuracy. If the nonlinearity of the planning model cannot be avoided and the solution time is not very important, one can consider using the Petkov method or the nonlinear approximation method as a comparison. However, because the Petkov method or the nonlinear approximation method introduces further nonlinearities into a planning model, we believe that the piecewiselinear approximation is still a better choice for a nonlinear planning model. 11.3. Case 3: Yearly Planning for a Real-World Refinery Complex. In case 3, the yearly planning model for the D refinery complex is developed. The configuration of the refinery complex is the same as in case 2. The yearly planning model divides the yearly processing amount of crude oil into 12 months optimally with consideration of the operating conditions, market conditions, and other factors over 12 months. A multiperiod site model is developed here to consider all of these factors. The model for the piecewise-linear approximation method involves 204 discrete variables, 5617 single continuous variables, and 4464 single constraints. To reduce the computing time, the relative optimality gap of the OSL solver for all methods is set to 0.05. In case 3, the market demands for 90# and 93# gasoline each month are assumed to conform to normal distributions N(60000,6000) and N(25000,2500), respectively. For the sake of simplicity, the market demands for all other products each month are assumed to be deterministic. The objective of this case study is to

6754


Table 7. Yearly Planning Results of Case 3

Table A.1. Simple Example for Service Levels

method

profit (106 Yuan)

solution time (s)

Petkov nonlinear approximation piecewise-linear approximation

1614.3 1605.9 1618.1

1601.2 333.0 35.8

determine the optimal yearly planning strategy for the D refinery complex under uncertain market conditions. The results are shown in Table 7. The customer satisfaction level is not considered in case 3. Instead, inventory is considered for each product. As for case 2, it can be seen from Table 7 that the piecewise-linear approximation method obtains final results in a much shorter solution time than does the Petkov method or the nonlinear approximation method. For this large-scale multiperiod problem, the advantage of the piecewise-linear approximation method in solution time becomes much more significant. We notice that the Petkov method obtained a worse solution in a much longer solution time. The possible reason for this is that it is more difficult to obtain a low relative optimality gap for a large MINLP model such as the Petkov method. It will take an extremely long time to increase the quality of its solution. The nonlinear approximation method determines the solution in a shorter solution time than does the Petkov method. The possible reason for this difference could be that the nonlinear approximation method has lower nonlinearity. However, the nonlinear approximation method has a worse solution quality. Compared to the performance of the piecewise-linear approximation method, the solution times and the solution statuses of the two nonlinear methods are unstable. Changing any part of the model might incur significantly different solution times. Sometimes, the Petkov method even fails to find a solution when the two approximation methods can. We conclude that, for a large-scale linear industrial planning problem, the piecewise-linear approximation is the best among these three methods.

scenario

product available

product demand

stock-outs

1 2 3 4 5 6 7 8

3000 3000 3000 3000 3000 3000 3000 3000

1200 2000 850 3100 3300 900 2500 1500

0 0 0 100 300 0 0 0

Appendix I. Derivation of the Loss Function for an Exponential Distribution For an exponential distribution with lambda λ, mean µ, and location parameter θ, the probability density function, F(x), and cumulative distribution function, D(x), are written as

F(x) ) λe-λ(x-θ)

D(x) ) 1 - e-λ(x-θ) Then, the loss function, ∫∞P (x - P)F(x) dx, can be expressed as

∫P∞(x - P)F(x) dx ) ∫P∞(x - µ)F(x) dx + ∫P∞(µ - P)F(x) dx )

∫P∞xλe-λ(x-θ) dx - µ∫P∞λe-λ(x-θ) dx +

(µ - P)[1 - D(P)]

)λ

∫P∞xe-λ(x-θ) dx - µ∫P∞λe-λ(x-θ) dx +

(µ - P)e-λ(P - θ)

given that

λ

∫P∞xe-λ(x-θ) dx ) λeλθ∫P∞xe-λx dx 1 ) Pe(λθ-λP) + e(λθ-λP) λ

12. Conclusions In this paper, a novel approximation-based approach for refinery planning under uncertainty is proposed. Compared to methods used in the literature, the proposed approach can obtain good accuracy with a better solution speed. A general formulation for revenue and cost calculation is proposed. The fill rate or type 2 service level objective, which is a greater concern of most managers and often ignored in current chanceconstrained programming, is considered. The tradeoff between the plant service level objectives and the plant profit is captured. Real-world large-scale refinery planning case studies are used to illustrate the effectiveness of the approach proposed in this paper. After discussing the approach proposed in this paper, the engineers in the D refinery complex are very interested in and satisfied with the results.

for x g 0

and

µ

∫P∞λe-λ(x-θ) dx ) µeλθe-λP

Thus, the loss function for the exponential distribution becomes

∫P∞(x - P)F(x) dx ) Pe(λθ-λP) + 1λe(λθ-λP) - µeλθe-λP + (µ - P)e-λ(P-θ) )

eλ(θ - P) λ

Acknowledgment

Appendix II. Simple Example To Illustrate the Difference between Type 1 and Type 2 Service Levels

The authors acknowledge financial support from the Research Grants Council of Hong Kong (Grant HKUST6014/99P & DAG00/01.EG05), the National Science Foundation of China (Grant 79931000), and the Major State Basic Research Development Program (G2000026308).

Table A.1 lists the product demand and stock-outs in eight scenarios for a plant. The amount of product available in the plant is assumed to be 3000 in each scenario. In Table A.1, stock-out occurs in scenarios 4 and 5. The fraction of scenarios in which there is no stock-out


is 6/8 ) 0.75; that is, the probability that all demands are met in a single scenario is 0.75. Thus, the customer satisfaction level is 0.75 according to the definition of the type 1 service level. In Table A.1, the total demands over eight scenarios number 15350. The total stock-outs over eight scenarios number 400. Thus, the proportion of the satisfied demand (fill rate) is (15350 - 400)/15350 ) 0.97. The customer satisfaction level is 0.97 according to the definition of the type 2 service level. We can see that there is a significant difference between these two service levels. A type 1 service level cannot accurately approximate a type 2 service level and should not be used to replace the type 2 service level.27 Literature Cited (1) Li, W.; Hui, C.-W.; Hua, B.; Tong, Z. Scheduling Crude Oil Unloading, Storage, and Processing. Ind. Eng. Chem. Res. 2002, 41, 6723-6734. (2) Clay, R. L.; Grossman, I. E. A disaggregation algorithm for the optimization of stochastic planning models. Comput. Chem. Eng. 1997, 21, 751. (3) Vin, J. P.; Ierapetritou, M. G. Robust Short-Term Scheduling of Multiproduct Batch Plants under Demand Uncertainty. Ind. Eng. Chem. Res. 2001, 40, 4543 (4) Petkov, S. B.; Maranas, C. D. Multiperiod planning and scheduling of multiproduct batch plants under demand uncertainty. Ind. Eng. Chem. Res. 1997, 36, 4864. (5) Dantzig, G. B. Linear programming under uncertainty. Manage. Sci. 1955, 1, 197-206. (6) Beale, E. M. L. On minimizing a convex function subject to linear inequalities. J. R. Stat. Soc. B 1955, 17, 173-184. (7) Horner, P. Planning under uncertainty. Questions & answers with George Dantzig. OR/MS Today 1999, 26, 26-30. (8) Neiro, S. M. S.; Pinto, J. Supply Chain Optimization of Petroleum Refinery Complexes. Presented at the Foundations of Computer Aided Process Operations Conference (FOCAPO 2003), Coral Springs, FL, Jan 12-15, 2003. (9) Kall, P.; Wallace, S. W. Stochastic Programming; Wiley: New York, 1994. (10) Li, P.; Wendt, M.; Wozny, G. Optimal Operations Planning under Uncertainty by Using Probabilistic Programming, Presented at the Foundations of Computer-Aided Process Operations (FOCAPO 2003), Coral Springs, FL, Jan 12-15, 2003. (11) Liu, M. L.; Sahinidis, N. V. Optimiztion in Process Planning under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 4154. (12) Ierapetritou, M. G.; Pistikopoulos, E. N. Batch Plant Design and Operations under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 772. (13) Chen, X. Newton-type methods for stochastic programming. Math. Comput. Modell. 2000, 31, 89-98. (14) Wellons, H. S.; Reklaitis, G. V. The design of multiproduct batch plants under uncertainty with stage expansion. Comput. Chem. Eng. 1989, 13, 115-126.

(15) Rooney, W. C.; Biegler, L. T. Optimal process design with model parameter uncertainty and process variability. AIChE J. 2003, 49, 438. (16) Lee, Y. G.; Malone, M. F. A General Treatment of Uncertainties in Batch Process Planning. Ind. Eng. Chem. Res. 2001, 40, 1507. (17) Uryasev, S. Probabilistic Constrained Optimization: Methodology and Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2000. (18) Henrion, R.; Li, P.; Mo¨ller, A.; Steinbach, M.; Wendt, M.; Wozny, G. Stochastic Optimization for Chemical Processes Under Uncertainty. In Online Optimization of Large Scale Systems; Gro¨tschel, M., Krumke, S. O., Rambau, J., Eds.; Springer-Verlag: New York, 2001; pp 455-476. (19) Gupta, A.; Maranas C. D. A Two-Stage Modeling and Solution Framework for Multisite Midterm Planning under Demand Uncertainty. Ind. Eng. Chem. Res. 2000, 39, 3799. (20) Gupta, A.; Maranas, C. D. Managing demand uncertainty in supply chain planning. Comput. Chem. Eng. 2003, 27, 1219. (21) Applequist, G. E.; Pekny, J. F.; Reklaitis, G. V. Risk and uncertainty in managing chemical manufacturing supply chains. Comput. Chem. Eng. 2000, 24, 2211-2222. (22) Kim, K. J.; Diwekar, U. M. Efficient Combinatorial Optimization under Uncertainty. 1. Algorithm Development. Ind. Eng. Chem. Res. 2002, 41, 1276. (23) Hene, T. S.; Dua, V.; Pistikopoulos, E. N. A hybrid parametric/stochastic programming approach for mixed-integer nonlinear problems under uncertainty. Ind. Eng. Chem. Res. 2002, 41, 67. (24) Balasubramanian, J.; Grossmann, I. E. Scheduling optimization under uncertainty-an alternative approach. Comput. Chem. Eng. 2003, 27, 469. (25) Wendt, M.; Li, P.; Wozny, G. Nonlinear Chance Constrained Process Optimization under Uncertainty. Ind. Eng. Chem. Res. 2002, 41, 3621-3629. (26) Hopp, W. J.; Spearman, M. L. Factory Physics: Foundations of Manufacturing Management; Irwin/McGraw-Hill: New York, 2000. (27) Nahmias, S. Production and Operations Analysis; McGrawHill: New York, 2001. (28) Brooke, A.; Kendrick, D.; Meeraus, A. GAMSsA User’s Guide, release 2.25; The Scientific Press: San Francisco, CA, 1992 (29) Li, W.; Hui, C.-W.; Li, A.-X. Integrating CDU, FCC and Product Blending Models Into Refinery Planning, Comput. Chem. Eng., in review. (30) Reid, E. B.; Allen, H. L. Estimating Pour Points of Petroleum Dist. Blends. Pet. Refin. 1951, 30 (5), 93-95. (31) Jiadong, L. Optimize diesel oil blending using PC-1500. Shiyou Lianzhi 1985, 16 (1), pp 15-18.

Received for review April 1, 2004 Revised manuscript received June 28, 2004 Accepted July 12, 2004 IE049737D