Mixed-Integer Nonlinear Programming Problem Process Synthesis

2680

Ind. Eng. Chem. Res. 1999, 38, 2680-2698

Mixed-Integer Nonlinear Programming Problem Process Synthesis under Uncertainty by Reduced Dimensional Stochastic Optimization Zorka Novak and Zdravko Kravanja* Faculty of Chemistry and Chemical Engineering, University of Maribor, Smetanova 17, P.O.Box 219, 2000 Maribor, Slovenia

This paper presents the developments of a novel simultaneous method for the optimization and the synthesis of complex chemical problems under uncertainty with a fixed degree of flexibility. The approach approximates the stochastic method in using a weighted objective function calculated over a reduced set of the extreme points (vertices). The feasibility of the design is ensured simultaneously by the feasibility constraints at critical vertices. The main part of the proposed method which was called the method for reduced dimensional stochastic optimization (the RDS method) is a special setup procedure for determining the reduced set of vertices and their weights for the approximation of the expected value of the objective function. A very attractive feature of this method is that the sizes of mathematical models and the computational times are reduced by 1 or 2 orders of magnitude when compared to those of the stochastic methods. On the basis of the RDS method, a robust strategy for the synthesis of complex models under uncertainty is proposed (the RDS strategy). The steps of the RDS method and the RDS synthesis strategy are illustrated in detail by one simple and two more extensive example problems, respectively. Introduction

thesis under uncertainty can be given in the following form:

Process synthesis at nominal (certain) conditions is quite well-known and has been extensively studied in numerous papers. However, in practice several pieces of data are not known exactly because they can vary in some range of values, for example, feed conditions (flow rate, composition, and temperature), cost data, product demands, and different parameters in the process (for example, heat-transfer coefficients, kinetic data, the efficiency of the catalyst, etc.). The ranges of possible values of the uncertain parameters can be forecasted by making use of previous experiences and good engineering practices. Furthermore, the probability distribution functions can also be estimated by means of which the probabilities of different combinations of the values of uncertain parameters can be calculated. The main concern of process engineers is to determine the structure of the process on one hand and the values of the operating and design variables on the other hand and to produce a design which would operate with the optimal value of economic criteria (for example, maximal profit, minimal cost, etc.) even if the values of some parameters change during the operation. Numerous approaches and strategies have been developed for designing, planning, and optimizing under uncertain conditions, while the problem of process synthesis under uncertainty has not yet been sufficiently defined. Despite the good progress which has been reached in the area in the past decade and the great practical importance of flexibility, a more general algorithmic strategy for the flexible synthesis of complex industrial problems has not yet been successfully developed. A general representation of the mixed-integer nonlinear programming problem (MINLP) of process syn-

max P(y,x,d,θ) y,x,d

s.t.

h(y,x,d,θ) ) 0

(P)

g(y,x,d,θ) e 0 x ∈ X, d ∈ D, y ∈ {0,1}m, θ ∈ TH where P represents the economic objective function (for example, profit), y represents the vector of 0-1 binary variables used to denote the existence or nonexistence of process units in the design, x represents the vector of continuous control and state variables, d is the vector of design variables, θ is the vector of uncertain parameters, and h and g are the vectors of process equations and inequalities that describe the problem mathematically. Note that the problem (P) corresponds to a MINLP infinite program since it is defined over an infinite set of values θ ∈TH. Few papers discuss the problem of process synthesis under uncertainty. Floudas and Grossmann1 presented a systematic procedure for the synthesis of flexible heatexchanger networks with uncertain flow rates and temperatures. Paules and Floudas2 presented a stochastic programming approach with MINLP recourse for the synthesis of heat-integrated distillation sequences. Papalexandri and Pistikopoulos3 proposed a multiperiod MINLP model for the synthesis of flexible heat- and mass-exchange networks. Pistikopoulos and Ierapetritou4 presented a process synthesis/planning example which was solved by means of the proposed approach for process design under uncertainty. Acevedo and Pistikopoulos5-7 presented algorithms for the

10.1021/ie980629z CCC: $18.00 © 1999 American Chemical Society Published on Web 05/28/1999

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2681

(multi)parametric solutions of MINLP and mixedinteger linear programming (MILP) models in process synthesis under uncertainty. Stability criteria was introduced at the synthesis/design stage of dynamic systems under uncertainty by Mohideen et al.8 Dua and Pistikopoulos9 described the techniques for process synthesis and material design under uncertainty. The problem of process synthesis under uncertainty is regarded as very difficult. However, the problems of process design, optimization, and planning under uncertainty are well-defined and numerous methods and approaches developed. In the last 20 years a lot of work has been done regarding the problem of process optimization under uncertainty. Several approaches and methods have been reported in the literature and it is beyond the scope of this paper to mention all of them. Grossmann et al.10 gave an overview of the optimization strategies for designing chemical processes with a fixed degree of flexibility and with an optimal degree of flexibility using rational methods. In the following years two main approaches for considering flexibility emerged (Grossmann and Straub11): the multiperiod deterministic approach and the two-stage stochastic approach: (a) Deterministic Approach. This approach assumes that the parameters vary within bounded ranges of values and approximates the original problem by using a number of uncertain parameter realizations known in advance; therefore, it is often called the “multiperiod” problem. It assumes that each uncertain parameter θ is specified by a nominal value θN and by a maximal positive and a maximal negative deviation from the nominal point which define the lower bound θLO and the upper bound θUP. The optimal design problem can then be approximated as follows: N

max d,xi

P(xi,d,θi) ∑ i)1

hi(xi,d,θi) ) 0 s.t. gi(xi,d,θi) e 0

i ) 1, ..., N

}

(P-D)

xi ∈ X, d ∈ D, θLO e θi e θUP where N is the number of uncertain parameter realizations (periods) and NP is the number of uncertain parameters. θi is the vector of uncertain parameters in the ith period. Petracci et al.12 studied the flexibility of an ethylene plant by using the deterministic approach. Hoch and Eliceche13 presented a sizing strategy for distillation columns at the design stage based on this approach. (b) Stochastic Approach. This approach uses the probabilistic distribution functions of the uncertain parameters θ to determine the design that maximizes (or minimizes) the expected value of the objective function. This problem is usually transformed into the deterministic one by means of the discretization of the uncertain parameter values followed by the approximation of the multiple integral over the feasible region for the expected value evaluation. Many authors have approximated this multiple integral through the Gaussian quadrature formula, for example, Pistikopoulos and Ierapetritou.4

N

EP ) max d,xi

piP(xi,d,θi) ∑ i)1

hi(xi,d,θi) ) 0 s.t. gi(xi,d,θi) e 0

i ) 1, ..., N

}

(P-S)

xi ∈ X, d ∈ D, θLO e θi e θUP where EP presents the expected value of the objective function and pi are the probabilities corresponding to the quadrature points (i ) 1, ..., N) regarding the joint distribution function of the uncertain parameters. Problem (P-S) is usually solved by a two-stage design strategy, where a design is selected at the first stage and the optimal vector of the continuous variables for every possible realization of the uncertain parameters at the second stage. Note that the mathematical model (P-D) is a special form of a more general (P-S) model where the probabilities of all realizations of uncertain parameters are equal. Many authors have developed methods and strategies based on the stochastic approach; for example, Straub and Grossmann14 developed a quantitative measure for the flexibility of a design denoted as the expected stochastic flexibility, Ierapetritou et al.15 presented an approach for process optimization problems involving deterministic and stochastic parameters, Ierapetritou and Pistikopoulos16 presented a framework to address the problem of batch plant design and operations under uncertainty, Liu and Sahinidis17 developed a two-stage stochastic approach to the problem of planning in process industries, Epperly et al.18 presented two global optimization algorithms for the solution of large-scale batch plants, and Ostrovsky et al.19 presented an algorithm for calculation of the flexibility function of complex systems. Despite all this, it is more or less impossible to solve complex engineering problems by using algorithmic stochastic methods because the number of quadrature points grows exponentially with the number of uncertain parameters.

Problem Formulation In this paper we attempt to develop a robust strategy for the MINLP synthesis of large-scale engineering problems under uncertainty for a fixed degree of flexibility. The objective is to optimize the structure (discrete decisions) and parameters (continuous decisions) simultaneously by considering the uncertainty of parameters, thus producing designs which satisfy all specifications and restrictions in the feasible region at the optimal value of the selected economic criteria. The main idea of this paper is to evolve the strategy toward the stochastic approach that is by no means attainable, but can successfully be approached by the proposed strategy. The following assumptions and requirements have been applied: (1) Our approach should handle the process topology, the design variables, and the operating variables simultaneously. (2) The uncertain parameters are assumed to be mutually independent.

2682 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999

(3) It is assumed that the feasible region coincides with the whole uncertainty space at a fixed degree of flexibility. (4) Although the real complex problems are often highly nonlinear and nonconvex, the derivation of our approach will rely on the assumption that the mathematical models are convex and that only the points at the extreme vertices can be the critical ones. The flexibility conditions that guarantee the feasibility of the design are then performed simultaneously only in the critical vertices. (5) We assume that the expected value of the economic criteria can be reasonably defined over a reduced set of vertices which we called basic vertices, rather than over an extensive set of quadrature points. This paper is organized in the following manner. First, the formulation of our method for the reduced dimensional stochastic optimization (RDS method) is formally introduced. An illustrating example of the optimization under uncertainty is then solved with the stochastic method and with the proposed RDS method. Then, on the basis of this method, a RDS strategy is derived for the MINLP synthesis of complex problems under uncertainty and illustrated with two examples: (i) the synthesis of a flexible heat-integrated heatexchanger network (HEN) and (ii) the synthesis of a flexible heat-integrated distillation sequence with its HEN.

Figure 1. Normal probability density function of uncertain parameters θ.

Figure 2. Probability density function of the economic objective function P(θ).

Reduced Dimensional Stochastic Optimization The deterministic and stochastic methods are based on utilizing numerous quadrature points and cannot directly be used in process synthesis in which the mathematical models are very extensive and highly nonlinear. Therefore, there is a need to develop a simplified approach for process synthesis under uncertainty. The main idea of our strategy is if the stochastic approach cannot be carried out, then it should be approximated through a simplified linear correlation which estimates the expected value of the objective function using only few basic vertices instead of a large number of quadrature points. In this way, the probabilities of the quadrature points should be replaced by weights of the basic vertices. Here, three questions arise: (1) how to decide which vertices to use (the number of vertices should be as small as possible to facilitate the computational effort, but the larger number of points may guarantee a more accurate result and feasibility over the whole region); (2) how to determine the corresponding weights of the vertices to approximate the expected value as close as possible to the true stochastic expected value; (3) how to ensure the feasibility of the design for all possible realizations of uncertain parameters. In the following subsections a novel approach to approximate the stochastic methods will be developed. The ideas of this approach will be introduced on a small one-dimensional (1-D) example, followed by the twodimensional (2-D) example, and finally expanded to an example of higher dimensions. The steps of the proposed approach will be illustrated in detail through the HEN optimization problem. (a) 1-D Example. Let us start with a simple 1-D example where the continuous uncertain parameter θ

