3056
Ind. Eng. Chem. Res. 1999, 38, 3056-3068
Integration and Computational Issues in Stochastic Design and Planning Optimization Problems Fernando P. Bernardo,† Efstratios N. Pistikopoulos,*,† and Pedro M. Saraiva‡ Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College, London SW7 2BY, U.K., and Department of Chemical Engineering, University of Coimbra, 3000 Coimbra, Portugal
In stochastic process design and planning optimization problems, the expected value of the objective function in face of uncertainty is typically evaluated through an n-dimensional integral, where n is the number of uncertain parameters. In this paper, suitable integration techniques are presented and computational issues are discussed in relation to the number of uncertain parameters and the uncertainty model considered. A specialized cubature technique, suitable to integrate normally distributed uncertainties, is introduced, which for n < 10 can reduce significantly the computational effort required when compared to other strategies, such as product Gauss rules or efficient sampling techniques. The computational performance of the different integration techniques and their applicability are discussed through two process-engineering examples. 1. Introduction The inclusion of uncertainty in process design and planning problems has been widely recognized as a key issue to study the effect of variability on decisions related to process performance and quality. Such problems can be posed as the following general mathematical formulation:
max f(d,z,x,θ) d,z,x
s.t. h(d,z,x,θ) ) 0
(1)
g(d,z,x,θ) e 0 d ∈ D, z ∈ Z, x ∈ X, θ ∈ Θ where d and z are the vectors of design and control variables, respectively, x is the vector of state variables, and θ is the vector of uncertain parameters over the domain Θ. The performance metric to be optimized is defined by a scalar function f and the model equality and inequality constraints are represented by the set of equations h and g, respectively. Several approaches have been suggested to formulate and solve this problem, differing in how uncertainty is handled as well as in the design objectives that may include process flexibility, profitability, and/or robustness. In respect to the way uncertainty is handled, three different approaches can be stated: (1) scenario-based approach, where the uncertainty is described by a set of scenarios (periods) with a given probability, resulting in a multiperiod optimization problem;1-5 (2) stochastic approach, where the uncertain parameters are described through a joint probability density function (PDF), resulting in a stochastic optimization problem;6-12 and (3) parametric approach, where the problem is solved parametrically in the space of the uncertain para* To whom all correspondence should be addressed. Tel.: +44-171-5946604. Fax: +44-171-5946606. E-mail:
[email protected]. † Imperial College. ‡ University of Coimbra.
meters.13-16 The present work focuses on the numerical and computational issues of the stochastic approach, as described in the next paragraphs. In the stochastic approach to process design and planning problems under uncertainty, the design objective assumed is typically the optimization of an average process performance value.17 This requires the evaluation of the expected value for a performance metric f, in the face of uncertainty. If the uncertain parameters, θ, follow a joint probability distribution (PDF), j, the expected value of f is given by the following n-dimensional integral, where n is the number of uncertain parameters:
ERn(f) )
∫R
n
f(θ) j(θ) dθ
(2)
The n-dimensional region Rn is a subset of the uncertainty space, Θ, and for a particular design is delimited by the problem constraints. In this paper, we consider a constraint relaxation strategy,18,19 as a result of which the integration region corresponds to the entire uncertainty space:
EΘ(f) )
∫Θ f(θ) j(θ) dθ
(3)
Obtaining accurate estimates for (3) is a considerably intensive computational task, that being one of the main reasons why fully deterministic process design and planning approaches have prevailed for quite a long time. Two groups of approaches have been suggested for obtaining accurate estimates for (3) with less computational effort. On one hand, it is possible to derive easy-to-compute functional approximations for f(θ), especially useful when the original function is complex, and thus demanding considerable effort for each computation of f(θ), for a particular value of θ. The probabilistic collocation method proposed by McRae and co-workers20 in the context of uncertainty analysis falls in such a group. In this method, the model response surface is approximated by polynomial functions of the input uncertain
10.1021/ie9807001 CCC: $18.00 © 1999 American Chemical Society Published on Web 07/03/1999
Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 3057
parameters, and an integration technique (usually a sampling technique) is then applied to this approximation. On the other hand, without introducing any sort of approximations for f(θ), two different integration techniques to estimate (3) can be stated: Gaussian quadrature9 and pseudorandom sampling techniques.8,21 In the first case, a quadrature formula is applied to each dimension of the uncertainty space, resulting in a product Gauss rule with a total number of points that increases exponentially with the number of uncertain parameters. Therefore, this strategy is not computationally efficient to solve large problems with a realistic number of uncertain parameters. In these cases, sampling techniques are more adequate, since the number of points required does not necessarily increase with the integral dimension (the number of uncertain parameters). However, even the most efficient sampling techniques, such as the Hammersley sequence sampling (HSS) introduced by Diwekar and Kalagnanam,21,22 require some hundreds of points to achieve a reasonable accuracy. The two groups of techniques (involving or not f(θ) functional approximations) are not contradictory and may even complement each other, since in the first case one tries to reduce the amount of time needed to compute each f(θ) value, while in the other one tries to come up with accurate estimates for (3) with as few f(θ) values needed as possible. In this article, we develop and compare different methodologies that fall under the second category (where the exact f(θ) function is used, without introducing any approximation to it). In particular, we will focus on how different integration techniques to estimate the multiple integral (3) can be applied to the two-stage stochastic approach to process design and planning problems under uncertainty. The main objective is to identify the most suitable integration technique for a problem with n uncertain parameters described by a given uncertainty model. We introduce a cubature technique,23 suitable to integrate normally distributed uncertainties, which for n < 10 can significantly reduce the computational effort required when compared to other integration strategies, such as product Gauss rules or efficient sampling techniques. The remaining parts of this paper are structured as follows. First, the application of the different integration techniques is illustrated through a simple example. Then, the fundamentals of the cubature integration technique and the set of integration formulas selected to compute (3) are presented. In section 4, the two-stage stochastic formulation is briefly presented, together with the description of the constraint relaxation strategy. Next, we define a heuristic procedure to select the most suitable integration technique to solve a given process design problem. Finally, the proposed procedure is illustrated through two process-engineering examples: a simplified linear planning model of a refinery plant and a nonlinear process synthesis/planning problem. 2. Illustrative Example Let us consider the function f(θ) ) θ12 + θ22, where θ1 and θ2 are uncertain parameters described by a joint normal PDF with mean
µ)
[] 0 0
Figure 1. Grid of points θi for different integration techniques.
and covariance matrix:
Σ)
[ ] 1 0 0 1
The expected value of f in face of the uncertainty is given by the following two-dimensional integral:
EΘ(f) )
1 × ∫-∞+∞∫-∞+∞(θ12 + θ22)2π
[
exp -
1 (θ 2 + θ22) dθ1 dθ2 2 1
]
To approximate the integral above, one needs to evaluate the integrand function for a specific grid of points θi that depends on the integration technique used. Figure 1 shows this grid when four different integration techniques are used: a product Gauss formula, a specialized product Gauss formula, a specialized cubature formula, and an efficient stratified sampling technique. The product Gauss formula (a particular case of the cubature technique) is constructed using a 5-points Gaussian quadrature in each integral dimension, resulting in a 25-points grid which estimates the integral value with a relative error of -8.3% (see Figure 1a). Both specialized formulas are especially constructed for the case where all uncertain parameters are normally distributed and have only 9 points (see Figure 1b,c). The solution obtained by both techniques is exact since the formulas degree (5) is greater than the degree of the function f. The HSS sampling technique leads to a 4.1% error solution when a sample with 50 observations is taken from the multivariate probability space (see Figure 1d). For this example, it is rather clear that the specialized formulas are more accurate and efficient, since they provide the exact solution with less computational effort. However, in the general case, the problem of finding the most suitable technique is not trivial and depends on the integral dimension and the uncertainty model considered. Before addressing this problem in section 5, the fundamentals of the cubature integration technique will be presented in the next section.
3058 Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999
3. Cubature Integration Technique The cubature technique is a numerical integration technique that generalizes the principles of one-dimensional quadrature to multidimensional integration. On the basis of the definition of the quadrature formula given by Engels,24 we define the cubature formula as follows. Definition. A cubature formula is a numerical rule whereby the value of a multidimensional definite integral is approximated by the use of information about the integrand only at a set of discrete points where the integrand is defined. The mathematical notation of the definition above is as follows:23
∫R
Np
n
f(u) w(u) du ≈
Bif(ui) ∑ i)1
(4)
Here, f is a scalar function of the vector of independent variables u, w a given weight function, and Rn the integration region in n-dimensional Euclidean space En. The cubature formula has Np points ui which lie in En and have corresponding coefficients or weights Bi. To approximate integral (3) using a formula like (4), we need to map the region Rn into the uncertainty space Θ. This is done through a general transformation of the form25
θ ) φ(u)
(5)
where u ∈ Rn and θ ∈ Θ. Each point ui in the Rn space has a corresponding point θi in Θ space. The transformed weights are given by
B/i ) Bi|det Jac(ui)|,
i ) 1, ..., Np
(6)
where Jac(ui) is the Jacobian matrix of function φ(u) evaluated at the point ui. Two different classes of cubature formulas can be distinguished: product formulas constructed using combinations of formulas for regions of dimensions less than n and nonproduct formulas constructed by other methods. Referring to the example of the previous section (see Figure 1), the product Gauss formula used, on one hand, is the “product” between two one-dimensional Gaussian quadratures: for each of the 5 quadrature points in the first dimension, there are 5 points across the second dimension, resulting in a total number of points of 5 × 5 ) 25. On the other hand, the specialized cubature is a nonproduct formula, specifically constructed to integrate over a two-dimensional region. In the remaining parts of this paper, we will refer to nonproduct cubatures as only cubatures and to both product and nonproduct cubatures as integration formulas. A cubature formula has a degree d (or degree of exactness d) if it is exact for any linear combination of monomials u1i1u2i2...unin with i1 + i2 + ... + in e d and there is at least one monomial of degree d + 1 for which the formula is not exact.23 This definition applies to a product Gauss formula with the extension of exactness for i1, i2, ..., in e d, since this formula is the product of one-dimensional Gaussian quadratures, each one of degree d. The comparison between product Gauss formulas and cubatures relies on a trade-off between accuracy and computational effort. On one hand, for equal degree and
Table 1. Abcissas and Weights of the One-Dimensional Quadrature N
abcissas (ri)
weights (Ai)
4
(0.339 98 (0.861 14 0.000 00 (0.538 47 (0.906 18
0.652 15 0.347 85 0.588 89 0.478 63 0.236 93
5
weight function, product Gauss formulas are always equally or more accurate than cubatures, but can be numerically very expensive or even impossible to apply in practice to large real problems, since the number of points increases exponentially with the integral dimension. Cubatures, on the other hand, seem to be very promising in terms of efficiency, since for some particular formulas the number of points only increases polynomially with the dimension, but their accuracy must be taken into account. The reduced accuracy of cubature formulas when compared to product Gauss rules may be attributed to three different reasons. First, unlike product Gauss formulas, cubatures may have some points that lie outside the region of integration Rn, which may clearly decrease the overall accuracy. The second cause is related to the coefficients of both formulas: a product Gauss formula has all the coefficients positive, while most of the known cubature formulas, especially those of higher degree, have some points with negative coefficients. Finally, as stated before, a product Gauss formula of degree d is exact for any linear combination of monomials u1i1u2i2...unin with i1, i2, ..., in e d, while a cubature formula of degree d is only exact when i1 + i2 + ... + in e d. This situation gets worse as the dimension grows, and we call it the dimension effect. The first two factors are critical and may reduce significantly the cubature accuracy, while the dimension effect is only significant for some particular highly nonlinear integrand functions. The dimension effect will be discussed in the context of the process examples of section 6. In the subsequent subsections, we select a set of integration formulas suitable to estimate the multiple integral (3). To overcome some of the potential limitations outlined above, the selected formulas must satisfy the following three conditions: (A) all the points must lie inside the integration region; (B) all the coefficients must be positive, and (C) the formula degree must be sufficiently high to integrate highly nonlinear integrands, such as a performance function f times a joint normal PDF. We consider two types of formulas: (1) formulas to integrate over a n-dimensional cube with unit weight function and (2) formulas to integrate the entire n-dimensional space with a weight function of the form w(u) ) e-u2. The first ones can be applied to any uncertainty model (once the joint PDF j is defined in the considered n-cube region), while the second formulas can only be applied for the case where the uncertain parameters are normally distributed. The simple example introduced in section 2 will be used to illustrate the application of the different integration formulas. 3.1. General Integration Formulas. Product Gauss Formula. One of the most common integration formulas is the product Gauss formula for the ndimensional cube, Cn ) {uj: -1 e uj e 1, j ) 1, ..., n}, constructed using a Gaussian quadrature over the interval [-1, 1] in each integral dimension. The Gaussian quadrature formula26 with N abcissas ri and weights Ai given in Table 1 has a degree of exactness, d ) 2N -
Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 3059
1. The resulting n-dimensional product formula has Np ) Nn points and the same degree d. The grid of points ui and corresponding weights Bi are of the following form:
(ri1, ri2, ..., rin) Ai1Ai2...Ain with each of the subscripts i1, i2, ..., in ranging independently over the integers 1, 2, ..., N. Note that the points constructed in this way always lie inside Cn (condition A satisfied) and always have positive weights (condition B satisfied). To satisfy condition C, we only consider formulas with N g 4. To construct an equivalent formula over the uncertainty space Θ ) {θj: θjL e θj e θjU, j ) 1 ..., n}, we perform the following transformation:
θj(uj) )
L L θU θU j + θj j - θj + uj , 2 2
j ) 1, ..., n (7)
The transformed product Gauss formula, labeled as PGd, is then n
EΘ(f) ≈
∏ j)1
L θU j - θj
2
abcissas (ri)
weights (Ai)
2 3
(1 0 (x3
1/2 2/3 1/6
respect to the performance function f(θ) and not with respect to the entire integrand f(θ) j(θ). Specialized Product Gauss Formula. Specialized one-dimensional quadratures can be constructed30 to integrate over the interval (-∞, +∞) with the weight function
wg(r) )
Bif[θ(ui)]j[θ(ui)] ∑ i)1
∫-∞+∞...∫-∞+∞ f(θ) jN(θ) dθ
(8)
(9)
1 1 exp - (θ - µ)TΣ-1(θ - µ) 2 (2π)n/2(det Σ)1/2 (10)
[
]
We shall denote this joint distribution by N(µ,Σ) and the corresponding uncertainty space by Θ ) {θ: θ ∼ N(µ,Σ)}. For this particular case, specialized formulas can be very advantageous.29 We present a specialized product Gauss formula and specialized cubatures. These formulas are constructed to integrate the entire n-dimensional space and thus condition A is trivial. To satisfy condition C, we only consider formulas with d g 3. It must be noted that the degree of these formulas is defined with
( )
1 r2 exp 2 x2π
(11)
The formula with N abcissas ri and corresponding weights Ai (see Table 2) has a degree of exactness d ) 2N - 1. On the basis of this specialized quadrature, a specialized product Gauss formula can be constructed with the weight function
wpg(u) )
where jN(θ) is the joint normal distribution28 of the random vector θ, with vector of means µ and covariance matrix Σ:
jN(θ) )
N
Np
n-Cube Cubatures. As stated before, product Gauss formulas such as the formula PGd can be very accurate but computationally expensive when applied to large problems. In these cases, cubatures can provide a more efficient way for obtaining multiple integrals estimates. We have tested the performance of several n-cube cubatures presented by Stroud23 and Genz and Malik.27 Unfortunately, none of them satisfies simultaneously both conditions (B) and (C): the formulas with all the coefficients positive do not have a sufficiently high degree (only formulas with degree d e 5 have all the coefficients positive) and the higher degree formulas always have some points with negative coefficients. This makes the n-cube cubatures not very suitable to approximate integral (3). 3.2. Specialized Integration Formulas. When all the uncertain parameters are normally distributed, the integral (3) over the entire probability space takes the following form:
EΘ(f) )
Table 2. Abcissas and Weights of the One-Dimensional Specialized Quadrature
1 1 exp - uTu n/2 2 (2π)
(
)
(12)
The formula has a total number of points given by Np ) Nn and the same degree d. The grid of points ui and corresponding weights Bi are of the following form:
(ri1, ri2, ..., rin) Ai1Ai2...Ain with each of the subscripts i1, i2, ..., in ranging independently over the integers 1, 2, ..., N. Once more, the points constructed in this way satisfy condition B. To estimate integral (9) using this product formula, we perform the transformation θ(u) ) µ + Σ1/2u, which substituted in (9) leads to
∫
EΘ(f) )
+∞
∫-∞
... -∞
+∞
Np
f[θ(u)] wpg(u) du ≈
Bif[θ(ui)] ∑ i)1
(13)
We label this specialized product Gauss formula as SPGd. Specialized Cubature Formulas. Specialized cubature formulas can be obtained based on the rules reported by Stroud23 and labeled as Enr2. These rules are constructed to integrate the entire n-dimensional space with the weight function
wc(u) ) exp(-uTu)
(14)
To integrate over the uncertainty space Θ ) {θ: θ ∼ N(µ,Σ)}, we use the transformation θ(u) ) µ + Ix2Σ1/2u, where Ix2 is a diagonal matrix with all the diagonal elements equal to x2. The substitution of this transformation in (9) leads to
EΘ(f) )
1 π
+∞ +∞ ...∫-∞ f[θ(u)] wc(u) du ≈ ∫ -∞ n/2
1
Np
∑ Bif[θ(ui)]
πn/2 i)1
(15)
We present specialized cubatures of degrees 3, 5, and 7.23
3060 Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999
The third degree formula is labeled as SC3 and it has only Np ) 2n points. The grid of points ui and corresponding weights Bi are of the following form:
(r, 0, ..., 0)FS B0 where FS denotes a fully symmetric set of points generated by permutation of the elements and its signs. The parameter r and the coefficient B0 are functions of the dimension n and the n-dimensional volume, V ) πn/2:
n r ) 2 2
B0 )
1 V 2n
This formula only satisfies condition B for n e 4. Finally, the seventh degree formula, SC7, defined for n g 3, has a number of points given by Np ) 2n+1 + 4n2. The complete grid of points ui and corresponding weights Bi are of the form
(r.r1, 0, 0, ..., 0)FS B.A1 (r.r2, 0, 0, ..., 0)FS B.A2 (s.r1, s.r1, s.r1, ..., s.r1)FS C.A1 (s.r2, s.r2, s.r2, ..., s.r2)FS C.A2
(16)
(t.r1, t.r1, 0, ..., 0)FS D.A1 (t.r2, t.r2, 0, ..., 0)FS D.A2
This formula has positive weights for all n (condition B satisfied). Two fifth-degree formulas are presented. The first one, labeled as SC5,1, is only defined for n g 3 and has a total number of points given by Np ) 2n + 2n. The grid of points ui and corresponding weights Bi are
The parameter values are given by the following expressions:
r)1 s2 )
1 n
t2 )
1 2
(r, 0, ..., 0)FS B0 (s, s, ..., s)FS B1 The parameters r and s and coefficients B0 and B1 are given by
r2 ) s2 ) B0 ) B1 )
n+2 4
n+2 2(n - 2)
(17)
4 V (n + 2)2 (n - 2)2
C)
2-nn3 V n(n + 2)(n + 4)
D)
4 V n(n + 2)(n + 4)
V 2n(n + 2)2
(0, 0, 0, ..., 0) B0 (r, 0, 0, ..., 0)FS B1 (s, s, 0, ..., 0)FS B2 with parameters r and s and coefficients Bi given by
r2 )
n+2 2
s2 )
n+2 4
B0 )
2 V n+2
4-n V 2(n + 2)2
B2 )
1 V (n + 2)2
n + 2-x2(n + 2) 2 (19)
8-n B) V n(n + 2)(n + 4)
This formula has positive weights for all dimensions (condition B satisfied). The second formula is labeled as SC5,2 and it has a number of points given by Np ) 2n2 + 1. The complete grid of points ui and corresponding weights Bi are of the form
B1 )
r12, r22 )
(18)
A1, A2 )
n + 2 ( x2(n + 2) n Γ 2 4(n + 2)
()
∫
∞
The Γ function is given by Γ(z) ) 0 tz-1e-t dt or simply Γ(z) ) (z - 1)!, when z is a positive integer. This formula only satisfies condition B for n e 8. 3.3. Illustrative Example. Here, three different integration formulas are applied to the motivating example of section 2: the product Gauss formula PG9, the specialized product Gauss formula SPG5, and the specialized cubature SC5,2. In the first case, the uncertainty space is truncated by considering a band of 3.09σ around each mean (99.8% probability of occurrence). In this way, θL1 ) θL2 ) -3.09 U and θU 1 ) θ2 ) 3.09. Next, the grid of 25 points ui and corresponding weights Bi are constructed in the following way: ui ) (ri1, ri2) and Bi ) Ai1Ai2, with the subscripts i1 and i2 ranging independently over the integers 1, 2, ..., 5 (the abcissas ri and weights Ai are given in Table 1, for N ) 5). Then, the grid of points θi shown in Figure 1a is obtained using the transformation (7). Finally, the integral approximation is given by eq 8, with EΘ(f) ≈ 1.834. The specialized product Gauss formula SPG5 has 9 points ui and corresponding weights Bi generated in the following way: ui ) (ri1, ri2) and Bi ) Ai1Ai2, with the subscripts i1 and i2 ranging independently over the
Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 3061
integers 1, 2, 3 (the abcissas ri and weights Ai are given in Table 2, for N ) 3). The transformation θ(u) ) µ + Σl/2u gives us the grid of points θi represented in Figure 1b. Using the formula in (13), we obtain the exact solution EΘ(f) ) 2. To apply the specialized cubature SC5,2, we begin by computing the expressions in (18). For n ) 2 and V ) π, we have r ) x2, s ) 1, B0 ) π/2, B1 ) π/16, and B2 ) π/16. Next, we generate the grid of 9 points ui: (0,0), (r,0), (-r,0), (0,r), (0,-r), (s,s), (-s,s), (s,-s), and (-s,s). Then, we perform the transformation θ(u) ) µ + Ix2Σ1/2u, obtaining the points θi represented in Figure 1c. Finally, the formula in (15) gives us the exact solution EΘ(f) ) 2.
4. Two-Stage Stochastic Optimization Formulation The two-stage stochastic optimization formulation is widely recognized as an effective approach to processengineering problems under uncertainty, differentiating between design and control variables. The design variables are considered “here and now” decisions, while the control variables are “wait and see” decisions. The assumption is that during operation there is enough information about the uncertain parameters to optimally adjust the control variables. This suggests a twostage strategy,9 as follows:
the design variables, this being the major difficulty for solving problem 20.31 In this work, partial feasibility is considered by relaxing the constraints that, for some realizations of the uncertain parameters, are responsible for infeasibility in the inner optimization problem and penalizing partially feasible scenarios.18,19,31 From an engineering point of view, this is equivalent to allowing process operation for every possible realization of the uncertain parameters penalizing process performance when a specification is not met. The exact format for the penalty function depends on the design objective considered. In the process examples of section 6, a linear penalty term is considered; when process quality is a design objective, the penalty function may rather refer to a quality cost, usually a quadratic term.12,19 With this relaxation strategy, the integration region Rn is the entire uncertainty space Θ and the expected value of f′ is now given by the following n-dimensional integral (equivalent to integral 3):
EΘ{f′(d,θ)} )
∫Θ f′(d,θ) j(θ) dθ
(23)
The integral above can be approximated either by one of the integration formulas presented in the last section or by use of a sampling technique. In both cases, the formulation 20 can be recast as a single-level optimization problem: Np
Design stage:
max max ERn(d){f′(d,θ)}
d,zi,xi
d
wi f(d,zi,xi,θi) ∑ i)1
s.t. h(d,zi,xi,θi) ) 0
d ∈ D, θ ∈ Θ
g(d,zi,xi,θi) e 0 Operating stage:
(20)
f′(d,θ) ) max f(d,z,x,θ)
(24)
d ∈ D, zi, ∈ Z, xi ∈ X i ) 1, ..., Np
z,x
s.t. h(d,z,x,θ) ) 0 g(d,z,x,θ) e 0 z ∈ Z, x ∈ X The design variables are selected in the first (design) stage and remain fixed in the second (operating) stage. The objective of the operating stage is to determine an optimal vector of control variables z for each possible realization of the uncertainties θ lying within the corresponding feasible region Rn(d):
Rn(d) ) {θ ∈ Θ| ∃ (z,x): h(d,z,x,θ) ) 0 ∧ g(d,z,x,θ) e 0} (21) If the uncertain parameters are described by a joint PDF j, the expected value of f′ is given by the following n-dimensional integral (n is the number of uncertain parameters) over the inner optimization problem:
ERn(d){f′(d,θ)} )
∫R
n(d)
f′(d,θ) j(θ) dθ
(22)
Since the model constraints are considered as hard constraints, the integration region Rn is a function of
Here, θi is the point i in the integration formula with corresponding weight wi (which includes the joint probability j(θi)). If a sampling technique is used, the points θi designate the No observations in the sample considered, all of them with weight wi ) 1/No. The objective function f includes the penalty terms corresponding to partial feasible scenarios. This formulation has the advantage of keeping the same structure as the corresponding deterministic problem; therefore, if the deterministic problem is convex, the resulting stochastic problem of the form (24) will also be convex. The number of points needed in the integration formula, Np, is the critical factor that determines the computational effort required to solve problem 24. In the next section, we provide some advice on how to select an integration technique with the minimum number of points.
5. Choice of a Suitable Integration Technique Here, we try to address the following problem: given a process design problem with n uncertain parameters described by a given uncertainty model, what is the most suitable integration technique or combination of
3062 Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 Table 3. Integration Formulas PG9 product Gauss formula formula
EΘ(f) ≈ integration region degree no. of points points ui and weights Bi transformation θ ) φ(u)
n
L θU j - θj
j)1
2
∏
Np
∑ B f [θ(u )]j[θ(u )] i
i
i
i)1
Θ ) {θj: e θj e j ) 1, ..., n} d ) 9 with respect to f(θ)j(θ) Np ) 5n (ri1, ri2, ..., rin) Ai1Ai2...Ain (see Table 1) with i1, i2, ..., in ranging independently over 1, 2, ..., 5 L U L θj(uj) ) (θU j + θj )/2 + uj(θj - θj )/2, j ) 1, ..., n θLj
θU j ,
SPG5 specialized product Gauss formula formula integration region degree no. of points points ui and weights Bi transformation θ ) φ(u)
Np EΘ(f) ≈ ∑i)1 Bif[θ(ui)] Θ ) {θ: θ ∼ N(µ,Σ)} µ, vector of means; Σ, covariance matrix d ) 5 with respect to f(θ) Np ) 3n (ri1, ri2, ..., rin) Ai1Ai2...Ain (see Table 2) with i1, i2, ..., in ranging independently over 1, 2, 3 θ(u) ) µ + Σ1/2u
SC5,1 specialized cubature formula formula
EΘ(f) ≈ integration region degree no. of points points ui and weights Bi
1 πn/2
Np
∑ B f [θ(u )], i
i
ng3
i)1
Θ ) {θ: θ ∼ N(µ,Σ)} µ, vector of means; Σ, covariance matrix d ) 5 with respect to f(θ) Np ) 2n + 2n (r, 0, ..., 0)FS B0 (s, s, ..., s)FS B1
n+2 4 n+2 2 s ) 2(n - 2) 4 B0 ) V (n + 2)2 (n - 2)2 B1 ) n V 2 (n + 2)2 r2 )
transformation θ ) φ(u)
techniques to solve problem 24? The set of options considered in this study includes the HSS sampling technique and the integration formulas PG9, SPG5, and SC5,1. We have chosen these particular formulas, because their degree is in principle adequate for most process systems applications. Table 3 shows these integration formulas together with their main properties. The criterion employed to select the most suitable integration technique is the computational effort required to achieve a satisfactory accuracy. In the case of the integration formulas, this criterion is directly related to the number of points in the formula, which is a function of the problem dimension n. When choosing between the formulas SPG5 and SC5,1, we assume that they are equally accurate, i.e., the dimension effect is neglected. In the case of the HSS technique, the computational effort is related to the number of observations, No, required to converge the solution to within a band error. Unfortunately, there are no analytical approaches to calculate No for stratified techniques, as is the case of the HSS technique, and thus No must be evaluated solving the problem several times with an increasing number of observations, until the solution is within a desired band error. Literature results12,21,22 indicate that
V ) πn/2 θ(u) ) µ + Ix2Σ1/2u Ix2, diagonal matrix with all the diagonal elements equal to x2
the HSS technique is much more efficient than other sampling techniques (Monte Carlo or Latin Hypercube sampling) and that good approximations (relative error less than 1%) can be obtained with values of No on the order of 200-600. On the basis of these results and on the example problems reported in section 6, the following heuristic rules are proposed to choose between the most efficient integration formula and the HSS technique: (1) if the integration formula has a number of points below 200, then the integration formula is probably the best choice; (2) if the integration formula has a number of points between 200 and 600, then the use of the HSS technique should be discussed; (3) if the number of points is greater than or equal to 600, then the HSS technique is probably the best choice. Rule 2 corresponds to cases in which a consistent choice becomes difficult. For instance, if the integration formula is the specialized cubature SC5,1, then the problem dimension is 8 or 9. For these dimensions, if the problem is particularly nonlinear, the dimension effect can be somewhat significant and then the sampling technique may become preferable. However, there are no guarantees for No to be smaller than 600, and even if it is, there is always a considerable computational effort associated with its evaluation. In these cases, the cubature technique may become more ef-
Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 3063 Table 4. Number of Points for the Product Cubature Formula n/n1 1 2 3 4 5 6 7 8 9 10
0
1
2
3
4
5
6
7
8
9
10
5 25 125 HSS HSS HSS HSS HSS HSS HSS
3 15 75 375 HSS HSS HSS HSS HSS HSS
9 45 225 HSS HSS HSS HSS HSS HSS
14 70 350 HSS HSs HSS HSS HSS
24 120 HSS HSS HSS HSS HSS
42 210 HSS HSS HSS HSS
76 380 HSS HSS HSS
142 HSS HSS HSS
272 HSS HSS
530 HSS
HSS
ficient. We will discuss these issues in the context of the process examples presented in section 6. Let us now consider the particular case where all uncertain parameters are described by independent or correlated normal distributions. The product Gauss rule PG9 is eliminated a priori due to its low efficiency when compared to the specialized product Gauss formula SPG5 (5n points compared to 3n). For n ) 1 and n ) 2 the specialized cubature SC5,1 is not defined and thus the formula SPG5 is chosen. As the number of uncertain parameters increases, the specialized cubature becomes more efficient. For large values of n, the use of the HSS technique should be discussed based on the heuristic rules stated above. In the general case, the uncertain parameters are described by independent or correlated PDFs of any kind (normal, beta, uniform, log-uniform, etc.), and thus specialized formulas cannot be applied. However, we can construct a product cubature formula using an efficient specialized rule to integrate over the subset of uncertainties that are normally distributed. In this way, the number of points is significantly reduced when compared to the one obtained when Gaussian quadrature is applied to all dimensions. Let θ1 be the n1 uncertain parameters that are normally distributed. A product cubature formula can be constructed in the following way. First, a specialized formula is applied to the n1 dimensions. If n1 ) 1 or n1 ) 2, we use the specialized product Gauss SPG5; if n1 g 3 the best formula is the specialized cubature SC5,1. The remaining (n - n1) dimensions are integrated using a product Gauss formula such as PG9. The number of points for the resulting formula is the product of the number of points involved in the two initial formulas:
n1 ) 0 w Np ) 5n n1 ) 1 w Np ) 3 × 5n-1 n1 ) 2 w Np ) 9 × 5n-2 n1 g 3 w Np ) (2n1 + 2n1) × 5n-n1 The choice between the product formula and the HSS sampling technique is based on the heuristic rules stated above. In Table 4 we present the values of Np for 1 e n e 10 and 0 e n1 e n. When Np g 600, we write “HSS” instead of the number of points, since in these cases the HSS sampling technique is probably the best option. Now, we revisit the problem introduced in the beginning of this section. Given a process design problem with n uncertain parameters described by a given uncertainty model, the choice of the most suitable integration
Figure 2. Refinery inputs and outputs.
technique to solve problem 24 is based on the following heuristic procedure: 1. If the uncertain parameters are normally distributed: (a) If n ) 2, then use the formula SPG9. (b) If 3 e n e 7, then use the formula SC5,1. (c) If n ) 8 or n ) 9, then use the formula SC5,1 or the HSS sampling technique. (d) If n g 10, then use the HSS sampling technique. 2. If only n1 of the uncertain parameters are normally distributed, then construct a product cubature formula with Np points (Np given by Table 4): (a) If Np e 200, then use the product cubature formula. (b) If 200 < Np < 600, then use the product cubature formula or the HSS sampling technique. (c) If Np g 600, then use the HSS sampling technique. 6. Case Studies This section aims to illustrate the usefulness of the proposed procedure to select the most suitable integration technique. Two process-engineering examples are considered: a linear planning problem and a nonlinear synthesis/planning problem. 6.1. Linear Planning Problem. The problem considered here is a variation of the refinery planning model presented in Edgar and Himmelblau,32 schematically shown in Figure 2. With a division in the time horizon of operation into two stages, the objective is to maximize profit over both stages. The profit function is defined as the difference between the revenue due to product sales and overall cost (cost of raw materials plus operating costs). We also include a penalty term for unfilled product demands. The mathematical model, together with its nomenclature, are presented in Table 5. The yield matrix is reported in Table 6 and the parameter values in Table 7. During the first stage, parameter values are assumed to be known, while over the second stage they can only be forecasted and given by a PDF. If we assume that during the second stage there is enough information about the uncertain parameters to optimally adjust the planning variables, then the planning problem is formulated as a two-stage stochastic problem of the form (20), where the first decision variables (design variables) that have to be decided prior to resolution of uncertainty are d ) {Pi,1, Ri,1, Ii,1in, Ii,lf, Sj,1, Ij,1in, Ij,1f, Ii,2in, Ii,2in, Ij,2in};
3064 Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 Table 5. Planning Model material balances capacity limits inventory constraints
purchase limits sales limits objective function (k$) index sets i j s
Table 8. Uncertainty Model 1
in f Pi,s - Ri,s + Ii,s - Ii,s )0 in f ∑iyi,jRi,s + Ij,s - Ij,s - Sj,s ) U 0 e ∑iyi,jRi,s e Cj,s (j * 4) f f U 0 e Ii,s e I , 0 e Ij,s e IU in in Ii,1 ) 0, Ij,1 )0 f in f in Ii,1 ) Ii,2 , Ij,1 ) Ij,2 0 e Pi,s e Ai,s
0
0 e Sj,s e Dj,s profit ) ∑s∑jPrj,sSj,s - ∑s∑iCi,sPi,s ∑s∑iCoi,sRi,s - k∑s∑jPrj,s(Dj,s - Sj,s) crude oil, i ) 1, 2 product, j ) 1, ..., 4 stage, i ) 1, 2
variables (kbbl) Pi,s Ri,s in f Ii,s and Ii,s
parameters yi,j Cj,sU
yield of crude oil i in product j maximum production capacity of product j at stage s inventory limit (IU ) 2kbbl) availability of crude oil i at stage s demand of product j at stage s sale price of product j at stage s purchase price of crude oil i at stage s operating cost relative to crude oil i at stage s penalty coefficient (k ) 1)
IU Ai,s Dj,s Prj,s Ci,s Coi,s k
Table 6. Yield Matrix gasoline kerosene fuel oil residual
crude oil 1
crude oil 2
0.8 0.05 0.1 0.05
0.44 0.1 0.36 0.1
Table 7. Parameter Valuesa A (kbbl) C ($/bbl) Co ($/bbl)
CU
(kbbl) D (kbbl) Pr ($/bbl) a
PDF
D12 D22 D32 D42 A12 A22 Pr12 Pr22 Pr32 Pr42
N(20,1.118) N(2.5,0.5) N(5,1) N(2,0.5) N(20,1.5) N(15,1) N(36,2) N(24,2) N(21,1.5) N(10,1)
Table 9. Results for Uncertainty Model 1
amount of crude oil i purchased at stage s amount of crude oil i utilized at stage s initial and final inventory of crude oil i at stage s amount of product j sold at stage s initial and final inventory of product j at stage s
Sj,s in f Ij,s and Ij,s
parameter
crude oil 1
crude oil 2
20 24 0.5
15 15 1
gasoline
kerosene
fuel oil
residual
24 20 36
15 2.5 24
10 5 21
2 10
Equal values at both stages.
the second decision variables (control variables) that can be optimally adjusted to uncertainty are z ) {Pi,2, Ri,2, Ii,2f, Sj,2, Ij,2f}. If we choose to always meet the market demand, Sj,s ) Dj,s, then the process will possibly be infeasible for some realizations of the uncertain parameters and the integration region in (22) will be a function of the design variables. Instead, we will consider unfilled demands by relaxing the constraint Sj,s ) Dj,s to Sj,s e Dj,s and penalizing the objective function with the loss of opportunity associated with the unfilled scenarios, Prj,s(Dj,s - Sj,s) (penalty coefficient, k ) 1).18 The planning problem is then formulated as a single-level optimization problem of the form (24) which exhibits the same linear structure of the corresponding deterministic problem. First, we consider the case in which all the uncertain parameters are described by normal PDFs N(µj,σj), j )
n
technique
Np
CPU (s)
E(profit) (k$)
2 3 4 5 6 8 8 10
SPG5 SC5,1 SC5,1 SC5,1 SC5,1 SC5,1 HSS HSS
9 14 24 42 76 272 300 300
0.9 1.8 4.2 10.0 35.3 615.3 609.6 250.6
485.2 483.0 478.1 476.2 475.1 474.4 478.0 478.4
1, ..., n. The joint normal PDF is denoted by N(µ,Σ), where µ is the vector of means and Σ the covariance matrix. Assuming that all the uncertain parameters are independent, the matrix Σ becomes equal to the variance matrix, V (diagonal matrix with variances as diagonal elements). The problem is solved for an increasing number of uncertain parameters, from n ) 2 to n ) 10, according to the sequence of uncertain parameters listed in Table 8, and using GAMS/CONOPT2 as optimization routines.33 The samples corresponding to the HSS technique are simulated using a Fortran subroutine (SPARCLHS2) that implements this and other sampling techniques, allowing the user to consider several kinds of input PDFs (independent or correlated, continuous or discrete). Table 9 shows the results obtained when the integration technique is chosen following the procedure proposed in section 5. The optimal design is the same for all problem dimensions:
d* ) {P/1,2, P/2,2, R/1,2, R/2,2, S/1,2, S/2,2, S/3,2, S/4,2} ) {20, 15, 19.836, 13.934, 20, 2.385, 5, 2} To clarify the application of the cubature technique to the planning problem of the form (24), let us consider in more detail the case where n ) 5. For this problem dimension, the specialized formula SC5,1 has 2n + 2n ) 42 points ui with corresponding weights Bi. The θi points in the normally distributed uncertainty space are obtained using the transformation θ(u) ) µ + Ix2Σ1/2u. The planning problem is then formulated in the form of problem 24 and the cubature formula 15 is then used to construct the objective function. The total number of decision variables is given by dim(d) + 42 dim(z). For n ) 8, the most suitable integration technique cannot be convincingly obtained by the proposed procedure, since two different techniques are suggested: the specialized cubature SC5,1 and the HSS sampling technique. Figure 3 shows the solutions obtained by the HSS technique when different sample sizes are used and also by the specialized cubature formulas SC3 with 16 points, SC5,1 with 272 points, and SC7 with 768 points (actually, only 736 points since 32 points have weight ) 0). The band error shown is constructed considering a relative error of 0.5% around the HSS solution with
Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 3065 Table 11. Loss of Accuracy with Formula SC5,2
a
n
Np
errora (%)
3 4 5 6
19 33 51 73
-0.1 0.0 27.6 50.9
Relative to SPG5 solution.
Table 12. Uncertainty Model 2 parameter
PDF
D12 D22 D32 D42 A12
N(20,1.118) N(2.5,0.5) N(5,1) LU(1,3) U(17,23)
Figure 3. HSS vs cubatures (uncertainty model 1, n ) 8).
Figure 5. Process superstructure.
Figure 4. HSS vs cubatures (uncertainty model 1, n ) 10). Table 10. Dimension Effect on Formula SC5,1
a
n
errora (%)
Np(SPG5)/Np(SC5,1)
3 4 5 6
0.6 0.0 0.1 0.1
1.9 3.8 5.8 9.6
Relative to SPG5 solution.
1200 observations. As we can see, the choice between the formula SC5,1 and the HSS technique is in fact dubious: both techniques require approximately the same number of points to estimate a solution within the band error (300 points, in the case of the HSS technique, and 272 points, in the case of the cubature formula). For n ) 10 (see Figure 4), the choice is clearer: the HSS technique is in fact more efficient since it only requires 300 observations to converge the solution to within the band error, compared to 1044 points of the cubature formula SC5,1. The proposed procedure does not take into account the dimension effect when the formula SC5,1 is preferred to the formula SPG5. To predict the dimension effect, we would need information about the nonlinearity of the integrand function. In our case, since the integrand is an optimization problem (see eq 23), we would have to solve it parametrically as a function of the uncertain parameters (see Acevedo and Pistikopoulos,34 for a hybrid parametric/stochastic approach). However, for this particular example, we can show that the dimension effect is negligible, since the results obtained by the specialized cubature are essentially the same as the ones obtained when the product Gauss formula is used (see Table 10), and with significantly less effort. To illustrate the effect of negative weights on the accuracy of integration formulas, we have also compared the performance of the formulas SC5,1 and SC5,2. Formula SC5,2 is very promising, since its number of points only increases quadratically with dimension (Np ) 2n2 + 1). Unfortunately, as shown in Table 11, the formula
is not suitable for n g 5, resulting in high inaccuracies due to negative weights for some of its points. Let us now consider the uncertainty model in Table 12 with different probability distributions. The demands for the three first products in the second stage are described by independent normal PDFs (n1 ) 3), while the demand for product 4 follows a log-uniform distribution (LU) between 1 and 3, and the availability of the first crude oil is uniformly distributed between 17 and 23. We consider a log-uniform distribution with base 10 logarithms, given by
j4(θ4) )
1 1 ) L (ln 3)θ4 (ln 10)θ4(log10 θU log θ ) 4 10 4
In this case, none of the specialized formulas can be used. Instead, we construct a product cubature formula, using the specialized cubature SC5,1 to integrate the first three dimensions of the uncertainty space, while the other two dimensions are integrated using the product Gauss rule PG9. The total number of points is Np ) (2n1 + 2n1) × 5n-n1 ) 14 × 25 ) 350 (note that if we were using Gaussian quadrature in all dimensions, the number of points would be Np ) 55 ) 3125). Using this product formula to solve problem 24, a CPU time of 356.7 s is required to obtain the optimal solution E(Profit) ) 472.8 k$, with a corresponding design identical to the one obtained with the uncertainty model 1. We have also solved the problem using the HSS sampling technique. A sample of 200 observations is enough to converge the solution to within a 0.5% band error, revealing that, in this case, the HSS technique is more efficient. The optimal expected profit obtained is 474.6 k$, with a CPU time of 246.8 s. 6.2. Nonlinear Process Synthesis/Planning Problem. Now we will consider the process superstructure shown in Figure 5 for the production of C, in which the feedstock to process 1 can be produced by process 2, and either process 3 and/or 4.35 A restricted amount of the intermediate product B can also be obtained in the market and used to produce C. The mathematical model for the synthesis/planning problem and corresponding nomenclature are presented in Table 13, while Table
3066 Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 Table 13. Synthesis/Planning Model material balances process conversions
demand limitation availability limitations production limitations objective function
variables y3, y4
Table 16. Results for Uncertainty Model 1
A 2 + A3 + A 4 ) A B2 + B3 + B4 ) B C ) 0.9(Bf + B) B2 e MP2 ln(1 + A2/k2) B3 e MP3 ln(1 + A3/k3) B4 e MP4 ln(1 + A4/k4) C e DC A e AA Bf e A B A2 g 5 A3 e 20y3 A4 e 20y4 profit ) PrCC - (PrAA + PrBBf [FC1 + 15(Bf + B)] - (FC2 + 5A2) (FC3y3 + 15A3) - (FC4y4 + 5A4) kPrC(DC - C) binary structural variables representing the existence (or not) of units 3 and 4 utilization of A at processes 2, 3, and 4, respectively total amount of A production of B at processes 2, 3, and 4, respectively total production of B the amount of material B purchased production of final product
A2, A3, A4 A B2, B3, B4 B Bf C parameters DC, AA, AB MP2, MP3, MP4 k2, k3, k4 PrC, PrA, PrB FC1, FC2, FC3 k
demand of C and availability of A and B production constants kinetic constants prices of C, A, and B fixed cost constants penalty coefficient (k ) 1)
Table 14. Parameter Values process
MP
1 2 3 4 DC AA AB
18 20 15 20 15 15
k
FC
20 21 26 PrC PrA PrB
100 80 130 150 550 250 300
Table 15. Uncertainty Model 1
a
parameter
PDF
DC AA AB k2 k3 k4 PrA PrB
N(20,2) N(15,1) N(15,1) N(20,2) N(21,4) N(26,4) N(250,6)a N(300,5)a
Correlation coefficient ) 0.7.
14 provides the parameter values. The profit function includes a revenue term, an overall costs (costs of raw material + fixed costs + operating costs) and a penalty corresponding to the unfilled demands. Note that the nonlinear equations for process conversions are relaxed (B2 e MP2 ln(1 + A2/k2) instead of B2 ) MP2 ln(1 + A2/ k2)), to obtain an equivalent convex reformulation. With some of the model parameters uncertain, the problem can be formulated as a two-stage stochastic program in which the first decision variables are the binary decisions, d ) {y3, y4}, and the flow rates stand as control variables, z ) {A2, A3, A4, B2, B3, B4, Bf, C}. Once more, the relaxed demand constraint C e DC guarantees feasibility over the entire uncertainty space and the problem can be formulated as a single-level convex MINLP of the form (24).
n
technique
Np
CPU (s)
E(profit)
2 3 4 5 6 8 8
SPG5 SC5,1 SC5,1 SC5,1 SC5,1 SC5,1 HSS
9 14 24 42 76 272 300
4.1 5.8 9.8 20.6 32.2 673.3 755.0
3576 3560 3570 3635 3636 3630 3610
Table 17. Uncertainty Model 2
a
parameter
PDF
DC PrC AA AB k2 k3 k4
U(18,22)a U(536,564)a N(15,1) N(15,1) N(20,2) N(21,4) N(26,4)
Correlation coefficient ) 0.7.
Let us first consider an uncertainty model defined by a joint normal distribution N(µ,Σ). We solve the problem for an increasing number of uncertain parameters, according to the sequence given in Table 15, and using GAMS/(MINOS5,OSL) as optimization routines.33 For dimension n ) 8, the prices for A and C are considered to have a positive correlation of 0.7. In this case, the covariance matrix, Σ, is not equal to the variance matrix, V, but rather given by Σ ) V1/2FV1/2, where F is the correlation matrix:
[ ]
1 0 0 0 F) 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0
0 0 0 0 0 0 1 0.7
0 0 0 0 0 0 0.7 1
Following the procedure proposed in the last section, the results of Table 16 are obtained. The optimal structure is the same for all the problem dimensions: d ) {y3, y4} ) {1, 0}. For n ) 8, the HSS technique and the formula SC5,1 are equally efficient (see Figure 6), since both techniques require approximately the same number of points to provide a solution estimate within a 0.5% band error (300 points, in the case of the HSS technique, and 272 points, in the case of the cubature formula). It is interesting to note that, in this particular case, the cubature formula SC3 is more accurate than the higher degree formulas SC5,1 and SC7. To understand this fact, one should keep in mind that the integration formulas presented here are interpolative formulas; therefore, a higher degree formula does not necessarily give a better estimate than a lower degree one. The dimension effect on the accuracy of formula SC5,1 is once more evaluated comparing the results of this formula with the ones obtained using the product Gauss rule SPG5. The relative error is only 0.1% for dimensions 3 and 4 and 0.2% for dimensions 5 and 6, showing that the dimension effect is once again negligible. Let us now consider the uncertainty model of Table 17 in which the demand for C and its price are uniformly distributed, with a positive correlation of 0.7. In this
Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 3067
Acknowledgment We acknowledge all the support (including unpublished materials and code) provided by Urmila Diwekar (Carnegie Mellon University). Literature Cited
Figure 6. HSS vs cubatures (uncertainty model 1, n ) 8).
case, a specialized formula cannot be used and the product cubature formula with the minimum number of points has 1050 points (using the formula PG9 to integrate the first two dimensions and the formula SC1 to integrate the other dimensions). According to the heuristic procedure defined in section 5, the HSS sampling technique will probably be the best choice. In fact, No ) 200 observations are enough to converge the optimal solution to within a 0.5% band error. The CPU time for this number of observations is 330.3 s and the optimal expected profit is 3643 units. The corresponding optimal design is d ) {y3, y4} ) {1, 0}. 7. Conclusions This paper presents a computational study regarding the two-stage stochastic optimization approach to process design and planning problems under uncertainty. A suitable integration technique to compute the multiple integrals involved in stochastic optimization problems is identified from a set of options that includes three integration formulas (a product Gauss rule, a specialized product Gauss rule, and a specialized cubature formula) and an efficient stratified sampling technique (HSS technique). The main criterion used to select among such alternatives is the computational effort required to achieve satisfactory accuracy, which is directly related to the number of points used by each technique. For the particular case of all n uncertain parameters being normally distributed, the most suitable technique is identified as a function of n. In the general case, a product cubature formula is constructed using an efficient specialized rule to integrate over a subset of n1 uncertain parameters that are normally distributed. The choice between this integration formula and the HSS technique is made based on the values of n and n1. In both cases, we have identified integration formulas that can be more efficient than the product Gauss formula and the HSS sampling technique, providing significant computational savings, when the problem dimension is not very high (n < 10). These results are mainly due to the introduction of a specialized cubature formula, suitable to integrate normally distributed uncertainties. We have also identified classes of problems where the HSS sampling technique is more efficient than integration formulas. Two process-engineering examples (a planning problem for a refinery plant and a process synthesis/planning problem) have been provided to illustrate and discuss these computational trends.
(1) Grossmann, I. E.; Sargent, R. W. H. Optimum Design of Chemical Plants with Uncertain Parameters. AIChE J. 1978, 37, 517. (2) Halemane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AIChE J. 1993, 29, 425. (3) Varvarezos, D. K.; Grossmann, I. E.; Biegler, L. T. An OuterApproximation Method for Multiperiod Design Optimization. Ind. Eng. Chem. Res. 1992, 31, 1466. (4) Subrahmanyam, S.; Pekny, J. F.; Reklaitis, G. V. Design of Batch Chemical Plants under Market Uncertainty. Ind. Eng. Chem. Res. 1994, 33, 2688. (5) Ahmed, S.; Sahinidis, N. V. Robust Design Planning under Uncertainty. Ind. Eng. Chem. Res. 1998, 37, 1883. (6) Reinhart, H. J.; Rippin, D. W. T. Design of Flexible Batch Chemical Plants. Presented at the Annual AIChE Meeting, New Orleans, LA, 1986. (7) Pai, C. C. D.; Hughes, R. R. Strategies for Formulation and Solving Two-Stage Problems for Process Design under Uncertainty. Comput. Chem. Eng. 1987, 11, 695. (8) Diwekar, U. M.; Rubin, E. S. Parameter Design Methodology for Chemical Processes Using a Simulator. Ind. Eng. Chem. Res. 1994, 33, 292. (9) Pistikopoulos, E. N.; Ierapetritou, M. G. Novel Approach for Optimal Process Design under Uncertainty. Comput. Chem. Eng. 1995, 19, 1089. (10) Clay, R. L.; Grossmann, I. E. A Disaggregation Algorithm for the Optimization of Stochastic Planning Models. Comput. Chem. Eng. 1997, 21, 751. (11) Petkov, S. B.; Maranas, C. D. Multiperiod Planning and Scheduling of Multiproduct Batch Plants under Demand Uncertainty. Ind. Eng. Chem. Res. 1997, 36, 4864. (12) Bernardo, F. P.; Saraiva, P. M. A Robust Optimization Framework for Process Parameter and Tolerance Design. AIChE J. 1998, 44, 2007. (13) Acevedo, J.; Pistikopoulos, E. N. A Parametric MINLP Algorithm for Process Synthesis Problems under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 147. (14) Acevedo, J.; Pistikopoulos, E. N. A Multiparametric Programming Approach for Linear Process Engineering Problems under Uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717. (15) Pertsinidis, A.; Grossmann, I. E.; McRae, G. J. Parametric Optimization of MILP Programs and a Framework for the Parametric Optimization of MINLPs. Comput. Chem. Eng. 1998, 22 (Suppl.), S205. (16) Papalexandri, K. P.; Dimkou, T. I. A Parametric MixedInteger Optimization Algorithm for Multiobjective Engineering Problems Involving Discrete Decisions. Ind. Eng. Chem. Res. 1998, 37, 1866. (17) Pistikopoulos, E. N. Uncertainty in Process Design and Operations. Comput. Chem. Eng. 1995, 19 (Suppl.), S553. (18) Ierapetritou, M. G.; Pistikopoulos, E. N. Global Optimization for Stochastic Planning, Scheduling and Design Problems. In Global Optimization in Engineering Design; Grossmann, I. E., Ed.; Kluwer Academic Publ.: Boston, 1996; Chapter 8. (19) Bernardo, F. P.; Pistikopoulos, E. N.; Saraiva, P. M. Robustness Criteria in Process Design Optimization under Uncertainty. Comput. Chem. Eng. 1999, 23, (Suppl), S459. (20) Tatang, M. A.; Pan, W.; Prinn, R. G.; McRae, G. J. An Efficient Method for Parametric Uncertainty Analysis of Numerical Geophysical Models. J. Geophys. Res. 1997, 102, 21925. (21) Diwekar, U. M.; Kalagnanam, J. R. Efficient Sampling Technique for Optimization under Uncertainty. AIChE J. 1997, 43, 440. (22) Diwekar, U. M.; Kalagnanam, J. R. An Efficient Sampling Technique for Off-line Quality Control. Technometrics 1997, 39, 308.
3068 Ind. Eng. Chem. Res., Vol. 38, No. 8, 1999 (23) Stroud, A. H. Approximate Calculation of Multiple Integrals; Prentice Hall: London, 1971. (24) Engels, H. Numerical Quadrature and Cubature; Academic Press: London, 1980. (25) Davis, P. J.; Rabinowitz, P. Methods of Numerical Integration; Academic Press: New York, 1975. (26) Carnahan, B.; Luther, H. A.; Wilkes, J. O. Applied Numerical Methods; John Wiley & Sons: New York, 1969. (27) Genz, A. C.; Malik, A. A. An Imbedded Family of Fully Symmetric Numerical Integration Rules. SIAM J. Numer. Anal. 1983, 20, 580. (28) Johnson, R. A.; Wichern, D. W. Applied Multivariate Statistical Analysis; Prentice Hall: London, 1988. (29) Epperly, T. G. W.; Ierapetritou, M. G.; Pistikopoulos, E. N. On the Global and Efficient Solution of Stochastic Batch Plant Design Problems. Comput. Chem. Eng. 1997, 21, 1411. (30) Tørvi, H.; Hertzberg, T. Estimation of Uncertainty in Dynamic Simulation Results. Comput. Chem. Eng. 1997, 21 (Suppl.), S181.
(31) Ahmed, S.; Sahinidis, N. V.; Pistikopoulos, E. N. An Improved Decomposition Algorithm for Optimization under Uncertainty. Submitted to Comput. Chem. Eng. 1998. (32) Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw Hill: New York, 1988. (33) Brooke, A.; Kendrick, D.; Meerans, A. GAMS: A User’s Guide, Release 2.25; The Scientific Press: Redwood City, CA, 1992. (34) Acevedo, J.; Pistikopoulos, E. N. A Hybrid Parametric/ Stochastic Programming Approach for Mixed-Integer Linear Problems under Uncertainty. Ind. Eng. Chem. Res. 1997, 36, 2262. (35) Sahinidis, N. V.; Grossmann, I. E.; Fornari, R. E.; Chathrathi, M. Optimization Model for Long-Range Planning in Chemical Industry. Comput. Chem. Eng. 1990, 9, 967.
Received for review November 6, 1998 Revised manuscript received April 13, 1999 Accepted April 13, 1999 IE9807001