is described by a normal distribution function (Figure 1) with the nominal value θN and lower and upper bounds (θLO and θUP) which can be defined by using σ-bounds (for example, θN ( 3σ). Let us also denote P(θLO) and P(θUP) as the optimal values of the objective function at the lower and upper bounds of the uncertain parameter, respectively (Figure 2). The evaluation of the expected value of the objective function EP is not straightforward, because EP * P(θN). The idea here is to transfer the distribution function of θ to the distribution of the objective function P(θ). EP can then be represented by the integral over the feasible region of θ (Figure 2) which is usually approximated with some numerical method, for example, with the rectangular rule, the trapezoidal rule, the Simpson’s rule, or the Gaussian integration formula: ND

EP )

∫θp(θ) P(θ) dθ ≈ ∑vjP(θj)

(1)

j)1

where θj represents discretized points of uncertain parameters (e.g., in the Gaussian integration formula θj corresponds to the zeros of the Legendre poyinomial) and vj coefficients determined by the numerical method. The approximation of EP in (1) requires ND sequential optimizations of the original problem at different θj values between the lower and upper bounds. The selection of ND points can in general be equidistant or nonequidistant, depending on the method used for numerical integration. It should be noted that a higher number of points involves a more accurate estimation of the expected value EP, but it also increases the number of required sequential optimizations of the original problem.


On the other hand, EP can be approximated in the 1-D problem only by using the extreme values:

EP ) w1P(θLO) + w2P(θUP)

(2a)

w1 + w2 ) 1

(2b)

where w1 and w2 represent the positive weights corresponding to the P values at the lower and upper values of the uncertain parameter. By combining eqs 1, 2a, and 2b, the following set of two linear equations can be solved for w1 and w2. ND

vjP(θj) ) w1P(θLO) + w2P(θUP) ∑ j)1

Figure 3. Problem with two uncertain parameters.

(3)

w1 + w2 ) 1 When w1 and w2 are known, the value of EP in the simultaneous optimization procedure can simply be approximated by the linear combination of the objective function values P at the lower and upper bounds of θ (eq 2a). (b) 2-D Example. We shall now expand our reflections to a 2-D problem with two uncertain parameters, θ1 and θ2, described by normal distribution functions. Again, since the expected value of the objective function in general is not equal to the optimal value of the objective function at the nominal conditions, there is a need, in simultaneous optimization, to compute the adequate EP rather than to perform an optimization at the nominal value θN. Assuming the independence of the uncertain parameters, the expected values of the objective function regarding each uncertain parameter separately, EP(θ1) and EP(θ2), can also be treated as independent and can be obtained, for example, through Gaussian integration: ND

EP (θ1) )

∑ j)1

v1,jP(θ1,j,θN 2 ) and ND

v2,jP(θN ∑ 1 ,θ2,j) j)1

[ ] [

] [ ] [ ] [

N N EP(θ1) P(θLO P(θLO 1 ,θ2 ) 1 ,θ2 ) ) w1 + w 2 LO UP + EP(θ2) P(θN P(θN 1 ,θ2 ) 1 ,θ2 )

N P(θ1,j,θN 2 ) ) max P(xj,dj,θ1,j,θ2 ) xj,dj

h(xj,dj,θ1,j,θN 2) ) 0 g(xj,dj,θ1,j,θN 2) e 0

N N P(θUP P(θUP 1 ,θ2 ) 1 ,θ2 ) N LO + w4 N UP P(θ1 ,θ2 ) P(θ1 ,θ2 )

]

4

(4)

where vk,j represents the probability of point θk,j according to the distribution function of θk. The values P(θ1,j, θN 2 ) are obtained through the sequential solutions of the problems (P1)j at selected θ1,j points, while the other uncertain parameter is kept at its nominal value θN 2:

xj ∈ X, dj ∈ D

is to express EP(θ1) and EP(θ2) by using a linear combination of the objective function values calculated at the extreme points of each uncertain parameter:

w3

EP (θ2) )

s.t.

Figure 4. Values of the objective function in the 2-D example.

(P1)j

j ) 1, ..., ND

N P(θN 1 ,θ2,j) can be obtained similarly at the fixed θ1 . Values of EP(θ1) and EP(θ2) can then be calculated by using (4). For flexibility reasons, Swaney and Grossmann20 demonstrated that under convexity conditions only combinations of the extreme vertices of the uncertain parameters need to be analyzed. The set of extreme vertices is then i ) 1, ..., N ) 2NP. The idea in this paper

wi ) 1 ∑ i)1

(5)

wi, i ) 1, ..., 4, in eq 5 plays the role of the fractions by which the corresponding vertices contribute to the EP(θk); therefore, their sum is equal to 1. Alternatively, the system (5) can also be obtained when the deviation N of EP(θk) from P(θN 1 ,θ2 ) is expressed as a linear combination of deviations of four extreme vertices from the nominal point as shown in Appendix A. The procedure is illustrated in Figure 3 which represents the situation in θ space. To determine EP(θ1) and EP(θ2), we fix one dimension (θ2 and θ1, respectively) and optimize the problem by changing another parameter (θ1 and θ2, respectively) in positive and negative directions. This gives rise to the diamond shape inside the square in θ1 and θ2 space. The corners of the N diamond are used to determine the points P(θLO 1 ,θ2 ), UP N N LO N UP P(θ1 ,θ2 ), P(θ1 ,θ2 ), and P(θ1 ,θ2 ) (Figure 4). The sides of the diamond represent the directions of the linear combinations which are used in eq 5. Under the assumption of the independence of uncertain parameters, the set in eq 5 provides weights for the approximation of the expected value (EP) by the linear combination of P values at four extreme vertices

2684 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999

which are obtained by combining the extreme values of uncertain parameters:

2NP

EP ) max

∑wi P(xi,d,θk,i|

k)1,...,NP

d,xi i)1 LO LO UP EP ) w1P(θLO 1 ,θ2 ) + w2P(θ1 ,θ2 ), + LO UP UP w3P(θUP 1 ,θ2 ) + w4P(θ1 ,θ2 ) (6)

(c) Extension to Higher Dimensions. The system of equations in (5) can now be written in the following general form:

∑ i)1

}

(P-RDS)

xi ∈ X, d ∈ D where

θk,i )

{

2NP

EP(θk) )

hi(xi,d,θk,i) ) 0 gi(xi,d,θk,i) e 0 i ) 1, ..., 2NP d g fi(xi,θk,i)

s.t.

)

NP-k θLO < i e (2n - 1) × 2NP-k k , if (2n - 2) × 2 NP-k θUP < i e (2n) × 2NP-k k , if (2n - 1) × 2

wiP(θk,i,θN kk|kk)1,...,NP and kk*k)

where

k ) 1, ..., NP; n ) 1, ..., 2k-1; i ) 1, ..., 2NP

θk,i )

Note that for the purpose of our work, design variable d in the problem (P-RDS) is expressed explicitly by using the vector of the design function f. The design expressions are written as nonequalities in order to provide feasibility of the design variables simultaneously in all selected points. However, design variables can also appear in the (non)equalities g and h, as for example in the cost functions. The important advantage of the (P-RDS) model compared to the stochastic model (P-S) is that it has to be solved simultaneously only at 2NP extreme points (vertices), instead of at NDNP quadrature points, thus multiplying the size of the original problem “only” by 2NP. Although this fact significantly facilitates the solution of the problem, the number of vertices in more complex problems is still too large to be able to solve the optimization problems. Because the set of equations in (7) has a positive degree of freedom, the idea then is to find a solution with the minimum number of nonzero weights, thus minimizing the number of needed vertices. The reduction can be achieved in two steps: first, by determining the critical vertices from a complete set of vertices to guarantee the feasibility of the solution, and second, by performing a special MILP problem to define the minimum number of additional vertices that are needed for the approximation of the expected value. (d) Determination of the Critical Vertices. The formal definition of the critical point is given for a fixed d and y by the following minimax problem:

{

NP-k θLO < i e (2n - 1) × 2NP-k k , if (2n - 2) × 2 NP-k UP θk , if (2n - 1) × 2 < i e (2n) × 2NP-k

k ) 1, ..., NP; n ) 1, ..., 2k-1; i ) 1, ..., 2NP

(7)

and 2NP

wi ) 1 ∑ i)1 0 e wi e 1

i ) 1, ..., 2NP

Note that in the upper model the vector θ is decomposed to the components (θk). P(θk,i,θN kk) represents the value of the objective function in the ith vertex obtained at the lower or upper value of θk, denoted as θk,i, while the vector of the other uncertain parameters is kept at the nominal value (θN kk). The system in (7) for three uncertain parameters is presented in Appendix B. To set up the system in (7) in addition to P values, we also need information about the expected value of the objective function for each uncertain parameter EP(θk), k ) 1, ..., NP which is calculated as follows: ND

EP(θk) ) θUP k

vk,jP(θk,j, θN ∑ kk| j)1

e θk,j e

θLO k

kk)1,...,NP and kk*k

j ) 1, ..., ND

)

max min u θ

(8)

x

h(y,x,d,θ) ) 0

(9)

g(y,x,d,θ) e u P(θk,j,θN kk)

where the values are obtained through the ND sequential optimizations of the original problem in UP quadrature points θk,j between θLO k and θk , while the remaining uncertain parameters are kept at the nominal values. It is important to note that the number of evaluations increases linearly with the number of uncertain parameters NP. Hence, the original problem has to be solved sequentially for NP × ND times to estimate EP(θk) for k ) 1, ..., NP. ND is usually 5. The solution of the system (7 and 8) gives the weights wi. They are then used for the approximation of the expected objective function (EP) as a linear combination of objective functions (P) at vertices:

Halemane and Grossmann21 determined the critical vertices in the iterative multiperiod algorithm. For complex and large-scale problems we preferred to make use of the fact that some of the uncertain parameters affect the design variables monotonically. Such parameters can be recognized with a small computing effort using a simple procedure by which the influence of the change in the uncertain parameter on the design variables is examined. This actually means analyzing the signs of gradients ∂d/∂θ. Uncertain parameters with equal signs of gradients for all design variables can be regarded as monotonic parameters and can be considered in the optimization only at their critical valuess


at either the upper or the lower bound. If the number of such parameters is NM, the number of extreme points is then reduced to 2NP-NM. Another option, which will be presented in detail in the following illustrating example, is to solve the problem sequentially in all vertices by using the automated computational procedure. The vertices that require the largest value for at least one design variable are then determined as critical. Because the critical vertices are selected heuristically, a feasibility/flexibility test is necessary to check the actual flexibility. (e) MILP Problem for Minimizing the Number of Vertices. The expected value of the objective function EP in the problem (P-RDS) is expressed as a linear combination of the objective function values of P at all vertices. The idea now is to express EP by using a reduced set of vertices. The reduced set should include as many critical vertices as possible because they have to be included in the optimization anyway to ensure the feasibility of the design. The selected critical vertices together with the minimal number of additional noncritical ones comprise the reduced set of vertices for the approximation of the expected value of the objective function. A special MILP problem will be introduced below for the selection of the reduced set of points together with the corresponding weights. The following notations will be used in the formulating of the MILP problem:

and

wi ) 1 ∑ i∈I

The original set of the extreme points I should be divided into two subsets: the subset of critical vertices (CV) and the subset of noncritical vertices (I\CV). The following MILP problem consists of k linear equality constraints (11) in which the weights wi are unknown. The expected values of the objective function EP(θk) for each uncertain parameter in (11) are calculated using (8). Weights wi are defined as positive variables, the sum of which is equal to 1 (12). y is the vector of the binary variables associated with the noncritical vertices to denote the selection (yi ) 1) or rejection (yi ) 0) of the vertex. The sum of binary variables for noncritical vertices is minimized (10) in order to ensure the minimal number of noncritical vertices which are required besides the critical ones. The logical constraint (13) guarantees that the weight becomes zero if the corresponding noncritical vertex is rejected.

min

∑

yi

(10)

i∈I\CV

EP(θk) )

wiP(θk,i,θN ∑ kk| i∈I

kk)1,...,NP and kk*k

)

(11)

where

θk,i )

{

NP-k θLO < i e (2n - 1) × 2NP-k k , if (2n - 2) × 2 NP-k θUP < i e (2n) × 2NP-k k , if (2n - 1) × 2

k ) 1, ..., NP; n ) 1, ..., 2k-1; i ) 1, ..., 2NP

(P-W)

i ∈ I\CV

wi e y i 0 e wi e 1

i∈I

(13)

As can be seen from the above (P-W) problem, the solution yields the values of the weights wi (i ) 1, ..., 2NP). Only vertices with nonzero weights, denoted as basic vertices (BV), are then used for the calculation of the expected value of the objective function. (f) Reduced Model for the Reduced Dimensional Stochastic Optimization. To ensure the feasibility of the design and to approximate the expected value at the same time, the optimization model is defined over the union of the critical vertices and basic vertices (CV ∪ BV). Thus, the union (CV ∪ BV) contains all critical vertices CV and the selected noncritical ones that have yi ) 1. The reduced P-RDS model can now be formulated as

EP ) max

∑ wiP(xi,d,θk,i)

d,xi i∈BV

hi(xi,d,θk,i) ) 0 gi(xi,d,θk,i) e 0 i ∈ CV ∪ BV d g fi(xi,θk,i)

Sets K ) set of uncertain parameters, (k|k ) 1, ..., NP) I ) set of vertices, (i|i ) 1, ..., 2NP) CV ) subset of critical vertices, (cv|cv ) critical vertex defined by the user) I\CV ) subset of noncritical vertices, (i and i * cv)

(12)

s.t.

}

(P-RDS)

xi ∈ X, d ∈ D where

θk,i )

{

NP-k θLO < i e (2n - 1) × 2NP-k k , if (2n - 2) × 2 UP θk , if (2n - 1) × 2NP-k < i e (2n) × 2NP-k

k ) 1, ..., NP; n ) 1, ..., 2k-1; i ) 1, ..., 2NP It should be noted that the number of vertices in the set BV is found to be NP + 1. If the number of critical vertices in the set CV can be reduced to 2NP-NM, then the maximal number of vertices in the union CV ∪ BV is equal to 2NP-NM + NP + 1. For example, if the problem has six uncertain parameters (NP ) 6) from which three monotonically affect the design variables (NM ) 3), then the size of the (P-RDS) model would at the most be 15 times the size of the original optimization problem. In our examples it was found that usually more than half of the vertices in the set BV are critical ones because the (P-W) problem favors them. The actual size of the union CV ∪ BV is always less than 2NP-NM + NP + 1. In the above-mentioned example with NP ) 6 and NM ) 3 we have obtained a union CV ∪ BV with 11 vertices, which clearly indicates that the computational effort is significantly reduced compared to the stochastic approach. (g) Summarized Procedure for Reduced Dimensional Stochastic Optimization. The procedure described in the previous sections can be summarized as follows:

2686 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 Step 1: Set up the problem. (a) For each uncertain parameter, solve the original optimization problem sequentially in ND points, while keeping the other uncertain parameters at their nominal values. Calculate the EP(θk) for uncertain parameters by using one of the numerical integration formulas, for example, the Gaussian one, and the suitably chosen coefficients vk,j (eq 8). If the extreme points of uncertain parameters do not coincide with one of the selected ND points (the integration formula is an UP open formula), solve the original problem also at the extreme points for all uncertain parameters to obtain P(θLO k ) and P(θk ). (b) Define the set of extreme vertices (I), determine the subsets of critical (CV) and noncritical vertices (I\CV), and formulate the (P-W) problem. The solution will give the set of basic vertices (BV) and the corresponding weights for the approximation of the expected value. Step 2: Solve the P-RDS optimization problem simultaneously over the reduced set of vertices (CV ∪ BV) by using the expected objective function which is approximated by the linear combination of the objective functions at basic vertices BV. Step 3: If necessary, check the actual flexibility of the optimal solution. If the solution is not flexible, determine new critical vertices and go to step 1b.

The Illustrating Example Statement To illustrate the above-mentioned procedure, an example of a heat-exchanger network optimization will be presented. The example problem will first be solved with the stochastic method and then with the reduced dimensional stochastic method so that the results can be compared. The heat-exchanger network consists of two cold streams, two hot streams, and cold utility (Biegler et al.).22 We have extended this example with a steam heater (unit 5) to increase the transparency of the problem (Figure 5). The uncertain parameters of our example are described by normal distribution functions using the form θ N(θN,σ): T3 N(388,3.33) K, T5 N(583,3.33) K, and T9 N(298,1.67) K. The mathematical model (P2) consists of the objective function C, heatbalance equations, temperature specifications, and the expressions for heat-exchanger area calculations. Note that in the objective function the total cost C (investment cost + operating cost) was chosen as the economic criterion due to constant feed flow rates. The heatexchangers areas are treated as design variables.

Al )

φl Ul∆lnTl

l ) 1, 2, ..., 5

U1 ) U2 ) U3 ) U4 ) 700 W/(m2 K), U5 ) 1000 W/(m2 K) Cg0 l ) 1, 2, ..., 5

T1, T2, T3, T4, T5, T6, T7, T8, T9, T10 g 0

4

∑ l)1

Heat-exchanger area calculations:

Al g 0, φl g 0

Objective function: C ) min(1846

Figure 5. Heat-exchanger network for the illustrating example.

A0.65 l

+

2350A0.65 5

+ 0.02φ4 + 0.23φ5)

Heat-balance equations: exchanger 1: φ1 ) 1.5 (620 - T2) ) 2 (T4 - T3) exchanger 2: φ2 ) T5 - T6 ) 2(T8 - T4) exchanger 3: φ3 ) T6 - T7 ) 3(393-313) exchanger 4: φ4 ) 1500(T2 - 350) T10 ) φ4/6000 + T9 exchanger 5: φ5 ) 2000(600 - T8) Temperature specifications: exchanger 1: T2 g T3 + 1 exchanger 2: T6 g T4 + 1 T5 g T8 + 1 exchanger 3: T7 g 313 + 1 T6 g 393 + 1 T7 e 323 exchanger 4: T10 e 323 exchanger 5: 619 g T8 + 1

(P2)

The problem (P2) was solved at the nominal conditions (T3, T5, T9) ) (388, 583, 298) K yielding the total cost of 44230 $/year and the optimal values of the design variables (A1, A2, A3, A4, A5) ) (15.50, 2.12, 6.08, 1.53, 2.09) m2. It is quite easy to test that the optimal design is not able to satisfy the specifications and constraints in (P2) if variations in uncertain parameter values occur. For example, a small change of the T3 value to 389 K at a fixed design leads to an infeasible solution. Therefore, the design has to be treated in such a manner that it becomes feasible for all combinations of uncertain parameters within their specified regions. The objective then is to determine the optimal oversizing of the design variables together with the optimal operating conditions. (a) Solving the Illustrating Problem with the Stochastic Method. Our illustrating example is small enough to be solved with the simultaneous stochastic method, in numerous quadrature points over a feasible region of uncertain parameters. In contrast to the strategy of Pistikopoulos and Ierapetritou,4 the discretization of the uncertain parameters described by normal distribution functions in our illustrating example was done a priori, over five points j ) 1, ..., 5 (Table 1). The points θj are calculated through the conversion of the interval [-1, 1] where the zeros of the Legendre polynomial of the fifth degree (tj) are defined in the interval (-3σ, 3σ). These bounds ensure a closed region

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2687 Table 1. Quadrature Points, Weights, and Probabilities of the Motivating Example j

tj

1 2 3 4 5

-0.906 18 -0.538 47 0.0 0.538 47 0.906 18

Bj

θj

(T3)j

(T5)j

(T9)j

v(θj)

0.236 93 0.478 63 0.568 89 0.478 63 0.236 93

- 2.718σ θN - 1.615σ θN θN + 1.615σ θN + 2.718σ

378.938 382.615 388.000 393.385 397.062

573.938 577.615 583.000 588.385 592.062

293.469 295.308 298.000 300.692 302.531

0.0070 0.1545 0.6770 0.1545 0.0070

θN

with almost 99.7% of the uncertain parameter values:

tj(b - a) + b + a θj ) 2 a ) θN - 3σ, b ) θN + 3σ

Al g (14)

The coefficients v which represent the probabilities for θj are calculated by using the Gaussian integration formula:

v(θj) )

b-a BjfG(θj) 2

j ) 1, 2, ..., 5

(15)

where Bj represents the coefficients corresponding to the points tj, fG is a density function of normal distribution for each uncertain parameter, and θj values are calculated by using eq 14. Three uncertain parameters (NP ) 3) with five values over the feasible region (ND ) 5) form 125(NDNP) different points in which the original model (P2) has to be solved simultaneously. To accomplish the task, the model (P2) was reformulated in the stochastic (P-S) form. In the reformulated stochastic (P-S) model (P3) the temperatures and heat flow rates were indexed over the 125 quadrature points m ) 1, 2, ..., 125. The expressions for the design variables were written as nonequalities to achieve feasibility in all quadrature points simultaneously. The objective function of problem (P3) is expressed as the minimized expected value of the total cost (EC). 125

EC ) min

∑ m)1

4

A0.65 + 2350A0.65 + ∑ l 5 l)1

[vm(1846

0.02φ4,m + 0.23φ5,m)] φ1,m ) 1.5(620 - T2,m) ) 2(T4,m - T3,m) φ2,m ) T5,m - T6,m ) 2(T8,m - T4,m) φ3,m ) T6,m - T7,m ) 3(393-313) φ4,m ) 1500(T2,m - 350) T10,m ) φ 4,m/6000 + T9,m φ5,m ) 2000(600 - T8,m) T2,m g T3,m + 1 T6,m g T4,m + 1 T5,m g T8,m + 1 T7,m g 313 + 1 T6,m g 393 + 1 T7,m e 323 T10,m e 323

619 g T8,m + 1 φl,m Ul ∆lnTl,m

U1 ) U2 ) U3 ) U4 ) 700 W/(m2 K), U5 ) 1000 W/(m2 K) EC g 0 Al g 0, φl,m g 0 T1,m, T2,m, T3,m, T4,m, T5,m, T6,m, T7,m, T8,m, T9,m, T10,m, g 0 m ) 1,2, ..., 125; l ) 1,2, ..., 5

(P3)

T3,m, T5,m, and T9,m are the values of the uncertain parameters in the mth point that are obtained with the combination of (T3)j, (T5)j, and (T9)j in Table 1 while the other temperatures and heat flow rates are the optimization variables. For example, for m ) 1 the quadrature point (378.938, 573.938, 293.469) K is determined with (T3)1, (T5)1, and (T9)1; for m ) 2 we obtain indices (T3)1, (T5)1, and (T9)2 and the point is (378.938, 573.938, 295.308) K, etc. Coefficient vm in the objective function represents the joint probability of the mth point. The uncertain parameters are assumed to be mutually independent; hence their individual probabilities can be multiplied to provide joint probabilities of 125 quadrature points. For example, the probability of the first quadrature point v1 ) v(T3)1v(T5)1v(T9)1 ) 0.0073 ) 3.43 × 10-7. The illustrating problem was solved simultaneously in 125 points by using the (P3) model which yielded the result 45 300 $/year and the design variables (A1, A2, A3, A4, A5) ) (15.34, 2.37, 6.32, 2.04, 2.34) m2. As expected, the majority of the design variables exceed the values obtained at nominal conditions because of the flexibility conditions. (b) Solving the Illustrating Problem by Reduced Dimensional Stochastic Optimization. In the next step, the illustrating example was solved by using the proposed RDS method described in the previous section. Step (1a). In the first step the original problem (P2) was solved sequentially at five values (j ) 1, ..., ND ) 5) for each uncertain parameter (NP ) 3) considering the flexibility conditions at eight vertices, while the other two uncertain parameters were kept at nominal values. The following values of the uncertain parameters were used:

(T3)j: 378.9, 382.6, 388, 393.4, and 397.1 K N (TN 5 ) 583 K and T9 ) 298 K)



j ) 1, ..., 5

2688 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 Table 2. Vertices in the Motivating Example vertex (i)

uncertain parameter

1

2

3

4

5

6

7

8

T3 (K) T5 (K) T9 (K)

378 (LO) 573 (LO) 293 (LO)

378 (LO) 573 (LO) 303 (UP)

378 (LO) 593 (UP) 293 (LO)

378 (LO) 593 (UP) 303 (UP)

398 (UP) 573 (LO) 293 (LO)

398 (UP) 573 (LO) 303 (UP)

398 (UP) 593 (UP) 293 (LO)

398 (UP) 593 (UP) 303 (UP)

Table 3. Values of the Design Variables at the Vertices vertex (i)

A1 (m2)

A2 (m2)

A3 (m2)

A4 (m2)

A5 (m2)

1 2 3 4 5 6 7 8

16.739 16.780 16.348 16.388 14.656 14.674 14.267 14.285

1.095 1.098 2.446 2.451 1.505 1.508 3.093 3.099

6.071 6.071 6.071 6.071 6.071 6.071 6.071 6.071

1.148 1.342 1.161 1.356 1.655 1.909 1.668 1.924

2.307 2.306 1.974 1.972 2.193 2.192 1.841 1.840

and the following values of the objective function (C) were obtained by solving problem (P2) at these points: N C(T3)j at TN 5 , T9 : 45 773, 45 492, 45 100, 44 730, and 44 494 $/year N C(T5)j at TN 3 , T9 : 46 693, 46 049, 45 100, 44 117, and 43 429 $/year N C(T9)j at TN 3 , T5 : 45 100, 45 100, 45 100, 45 100, 45 100 $/year

which were then used for the estimation of EC(θ) through the Gaussian integration formula (8) with the coefficients v (θj) from Table 1 giving the following results:

at the extreme values of one uncertain parameter and at the nominal values of the remaining two parameters. The following results are obtained: N N C(TLO 3 , T5 , T9 ) ) C(378, 583, 298) K

N N C(TUP 3 , T5 , T9 ) ) C(398, 583, 298) K

LO N C(TN 3 , T5 , T9 ) ) C(388, 573, 298) K

C(TN 3,

TUP 5 ,

TN 9)

) C(388, 593, 298) K

C(TN 3,

TN 5,

TLO 9 )

) C(388, 583, 293) K

C(TN 3,

TN 5,

TUP 9 )

) C(388, 583, 303) K

) 44 436 $/year ) 46 869 $/year ) 43 252 $/year ) 45 100 $/year ) 45 100 $/year

and rearranged in Table 4 at the corresponding extreme values of the uncertain parameters. Using the data in Table 4 and the EC data obtained in step (1a), the following MILP program can be postulated:

EC(T3) ) 45 104 $/year, EC(T5) ) 45 094 $/year, and EC(T9) ) 45 100 $/year Note that 15 successive optimizations of the problem (P2) were required in this step (the actual number of optimizations in this example is 13 because the three optimizations for the point with all temperatures at the nominal conditions need to be done only once). Step (1b). The feasible region of uncertain parameters was bounded with the use of σ bounds (θN ( 3σ). The extreme values of three uncertain parameters obtained UP LO in such a way are TLO 3 ) 378 K, T3 ) 398 K, T5 ) 573 UP LO UP K, T5 ) 593 K, T9 ) 293 K, and T9 ) 303 K. If we assume that only extreme points are taken into account, we have to consider 8(2NP) vertices, i ) 1, 2, ..., 8 (Table 2). Because no uncertain parameter has a monotonic effect on the design variables in this example, another option for the identification of critical vertices was used: the problem was solved sequentially in all vertices while the values of the design variables were examined (Table 3). The largest values of design variables appear in vertices 1, 2, and 8. Hence, these three vertices are selected as critical. The set of critical vertices is then CV ) {1, 2, 8}. The remaining vertices compose the subset of noncritical vertices I\CV ) {3, 4, 5, 6, 7}. To calculate weights by solving the (P-W) problem, we need information about the objective function values at the extreme values of the uncertain parameters. Since the Gaussian integration in step (1a) is not performed at the extreme values (see Table 1), six additional optimization problems (P2) have to be performed sequentially

) 45 845 $/year

7

min

yi ∑ i)3

s.t. EC(T3) ) w1C(T3)1 + w2C(T3)2 + ... + w8C(T3)8 EC(T5) ) w1C(T5)1 + w2C(T5)2 + ... + w8C(T5)8 EC(T9) ) w1C(T9)1 + w2C(T9)2 + ... + w8C(T9)8 8

wi ) 1 ∑ i)1 wi e yi

i ) 3, 4, 5, 6, 7

(16)

The (P-W) problem (16) is then solved for weights giving the following result:

w1 ) 0.4738, w6 ) 0.0355, w8 ) 0.4906, w2 ) w 3 ) w4)w5)w7)0 In this case, nonzero weights are assigned to critical vertices 1 and 8 and to an additional noncritical one, that is, vertex 6. The set of basic vertices is then BV ) (1, 6, 8) and their corresponding weights for the linear approximation of the expected value are w1 ) 0.4738, w6 ) 0.0355, and w8 ) 0.4906. The union CV ∪ BV includes all basic vertices and critical vertex 2, that is, CV ∪ BV ) (1, 2, 6, 8). Step 2. Note that the (P-RDS) model needs to be solved only in 4 points, and not in 125, when using the objective function with the weights calculated in step

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2689 Table 4. Cost Values for Uncertain Parameters vertex (i) C ($/year)

1

2

3

4

5

6

7

8

C(T3) C(T5) C(T9)

45 845 (LO) 46 869 (LO) 45 100 (LO)

45 845 (LO) 46 869 (LO) 45 100 (UP)

45 845 (LO) 43 252 (UP) 45 100 (LO)

45 845 (LO) 43 252 (UP) 45 100 (UP)

44 436 (UP) 46 869 (LO) 45 100 (LO)

44 436 (UP) 46 869 (LO) 45 100 (UP)

44 436 (UP) 43 252 (UP) 45 100 (LO)

44 436 (UP) 43 252 (UP) 45 100 (UP)

(1b). The original model (P2) can now be reformulated in the (P-RDS) form: 4

EC ) min

A0.65 + 2350A0.65 + ∑ [wi(1846∑ l 5 i)1,6,8 l)1 0.02φ4,i + 0.23φ5,i)]

Table 5. Model Sizes and CPU Times on the IBM RISC 6000 Model 43P single equations single variables CPU (s)

P-S

P-RDS

2397 1072 42

73 50 0.27

φ1,i ) 1.5(620 - T2,i) ) 2(T4,i - T3,i) φ2,i ) T5,i - T6,i ) 2(T8,i - T4,i) φ3,i ) T6,i - T7,i ) 3(393-313) φ4,i ) 1500(T2,i - 350) T10,i ) φ 4,i/6000 + T9,i φ5,i ) 2000(600 - T8,i) T2,i g T3,i + 1 T6,i g T4,i + 1 T5,i g T8,i + 1 T7,i g 313 + 1

(P4)

T6,i g 393 + 1 T7,i e 323 T10,i e 323 619 g T8,i + 1 Al g

φl,i Ul∆lnTl,i

U1 ) U2 ) U3 ) U4 ) 700 W/(m2 K), U5 ) 1000 W/(m2 K) w1 ) 0.4738, w6 ) 0.0355, w8 ) 0.4906 EC g 0

Figure 6. RDS strategy for process synthesis under uncertainty.

Al g 0, φl,i g 0

variables from the (P-RDS) model. The solution yields the corrected expected value EC* ) 45 932 $/year which differs from the stochastic result by about 1.4%. As can be seen, the result of P-RDS is very close to the one of P-S which has motivated us to apply the idea of the RDS approach to the problem of process synthesis. In addition, the size of the P-RDS problem is reduced by more than 1 order of magnitude and the central processing unit (CPU) time by more than 2 orders of magnitude when compared to the P-S problem (Table 5).

T1,i, T2,i, T3,i, T4,i, T5,i, T6,i, T7,i, T8,i, T9,i, T10,i, g 0 i ∈ CV ∪ BV ){1, 2, 6, 8}; l ) 1, 2, ..., 5 The solution has the expected cost of 45 970 $/year and the design variables (A1, A2, A3, A4, A5) ) (14.40, 3.14, 6.16, 1.92, 2.38) m2. Comparison between the (P-S) and the (P-RDS) Solutions. The results of the P-S (45 300 $/year) and the P-RDS (45 970 $/year) model obtained in subsections (a) and (b) differ by 1.5%. Because of the linear approximation of the objective function in the (P-RDS) problem, the resulting design variables do not correspond correctly to the cost 45 970 $/year. If we assume that the stochastic model (P-S) gives a correct result, it should be interesting to verify the cost of the (P-RDS) solution by resolving the (P-S) problem at fixed design

The Strategy for Process Synthesis under Uncertainty The ideas that were introduced on the problem of the optimization under uncertainty in the previous sections can now be expanded to the problems of process synthesis under uncertainty.

2690 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 Step 1: Definition of the Problem. Define the superstructure of the problem. Select the initial flowsheet by fixing binary variables. Define the uncertain parameters and bound their feasible region by using σ bounds. Step 2: NLP Step. (a) Try to find out if some of the uncertain parameters affect the design variables monotonically and consider those parameters either at lower or upper extreme values only. From the set of extreme points, vertices (I) define the set of critical vertices (CV). (b) Set up the problem. Sequentially solve the NLP subproblems for a fixed structure by changing the values of the kth uncertain parameter ND times and keep the remaining uncertain parameters at nominal values. Repeat the ND calculations for each uncertain parameter and calculate the vector of the expected values of the objective function for each uncertain parameter EP(θ k). If the UP integration formula is an open formula solve the NLP subproblems at the extreme points to calculate P(θLO k ) and P(θk ). Formulate and solve the (P-W) problem to determine the reduced set of vertices BV and their weights for the approximation of the expected value. (c) Solve the (P-RDS) problem for a fixed structure with the objective function expressed as the linear combination of the functions at BV vertices with weights estimated in step (2b). Define the constraints at vertices from the union CV ∪ BV to ensure also the feasibility of the design. The result is the optimal expected value of the objective function EP for the particular structure. Step 3: Check the Termination of the OA/ER Algorithm. Compare the EP value with the best one obtained so far. If the values are not improving, then go to step 5; otherwise, save the current EP value as the best one and go to step 4. Step 4: MILP Step. (a) Derive the outer approximations of the P-RDS problem and construct the MILP master problem. (b) Solve the MILP master problem to predict the new topology and go to step 2. Step 5: Check the Feasibility. Check the feasibility of the optimal solution sequentially in all quadrature points with respect to the optimal values of design variables d. If the design is not feasible, find new critical vertices and go to step 2a.

Table 6. Data for Example 1a stream H1 H2 C1 C2 CW St a

TIN (K)

TOUT (K)

CF (kW/K)

R (W/m2‚K)

cost ($/kW‚year)

443 423 293 353 290 438

333 303 398 413 303

30 15 20 40

1000 1000 800 800 800 5000

20 80

Exchanger cost ($/year) ) 5500 + 150 × area (m2). Minimum approach temperature ) 5 K.

One of the most promising methods for solving MINLP synthesis problems is the outer approximation/ equality relaxation (OA/ER) algorithm which was postulated in its earlier version by Duran and Grossmann23 and then extended by Kocis and Grossmann24 and Viswanathan and Grossmann.25 The algorithm is based on the alternating solution of a sequence of a nonlinear programming (NLP) subproblem and a mixed-integer linear programming (MILP) master subproblem. The algorithm is automated together with its extensions and improvements (convexity test) in the computer package PROSYN-MINLP synthesizer by Kravanja and Grossmann.26 These facts have encouraged us to implement the reduced dimensional stochastic optimization through the OA/ER algorithm in the PROSYN-MINLP package, thus expanding its capabilities to be able to solve process synthesis problems under uncertainty. Figure 6 presents the proposed RDS strategy. Example Problems In this section, the above-mentioned strategy is applied to two example problems: the synthesis of the flexible heat-integrated heat-exchanger network (HEN) and the synthesis of the flexible distillation sequence for the separation of a four-component mixture together with its heat-integrated HEN. These examples cannot be solved with the simultaneous stochastic method; therefore, they were solved by the use of the deterministic method. The first example was also solved by the two-stage stochastic method with Monte Carlo simulation while this was not possible with the second example. Example 1: MINLP Synthesis of a Flexible HEN. Given are two hot streams (H1, H2), two cold streams (C1, C2), steam (St), and cooling water (CW) as utilities. Table 6 gives the supply and target temperatures (TIN and TOUT) of the streams, the heat capacity flow rates (CF), the film and fouling coefficients (R), and the utility costs. Supply temperatures of the hot stream H1, cold stream C2, cooling water, and steam will be treated as stochastic uncertain parameters described by normal

Table 7. Uncertain Parameters of Example 1 uncertain parameter

distribution function

nominal value

expected deviations

TIN,H1 (K) TIN,C2 (K) TIN,CW (K) TIN,St (K) r

N (443, 3.3) N (353, 3.5) N (290, 4.0) N (438, 8.4) uniform

443 353 290 438 0.5

(9.0 (9.5 (11.0 (23.0 (0.5

Table 8. Nominal Values of Film and Fouling Coefficients and Expected Deviations film and fouling coefficients (W/m2‚K)

nominal value

expected deviations

RH1, RH2 RC1, RC2, RCW RSt

1000 800 5000

(181 (145 (1133

distribution functions as shown in Table 7. The feasible regions of uncertain temperatures are determined by using σ bounds θN ( 2.7σ to coincide the endpoints of the intervals with two quadrature points and thus to reduce the number of computations. It will be assumed that film and fouling coefficients decrease linearly during the operation with the same rate of change (r) from the upper to the lower values that are known in advance (Table 8):

R ) RLO + r(RUP - RLO)

(17)

r can take any value from 0 to 1; hence, it can be treated as the fifth uncertain parameter with the uniform distribution function. The uncertain parameters are assumed to be independent. (a) Synthesis at Nominal Conditions. This was carried out by using the modular computer package PROSYN-MINLP. An MINLP model for heat integration and HEN synthesis based on the two-stage superstructure by Yee and Grossmann27 has been used within PROSYN. GAMS/MINOS was used for the solution of nonlinear subproblems and GAMS/OSL for the solution of MILP master problems. The mathematical model includes 102 single equations, 97 single variables, and 12 binary variables (8 of


Figure 7. Optimal HEN structure at the nominal conditions (example 1). Table 9. Quadrature Points and Probabilities of Uncertain Parameters quadrature points for uncertain temperatures θj

j 1 2 3 4 5

- 2.718σ θN - 1.615σ N θ θN + 1.615σ θN + 2.718σ

θN

quadrature points for uncertain film coefficients

vj

rj

vj

0.0070 0.1545 0.6770 0.1545 0.0070

0 0.2 0.5 0.8 1

0.1185 0.2393 0.2844 0.2393 0.1185

them for heat matches between process streams in two stages and 4 binary variables for matches between streams and utilities). The areas of heat-exchanger units are treated as design variables. The optimal HEN structure (Figure 7) is found in the third major iteration (5.1 s of CPU on IBM RISC 6000 Model 43P) with heat matches: H1-C2-1 (hot stream H1-cold stream C2stage1), H2-C1-1, H1-C1-2, and H2-CW. The total annual cost is equal to 129.4 k$/year, 117.4 k$/year from it being investment cost. Note that the solution is a threshold problem without any heating requirements. (b) Synthesis under Uncertainty by Using the RDS Synthesis Strategy. The objective is to determine a HEN structure with the lowest expected annual cost (EC) for the utility and investment of heat exchangers. Set I is defined first. It is composed of 32 (25) vertices, each given by a different combination of extreme points obtained by using expected deviations in Tables 7 and 8. For each substructure predicted by the MILP master problem, the EC(θk) values of uncertain parameters are calculated by using Gaussian integration formula: sequential optimizations of the original problem are performed for the structure at five values (j ) 1, ..., 5) of each uncertain parameter while the remaining four uncertain parameters are kept at the nominal values. The quadrature points and the corresponding probabilities of normally distributed uncertain temperatures (Table 9) are obtained similarly as shown in Table 2. The points represented by rj for the uniformly distributed parameter r are obtained by using eq 14 with a ) 0 and b ) 1. The corresponding coefficients represented by v are obtained by using the Gaussian integration formula:

vj )

b-a Bj fU(θj) 2

fU (θj) )

1 b-a

j ) 1, 2, ..., 5

(18)

where fU is the density function of the uniform distribution. Note that 21 simple sequential optimizations are required for each substructure to prepare the data for

Figure 8. Optimal flexible HEN structure (example 1).

the (P-W) problem. The average CPU time for each sequential optimization is only 0.16 s. The set of all vertices (I) is divided into subsets of critical (CV) and noncritical vertices (I\CV). The former are determined on the assumption that film and fouling coefficients monotonically affect the design variables (areas). The subset CV contains 16 vertices (24) with fixed “worst” values of film and fouling coefficients (r ) 0). The (P-W) problem is then performed to obtain the basic vertices (BV) and weights for the expected cost evaluation in the (P-RDS) NLP subproblem. The average CPU time for the (P-W) problem is 0.45 s. The (PRDS) subproblem is solved simultaneously at vertices from the union CV ∪ BV which usually contains 18 vertices (16 critical and 2 noncritical vertices). The size of the resulting (P-RDS) problem is 1820 equations, 1460 continuous variables, and 12 binary variables. Note that the solution of the (P-RDS) problem is used to derive outer approximations. Finally, the MILP master problem is solved to predict a new topology which is then processed in a new main iteration. The course of the OA/ER algorithm is shown in Table 10. The termination criterion by Viswanathan and Grossmann25 was used and modified in the following way: the optimization was terminated when the results of two successive NLP steps were worse than the best current solution. The solution with the lowest expected cost (Figure 8) is found in the seventh major iteration after 83 min of CPU time on the IBM RISC 6000 Model 43P (plus approximately 27 s for producing the data and solving 7 (P-W) problems). The expected total cost amounts to 185.4 k$/year (129.4 k$/year at nominal conditions). The investment cost of the flexible structure amounts to 160.0 k$/year, which is 36% more than the investment cost of the optimal structure at nominal conditions (117.4 k$/year), which is the price to be paid to attain a feasible operation under uncertain conditions. Note that the structures in iterations 4, 7, and 8 in Table 10 have very close results, which indicates many optimal or near-optimal solutions of the example. (c) Optimization of the Best structures by Using Known Methods. The structures 4, 7, and 8 were also optimized by using the deterministic method with the (P-D) problem and by using the two-stage stochastic model with Monte Carlo simulation. The deterministic optimization was performed simultaneously in 16 critical points with equal weights of 1/16 giving an average cost of flexible design. When the two-stage stochastic method was applied, the design variables were optimized in the first step by the direct search method DCA (Kravanja and Glavic28). In the second stage Monte Carlo simulations were performed sequentially at approximately 300 randomly generated points and the expected cost of flexible design was calculated. The number of points corresponds to a 95% confidence limit

2692 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 Table 10. Solution of Example 1 NLP

MILP master

iteration

structure

EC (k$/year)

CPU (s)

EC (k$/year)

CPU (s)

1 2 3 4 5 6 7 8 9

H1-C1-1, H2-C1-1, H1-C2-2, H1-CW, C2-St H1-C1-1, H1-C2-1, H2-C1-1, H1-C1-2, H1-CW, C2-St H1-C1-1, H1-C2-1, H2-C1-1, H1-C1-2, H1-CW, H2-CW, C2-St H1-C1-2, H1-C2-2, H2-C1-2, H1-CW, C2-St H1-C2-1, H2-C1-1, H1-C1-2, H2-C1-2, H1-CW, C2-St H1-C1-1, H1-C2-1, H2-C1-1, H1-CW, C2-St H1-C1-1, H1-C2-1, H2-C1-1, H1-C1-2, H2-C1-2, H1-CW, C2-St H1-C1-1, H1-C2-2, H2-C1-2, H1-CW, C2-St H1-C1-1, H1-C2-1, H2-C1-1, H1-C1-2, H1-CW, C1-St, C2-St

200.6 190.2 199.7 187.8 202.3 187.8 185.4 186.3 197.3

34 70 99 46 86 28 115 35 75

179.3 178.9 182.3 176.2 181.2 178.9 179.7 180.2

583 849 820 791 714 722 828 741

structure

P-RDS

P-D

P-S (Monte Carlo)

H1-C1-2, H1-C2-2, H2-C1-2, H1-CW, C2-St H1-C1-1, H1-C2-1, H2-C1-1, H1-C1-2, H2-C1-2, H1-CW, C2-St H1-C1-1, H1-C2-2, H2-C1-2, H1-CW, C2-St

187.8 185.4 186.3

195.9 200.9 192.8

187.1 187.8 184.5

Table 11. Best Solutions of the (P-RDS), (P-D), and (P-S) Models (Example 1) EC (k$/year) 4 7 8

Table 12. Optimal Design Variables (Example 1) AH1-C1-1 (m2) AH1-C2-1 (m2) AH2-C1-1 (m2) AH1-C1-2 (m2) AH2-C1-2 (m2) AH1-CW (m2) AC2-St (m2)

P-RDS

P-D

P-S (Monte Carlo)

26.9 262.6 105.3 15.3 168.8 69.7 161.5

26.5 273.7 110.2 6.6 163.7 69.7 159.7

35.1 276.7 107.8 25.2 166.4 69.7 159.0

Table 13. Model Sizes (Example 1) single equations single variables CPU (s)

P-RDS

P-D

P-S (Monte Carlo)

1820 1460 115

1618 1301 53

107 117 715

for the result being within an error of (1 k$/year of the true value. The results of the RDS method and the Monte Carlo stochastic method are rather close (Table 11), while the deterministic results overestimate these results up to 8.4%, which was expected since in the critical vertices only the worst (lower) extreme values for the film and fouling coefficients were considered. Table 12 presents the optimal design variables for structure 7 obtained by different methods. The sizes of the (P-RDS) model and the (P-D) model are of the same order of magnitude (Table 13); hence, the computational effort to solve the (P-RDS) model is not much larger than that in the case of the (P-D) model. The stochastic model is much smaller than the other two, but repetitive Monte Carlo simulations are computationally exhaustive. (d) Testing the Feasibility of the Solutions. Finally, the feasibility of RDS, deterministic, and Monte Carlo solutions was tested. Five points were taken for each uncertain parameter, as shown in Table 9, and combined into 3125 (55) quadrature points. The simultaneous stochastic optimization was not possible because the model was too large (approximately 304 000 equations and 253 500 variables). Hence, the feasibility of the solutions was tested by performing the stochastic optimization sequentially in 3125 quadrature points, while the design variables were fixed to the optimal values obtained by all three methods. It was established that the operating conditions can be selected in all quadrature points to achieve a feasible operation, which confirms the feasibility of all three designs. Besides, the optimal values of the objective function in all quadrature points were weighted according to the joint probability

Table 14. Results of (P-RDS) and (P-D) Models (Example 1) EC (k$/year)

EC*(k$/year)

P-RDS

P-D

185.4

200.9

fix dP-RDS

dfix P-D

185.6

187.8

function, giving the expected stochastic cost (EC*) for fix ) and the (P-D) solution the (P-RDS) solution (dP-RDS fix (dP-D) as shown in Table 14 for structure 7. The total CPU time needed to perform 3125 sequential optimizations was approximately 575 s. It is interesting to note that the relative error of the (P-RDS) model is only 0.1% while the error of the (P-D) model is 7.0%. Example 2: MINLP Synthesis of a Flexible HeatIntegrated Distillation Sequence and Its HEN. The synthesis of the heat-integrated distillation sequence for a four-component mixture at nominal conditions was extensively studied in previous work.29 In this work the synthesis of a slightly modified problem was first performed at nominal conditions which was followed by synthesis under uncertainty. The compact superstructure29 was used to present the alternative distillation sequences. The superstructure was modeled as an MINLP problem using GillilandUnderwood-Fenske’s model for distillation columns and the one-stage model for heat integration and HEN synthesis by Yee and Grossmann.25 The distillation superstructure and the HEN superstructure were merged together as shown in Figure 9. The condensing streams, Con1, Con2, and Con3 of columns A/BCD, AB/CD, and ABC/D, respectively, were defined as hot streams of the HEN. The boiling streams, Reb1, Reb2, and Reb3 of columns A/BCD, AB/CD, and ABC/D, respectively, were defined as cold streams, and the cooling water (CW) and the steam (St) as utilities. The design variables were the number of trays and the diameters of the distillation columns and the areas of reboilers and condensers. The objective was to maximize the annual profit (P) composed of the product incomes minus the raw materials cost, utility cost, and the annualized investment cost. The data of the problem are shown in Table 15. (a) Synthesis at Nominal Conditions. This was carried out first. The mathematical model included 497 equations, 542 continuous variables, and 25 binary variables (13 for separation paths and 12 for heat


Figure 9. Superstructure for example 2. Table 15. Data for Example 2a component

mole fraction

A. benzene B. toluene C. o-xylene D. diphenyl

0.386 0.292 0.220 0.102

utilities

T (K)

cost

cooling water (CW) steam (St) electricity

293 601

0.159 $/(106 kJ) 5.333 $/(106 kJ) 0.03 $/(kW‚h)

a Feed flow rate ) 0.5484 kmol/s; feed temperature ) 374 K; desired products: pure benzene, pure toluene, pure xylene, and pure diphenyl; recovery in all columns is specified to be 99.5%; capital charge factor is specified to be 0.33 year-1.

matches). The problem was solved by using the PROSYNMINLP synthesizer. The optimal structure at nominal conditions (Figure 10) had the separation sequence A/BCD, AB/CD, and ABC/D. The annual profit amounted to 10.65 M$/year. The investment costs of distillation columns and HEN units were 2.34 M$/year. Optimal values of design variables are shown in Table 16. The heat-integration scheme can be clearly seen on the T versus φ diagram (Figure 11): the reboiler of column A/BCD is completely covered by the condensers of columns AB/CD and ABC/D. The reboilers of later columns are heated by steam while cooling water is needed only in the condenser of the first column. (b) Synthesis under Uncertainty by Using the RDS Synthesis Strategy. The problem involves six

Figure 10. Optimal solution at the nominal conditions (example 2).

uncertain parameters (feed flow rate, feed temperature, cooling water temperature, and compositions of components A, B, and C in the feed) with probability distributions shown in Table 17. The set of extreme points (I) includes 64 (26) vertices. A simple parameter analysis was performed for each substructure predicted by the MILP master problem to analyze gradients ∂d/∂θ ≈ ∆d/∆θ and to find out if some of the uncertain parameters affect the design variables monotonically. Such parameters can be fixed to their critical values, either at the upper or at the lower bound, reducing the number of critical vertices (CV). The remaining vertices compose the subset of noncritical vertices (I\CV).

2694 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 Table 16. Design Variables of the Optimal Structure at Nominal Conditions (Example 2) column design variables

A/BCD

number of trays diameter (m) reboiler area (m2) condenser area (m2) area of heat matches (m2)

50 5.7 781

Figure 11. T vs φ diagram for optimal solution at the nominal conditions (example 2). Table 17. Data for Uncertain Parameters (Example 2) uncertain parameter

distribution function

nominal value

expected deviations

FF (kmol/s) TF (K) TCW (K) xA xB xC

N(0.5484,0.05484) N(374,6.67) N(293,3.33) N(0.386,0.0103) N(0.292,0.0078) N(0.220,0.0059)

0.5484 374 293 0.386 0.292 0.220

(0.1491 (18 (9 (0.0280 (0.0212 (0.0159

The expected value of the annual profit EP(θk) is determined for each uncertain parameter through repetitive sequential optimizations at five points (θN 2.7σ, θN - 1.6σ, θN, θN + 1.6σ, θN + 2.7σ) while the remaining five uncertain parameters are kept at the nominal values. The average CPU time for sequential optimization is 3 s. Twenty-five (6 × 4 + 1) sequential optimizations altogether have to be performed for each structure to provide data for the (P-W) problem which determines the basic vertices (BV) that are included in

AB/CD

ABC/D

68 4.0 218

24 2.9 98

condenser (AB/CD) - reboiler (A/BCD): 1279 condenser (ABC/D) - reboiler (A/BCD): 609

the (P-RDS) subproblem. The average CPU time for the (P-W) problem is 30 s. The (P-RDS) subproblem for each particular structure is solved simultaneously at the points from the union BV ∪ CV. The MILP master problem predicts a new topology for the next main iteration. The iterations of the OA/ER algorithm are shown in Table 18. The optimal structure (Figure 12) is found in the third major MINLP iteration after 303 min of CPU time on the IBM RISC 6000 Model 43P (plus approximately 225 s for the production of the data and 90 s for the execution of three (P-W) problems). The optimal solution has the expected annual profit of 10.34 M$/year. The optimal separation sequence of the flexible structure is the same as that of the nonflexible one, with the exception that the HEN structure has been changed: the condensing stream of the column AB/CD needs additional cooling; hence, a water cooler is added to the HEN structure (cooler 2 in Figure 12). The situation can clearly be seen in the T versus φ diagram (Figure 13) for the following combination of uncertain parameters: FF ) 0.6975 kmol/s, TF ) 392 K, TCW ) 302 K, xA ) 0.414, xB ) 0.313, xC ) 0.236, and xD ) 0.037. The investment cost of the flexible structure is 2.80 M$/year, that is, 0.46 M$/year (20%) higher when compared to the investment cost of the nonflexible structure. The optimal values of the design variables are shown in Table 19. As expected, the design variables of the flexible solution overestimate the design variables determined at the nominal conditions shown in Table 16. It is interesting to mention that the uncertain parameters feed flow rate, feed temperature, and cooling water temperature monotonically affect the design variables in structures 1 and 2. These two structures have only 8 (23) critical vertices and three additional noncritical vertices determined by (P-W) problems. The size of the first two (P-RDS) problems amounts to 5400 equations and 5600 continuous variables. The third and fourth structure have only two monotonic parameters (the feed flow rate and the cooling water temperature); hence, the number of critical vertices is 16 (24). Three (for structure 3) and two (for structure 4) noncritical vertices are determined by the (P-W) problems. The size

Table 18. Solution of Example 2 iteration

structure

1

separation sequence: A/BCD, AB/CD, ABC/D heat-integration scheme: all condensers-CW, all reboilers-St separation sequence: A/BCD, AB/CD, ABC/D heat-integration scheme: Con3-Reb2, Con2-Reb1, Con1-CW, Reb3-St, and Reb2-St separation sequence: A/BCD, AB/CD, ABC/D heat-integration scheme: Con3-Reb1, Con2-Reb1, Con1-CW, Con2-CW, Reb2-St, and Reb3-St separation sequence: A/BCD, AB/CD, ABC/D heat-integration scheme: Con3-Reb2, Con2-Reb1, Con3-Reb1, all condensers-CW, Reb2-St, and Reb3-St

2 3 4

EP (M$/year)

CPU (min)

NLP MILP NLP MILP

6.83 15.75 7.25 14.66

11.5 25.3 18.5 12.5

NLP MILP

10.34 15.80

235.0 77.2

NLP MILP

7.35

180


Figure 12. Optimal flexible structure (example 2).

tions of the uncertain parameters. The objective was to maximize the sum of the total annual profits in the vertices that are weighted with the same value of 1/16. The deterministic result is 12.90 M$/year which exceeds the (P-RDS) result by 25%. The main reason for the deviation is that the feed flow rate was fixed at the critical upper bound which overestimates the product income. The optimal values of the design variables are overestimated as well (see Table 20). The sizes of the (P-RDS) and the (P-D) model are rather similar (Table 21); hence, the computational effort to solve the (P-RDS) model is not much greater than that to solve the (P-D) model. However, the (P-RDS) model gives a more reliable solution. (d) Testing the Feasibility of the (P-RDS) and (P-D) Solutions. According to the suggested strategy in Figure 6, the feasibility of both solutions was tested at the end. Sequential optimizations were performed in 15 625 quadrature points (56) which were defined by using five values (ND ) 5) for each uncertain parameter (NP ) 6) with respect to the optimal values of the design variables that were fixed to the values determined by the use of the (P-RDS) model (Table 19) and the (P-D) model (Table 20). The feasibility of both designs was confirmed. The expected profits were also determined by using the joint probabilities of each quadrature point, giving the values of the stochastic expected profits (EP*) as shown in Table 22. As can be seen, the proposed reduced dimensional stochastic approach has a 4.6% error while the error of the deterministic result amounts to 31.6%, which confirms the fact that the (P-RDS) problem yields significantly more reliable solutions than the (P-D) problem, especially when feasibility constraints and the objective

Figure 13. T vs φ diagram for the optimal flexible solution (example 2).

of the third (P-RDS) problem (9320 equations and 9690 variables) entails a much larger CPU time for its optimization. (c) Optimization of the Optimal Structure by Using Known Methods. As mentioned before, it was not possible to solve this example by using two-stage Monte Carlo simulation with the existing hardware. Only the deterministic optimization of the optimal structure was performed simultaneously at 16 critical vertices without consideration of the distribution funcTable 19. Design Variables of Flexible Structure (Example 2)

column design variables number of trays diameter (m) reboiler area (m2) condenser area (m2) area of heat matches (m2)

A/BCD 60 6.3 1121

AB/CD

ABC/D

55 4.5 486 41 condenser (AB/CD) - reboiler (A/BCD): 1448 condenser (ABC/D) - reboiler (A/BCD): 602

28 3.2 105

Table 20. Design Variables of the Flexible Structure Using the Deterministic Method (Example 2) column design variables number of trays diameter (m) reboiler area (m2) condenser area (m2) area of heat matches (m2)

A/BCD 61 6.3 1174

AB/CD

ABC/D

57 4.6 466

28 3.3 124

condenser (AB/CD) - reboiler (A/BCD): 1717 condenser (ABC/D) - reboiler (A/BCD): 631

2696 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 Table 21. Sizes of (P-RDS) and (P-D) Models (Example 2) single equations single variables CPU (min)

P-RDS

P-D

9320 9690 236

7850 8160 262

Table 22. Results of (P-RDS) and (P-D) Models (Example 2) EP (M$/year)

EP* (M$/year)

P-RDS

P-D

10.34

12.90

fix dP-RDS

dfix P-D

9.89

9.80

function are posed only for the critical vertices. We could expect an even larger error of the (P-D) problem if the economic uncertain parameters were considered in the optimization. Conclusions This paper has presented the evolution of the robust strategy for the synthesis of the complex engineering problems under uncertainty. The strategy is based on the new method called the reduced dimensional stochastic optimization (RDS method). The main idea of this method is to treat complex optimization problems in the simultaneous procedure on a relatively small set of vertex points. The expected value of the economic objective function is approximated through the simplified linear expression using only few extreme points instead of many quadrature points, while the feasibility tests of the structure are performed simultaneously in the critical points. It was assumed in this work that the critical points can only be the combination of the extreme values of the uncertain parameters. This simplification was suitable in the sense of building the strategy for solving process synthesis problems under uncertainty which are very extensive and complicated for optimization. A special MILP (P-W) problem has been proposed to determine the vertices and weights for the evaluation of the expected value. The method was illustrated by a small example of a heat-exchanger network optimization which was also solved by using the stochastic method. On the basis of this new method, a synthesis strategy (RDS strategy) was proposed which offers the designers a way for handling the complex synthesis problems under uncertainty. The RDS strategy is based on the modified OA/ER algorithm for the MINLP problems optimization. The main NLP (P-RDS) problem is performed for the particular structure simultaneously only in few points, thus keeping the size of the mathematical model reasonable. The procedure has been implemented successfully in the PROSYN-MINLP synthesizer. Two example problems of the process synthesis under uncertainty have been presented to illustrate the potential of the proposed RDS strategy. The examples show the attractive features of our method; the results are oriented toward the stochastic results. However, the improvement of the accuracy of our method and the simplification of the setup procedure still remain as the challenging tasks for future work. Current experiences indicate that greater accuracy is connected with more extensive and more complicated setup procedures.

Acknowledgment The authors are grateful to the Slovenian Ministry of Science and Technology for partial financial support of this work. Nomenclature Indices cv ) critical vertices i ) vertices j ) points for the calculation of the expected value k ) uncertain parameters l ) heat exchangers m ) quadrature points Variables, Functions A ) area of heat exchanger, m2 C ) economic objective function (cost), k$/year d ) vector of design variables EC ) expected cost, k$/year EC* ) expected cost calculated with the stochastic method at fixed design variables, k$/year EP ) expected profit, M$/year EP* ) expected profit calculated with the stochastic method at fixed design variables, M$/year f ) vector of design functions FF ) feed flow rate, kmol/s fG ) density function of normal distribution fU ) density function of uniform distribution g ) vector of model inequalities h ) vector of model equations P ) economic objective function (profit), M$/year T ) temperature, K x ) vector of operating variables y ) vector of binary variables Parameters B ) weight corresponding to the zero of the Legendre polynomial CF ) heat capacity flow rate, kW/K N ) number of the uncertain parameters realizations for discretization of an infinite problem ND ) number of selected points for the expected value estimation NM ) number of points with monotonic influence on design variables NP ) number of uncertain parameters p ) probability TCW ) cooling water temperature, K TF ) feed temperature, K TIN and TOUT ) supply and target temperatures, respectively, K t ) zero of the Legendre polynomial U ) heat-transfer coefficient, W/(m2K) v ) coefficient in the Gaussian integration formula w ) weight of the vertex xA, xB, xC ) mole fractions of components A, B, C, respectively Abbreviations C1, C2 ) cold process streams Con1, Con2, Con3 ) condensers of distillation columns CW ) cooling water dfix) sequential stochastic optimization with fixed design variables H1, H2 ) hot process streams M1, ..., M5 ) single-choice mixers N ) normal distribution function P1, P2, P3 ) pumps P-D ) deterministic mathematical model

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2697 P-S ) stochastic mathematical model P-RDS ) mathematical model for the reduced dimensional stochastic optimization P-W ) MILP model for the determination of basic vertices and their weights Reb1, Reb2, Reb3 ) reboilers S1, ..., S5 ) single-choice splitters St ) steam RDS ) Method (strategy) for the reduced dimensional stochastic optimization V1, V2, V3 ) valves Greek Symbols R ) individual film and fouling coefficient, W/(m2 K) ∆lnT ) logarithmic temperature difference, K φ ) heat flow rate, kW θ ) vector of uncertain parameters θN, θLO, θUP ) nominal value and lower and upper bounds of the uncertain parameter θ σ ) standard deviation

k ) 1, n ) 1 k ) 2, n ) 1, 2 θ2,i )

{

The deviation of EP(θk) from P at the nominal conditions for k ) 1,2 can be expressed as the linear combination of deviations of four extreme vertices from the nominal point:

θLO 2 , if 0 < i e 2 or 4 < i e 6 θUP 2 , if 2 < i e 4 or 6 < i e 8

k ) 3, n ) 1, 2, 3, 4

θLO 3 , if 0 < i e 1 or 2 < i e 3 or 4 < i e 5 or 6 < i e 7 θ3,i ) UP θ3 , if 1 < i e 2 or 3 < i e 4 or 5 < i e 6 or 7 < i e 8 The system for the calculation of weights wi is then N N LO N N EP(θ1) ) w1P(θLO 1 ,θ2 ,θ3 ) + w2P(θ1 ,θ2 ,θ3 ) + N N LO N N w3P(θLO 1 ,θ2 ,θ3 ) + w4P(θ1 ,θ2 ,θ3 ) +

N N UP N N w7P(θUP 1 ,θ2 ,θ3 ) + w8P(θ1 ,θ2 ,θ3 ) LO N N LO N EP(θ2) ) w1P(θN 1 ,θ2 ,θ3 ) + w2P(θ1 ,θ2 ,θ3 ) + UP N N UP N w3P(θN 1 ,θ2 ,θ3 ) + w4P(θ1 ,θ2 ,θ3 ) + LO N N LO N w5P(θN 1 ,θ2 ,θ3 ) + w6P(θ1 ,θ2 ,θ3 ) +

N LO N N N EP(θ1) - P(θN 1 ,θ2 ) ) w1[P(θ1 ,θ2 ) - P(θ1 ,θ2 )] +

UP N N UP N w7P(θN 1 ,θ2 ,θ3 ) + w8P(θ1 ,θ2 ,θ3 )

N N N UP N w2[P(θLO 1 ,θ2 ) - P(θ1 ,θ2 )] + w3[P(θ1 ,θ2 ) N UP N N N P(θN 1 ,θ2 )] + w4[P(θ1 ,θ2 ) - P(θ1 ,θ2 )] (A1) N N LO N N EP(θ2) - P(θN 1 ,θ2 ) ) w1[P(θ1 ,θ2 ) - P(θ1 ,θ2 )] +

N LO N N UP EP(θ3) ) w1P(θN 1 ,θ2 ,θ3 ) + w2P(θ1 ,θ2 ,θ3 ) + N LO N N UP w3P(θN 1 ,θ2 ,θ3 ) + w4P(θ1 ,θ2 ,θ3 ) + N LO N N UP w5P(θN 1 ,θ2 ,θ3 ) + w6P(θ1 ,θ2 ,θ3 ) +

UP N N N LO w2[P(θN 1 ,θ2 ) - P(θ1 ,θ2 )] + w3[P(θ1 ,θ2 ) -

N LO N N UP w7P(θN 1 ,θ2 ,θ3 ) + w8P(θ1 ,θ2 ,θ3 )

N N UP N N P(θN 1 ,θ2 )] + w4[P(θ1 ,θ2 ) - P(θ1 ,θ2 )] (A2)

8

wi ) 1 ∑ i)1

Since the sum of wi (i ) 1, ..., 4) is equal to 1, (A1) and (A2) become N LO N LO N EP(θ1) - P(θN 1 ,θ2 ) ) w1P(θ1 ,θ2 ) + w2P(θ1 ,θ2 ) + N UP N N N w3P(θUP 1 ,θ2 ) + w4P(θ1 ,θ2 ) - P(θ1 ,θ2 ) (A3) N N LO N UP EP(θ2) - P(θN 1 ,θ2 ) ) w1P(θ1 ,θ2 ) + w2P(θ1 ,θ2 ) +

+

UP w4P(θN 1 ,θ2 )

-

N P(θN 1 ,θ2 )

] [ ] [ ] [

(A4)

+

N P(θUP 1 ,θ2 ) w4 UP P(θN 1 ,θ2 )

]

4

wi ) 1 ∑ i)1

Once the values of wi, i ) 1, ..., 8 are known, the expected value EP can be approximated in the simultaneous optimization by using eight combined vertices:

UP LO LO UP UP w3P(θLO 1 ,θ2 ,θ3 ) + w4P(θ1 ,θ2 ,θ3 ) + LO LO UP LO UP w5P(θUP 1 ,θ2 ,θ3 ) + w6P(θ1 ,θ2 ,θ3 ) + UP LO UP UP UP w7P(θUP 1 ,θ2 ,θ3 ) + w8P(θ1 ,θ2 ,θ3 ) (B2)

Literature Cited

N N EP(θ1) P(θLO P(θLO 1 ,θ2 ) 1 ,θ2 ) ) w1 + w 2 LO UP + EP(θ2) P(θN P(θN 1 ,θ2 ) 1 ,θ2 ) N P(θUP 1 ,θ2 ) w3 LO P(θN 1 ,θ2 )

(B1)

LO LO LO LO UP EP ) w1P(θLO 1 ,θ2 ,θ3 ) + w2P(θ1 ,θ2 ,θ3 ) +

N Thus, P(θN 1 , θ2 ) can be deleted, resulting in the following system of equations:

[ ] [

{

{

θLO 1 , if 0 < i e 4 θUP 1 , if 4 < i e 8

N N UP N N w5P(θUP 1 ,θ2 ,θ3 ) + w6P(θ1 ,θ2 ,θ3 ) +

Appendix A

LO w3P(θN 1 ,θ2 )

θ1,i )

(A5)

Appendix B In a three-dimensional example (NP ) 3 and k ) 1, 2, 3), we obtain the following values of the uncertain parameters at the vertices i ) 1, 2, ..., 8:

(1) Floudas, C. A.; Grossmann, I. E. Synthesis of Flexible Heat Exchanger Networks with Uncertain Flowrates and Temperatures. Comput. Chem. Eng. 1987, 11, 319. (2) Paules, G. E.; Floudas, C. A. Stochastic Programming in Process Synthesis: A Two-Stage Model with MINLP recourse for Multiperiod Heat-Integrated Distillation Sequences. Comput. Chem. Eng. 1992, 16, 189. (3) Papalexandri, K. P.; Pistikopoulos, E. N. A Multiperiod MINLP Model for the Synthesis of Flexible Heat and Mass Exchange Networks. Comput. Chem. Eng. 1994, 18, 1125. (4) Pistikopoulos, E. N.; Ierapetritou, M. G. Novel Approach for Optimal Process Design under Uncertainty. Comput. Chem. Eng. 1995, 19, 1089. (5) Acevedo, J.; Pistikopoulos, E. N. A Parametric MINLP Algorithm for Process Synthesis Problems under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 147.

2698 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 (6) Acevedo, J.; Pistikopoulos, E. N. A Multiparametric Programming Approach for Linear Process Engineering Problems under Uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717. (7) 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. (8) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Robust Stability Considerations in Optimal Design of Dynamic Systems under Uncertainty. J. Process. Control. 1997, 7 (5), 371. (9) Dua, V.; Pistikopoulos, E. N. Optimization Techniques for Process Synthesis and Material Design under Uncertainty. Trans. Inst. Chem. Eng. 1998, 76 (Part A), 408. (10) Grossmann, I. E.; Halemane, K. P.; Swaney, R. E. Optimization Strategies for Flexible Chemical Processes. Comput. Chem. Eng. 1983, 7, 439. (11) Grossmann, I. E.; Straub, D. A. Recent Developments in the Evaluation and Optimization of Flexible Chemical Processes. In Computer Oriented Process Engineering; Puigjaner, L., Espun˜a, A., Eds.; Elsevier: Amsterdam, The Netherlands, 1991; p 49. (12) Petracci, N. C.; Hoch, P. M.; Eliceche, A. M. Flexibility Analysis of an Ethylene Plant. Comput. Chem. Eng. Suppl. 1996, 20, S443. (13) Hoch, P. M.; Eliceche, A. M. Flexibility Analysis Leads to a Sizing Strategy in Distillation Columns. Comput. Chem. Eng. Suppl. 1996, 20, S139. (14) Straub, D. A.; Grossmann, I. E. Integrated Stochastic Metric of Flexibility for Systems with Discrete State and Continuous Parameter Uncertainties. Comput. Chem. Eng. 1990, 14, 967. (15) Ierapetritou, M. G.; Acevedo, J.; Pistikopoulos, E. N. An Optimization Approach for Process Engineering Problems under Uncertainty. Comput. Chem. Eng. 1996, 20, 703. (16) Ierapetritou, M. G.; Pistikopoulos, E. N. Batch Plant Design and Operations under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 772. (17) Liu, M. L.; Sahinidis, N. V. Optimization in Process Planning under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 4154. (18) Epperly, T. G.; Ierapetritou, M. G.; Pistikopoulos, E. N. On the Global and Efficient Solution of Stochastic Batch Plant Design Problems. Comput. Chem. Eng. 1997, 21, 1411. (19) Ostrovsky, G. M.; Volin, Yu. M.; Golovaskhin, D. V. Optimization Problem of Complex System under Uncertainty. Comput. Chem. Eng. 1998, 22, 1007.

(20) Swaney, R.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part I: Formulation and Theory. AIChE J. 1985, 31, 621. (21) Halemane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AIChE J. 1983, 29, 425. (22) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods for Chemical Process Design; Prentice Hall: Englewood Cliffs, NJ, 1997; Chapter 21. (23) Duran, M. A.; Grossmann, I. W. An Outer-Approximation Algorithm for a Class of Mixed-Integer Nonlinear Programs. Math. Prog. 1986, 36, 307. (24) Kocis, G. R.; Grossmann, I. E. Relaxation Strategy for the Structural Optimization of Process Flow Sheets. Ind. Eng. Chem. Res. 1987, 26, 1869. (25) Viswanathan, J.; Grossmann, I. E. A Combined Penalty Function and Outer-Approximation Method for MINLP Optimization. Comput. Chem. Eng. 1990, 14, 769. (26) Kravanja, Z.; Grossmann, I. E. New Developments and Capabilities in PROSYNsAn Automated Topology and Parameter Process Synthesizer. Comput. Chem. Eng. 1994, 18, 1097. (27) Yee, T. F.; Grossmann, I. E. Simultaneous Optimization Models for Heat IntegrationsII. Heat Exchanger Network Synthesis. Comput. Chem. Eng. 1990, 14, 1165. (28) Kravanja, Z.; Glavic, P. Cost Targeting for HEN through Simultaneous Optimization Approach: A Unified Pinch Technology and Mathematical Programming Design of Large HEN. Comput. Chem. Eng. 1997, 21, 833. (29) Novak, Z.; Kravanja, Z.; Grossmann, I. E. Simultaneous Synthesis of Distillation Sequences in Overall Process Schemes Using an Improved MINLP Approach. Comput. Chem. Eng. 1996, 20, 1425.

Received for review October 2, 1998 Revised manuscript received March 18, 1999 Accepted March 19, 1999 IE980629Z

Mixed-Integer Nonlinear Programming Problem Process Synthesis

Recommend Documents