2368
Ind. Eng. Chem. Res. 2000, 39, 2368-2377
A New Algorithm for Computing Process Flexibility Gennadi M. Ostrovsky, Luke E. K. Achenie,* and Yiping Wang Department of Chemical Engineering, University of Connecticut, U-222, Storrs, Connecticut 06269
Y. M. Volin Karpov Institute of Physical Chemistry, Moscow, Russia
In the design of a chemical process (CP), certain design specifications (for example, those related to process economics, process performance, safety, and the environment) must be satisfied. During the operation of the plant, because design models have uncertainties associated with them, we need to ensure that, within the region of uncertainty, all design specifications are satisfied. In recent years, research has focused on the investigation of the process flexibility (Biegler et al. Systematic Methods of Chemical Process Design; Prentice Hall: Englewood Cliffs, NJ, 1997) based on the feasibility function (Halemane and Grossmann AIChE J. 1983, 29 (3), 425) which is a measure of the CP’s ability to meet design specifications under uncertainty. Several researchers have proposed methods for calculating the process feasibility function, which involves solving a very complex multiextremal and nondifferentiable optimization problem. Current methods for calculation of the flexibility function use an enumeration procedure (explicit or implicit), which in the worst case can require a large number of iterations. To try to address this issue, in this paper, we have introduced an efficient approach, which avoids enumeration. Through examples, we have shown that the new method leads to a small number of iterations and has low CPU requirements. Introduction In the design of a chemical process (CP), certain specifications (for example, those related to process economics, process performance, process safety, and the environment) must be satisfied. These specifications can be modeled in terms of design and control variables and parameters. Nevertheless, the mathematical models always have uncertainties (parametric and structural) associated with them. In this paper, we will assume that structural uncertainties are not present and instead only focus on parametric uncertainties. Examples of these are kinetic constants, transfer coefficients, parameters related to thermophysical property models, a change in catalyst activity, fluctuations in the stream flow rate, temperature, pressure, and stream compositions. At the operational stage, we need to ensure that, within the region of uncertainty, all design specifications are satisfied. In recent years, research has focused on the investigation of the process flexibility1 via the feasibility function,2 which is a measure of the CP’s ability to satisfy design specifications under uncertainty. Analysis of the flexibility of the CP is especially important in assessing process safety and product quality. The efficacy of any process feasibility function evaluation strategy should be measured by (a) the reliability of the method, (b) the range of process constraints it can handle, and (c) the speed with which the method can calculate the feasibility function. Several researchers have proposed methods for calculating the flexibility function with the above attributes in mind. The calculation of process flexibility involves the solution of a very * To whom correspondence should be addressed. E-mail:
[email protected]. Phone: (860) 486 2756. Fax: (860) 486 2959.
Figure 1. Flowsheet of a chemical process.
complex multiextremal and nondifferentiable optimization problem. Current methods use an enumeration procedure, which in the worst case can require a large number of iterations. To try to address this issue, we have introduced an efficient approach, which avoids enumeration. This paper is organized as follows. In section 1, we will present a motivational example. In section 2, we will outline the mathematical basis for our approach. This is followed by the development of the new approach in section 3. The efficiency of the approach is demonstrated through the solution of three chemical process problems in section 4. Finally, in section 5, we will discuss the strengths of the approach and ongoing efforts to improve the algorithm. 1. Motivational Example Let us consider the optimization of a reactor and heat exchanger system as shown in Figure 1.2 The reaction is assumed to be first order and exothermal of the type, A f B. The flow rate through the heat exchanger is adjusted to maintain the reactor temperature below
10.1021/ie9905207 CCC: $19.00 © 2000 American Chemical Society Published on Web 05/10/2000
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2369
T1max and to get a minimum of 90% conversion. The material and energy balance equations for the reactor can be represented as
F0(CA0 - CAl)/CA0 ) VkR exp(-E/RT1)CA1 (-H)F0(CA0 - CA1)/CA0 ) F0CP(T1 - T0) + QHE For the heat exchanger, the energy balances are
QHE ) F1CP(T1 - T2) ) Cpw(Tw2 - Tw1)W QHE ) AU
(T1 - Tw2) - (T2 - Tw1) ln{(T1 - Tw2)/(T2 - Tw1)}
where F0, T0, and CA0 are the feed flow rate (kmol h-1), temperature of the feed (K), and the concentration of the reactant in the feed (kmol m-1), respectively, V, T1, and CA1 are the values of the reactor volume (m3), reaction temperature (K), and the concentration of the reactant A in the product (kmol m-3), H is the heat of reaction (kJ kg-1 mol-1), F1 is the flow rate of the recycle (k mol h-1), T2 is the recycle temperature, Cp and Cpw (kJ kg mol-1 ) are the heat capacity of the recycle mixture and the cooling water, respectively, Tw1, Tw2, and W are the inlet and outlet temperatures and the flow rate (kg h-1) of the cooling water, respectively, A is the heat-transfer area of the heat exchanger (m2); and U is the overall heat-transfer coefficient (kJ m-2 h-1 k-1). In this problem, V and A are the design variables and T1 (311 e T1 e 389) and Tw2 (301 e Tw2 e 355) are the control variables. The state variables are [CA1, T2, F, W] and the uncertain parameters are [F0, T0, Tw1, kR, U] with nominal values [45.36, 333, 300, 9.8, 1635]. After eliminating the state (i.e., dependent) variables, we have
conv )
T2 )
VkRcA0 exp(-E/RT1) F0 + VkRcA0 exp(-E/RT1)
2(-H)F0 conv 2F0cp(T1 - T0) AU AU (T1 - Tw2) + Twl
where “conv” represents the conversion of reactant A. Using the above notation, we can express the process constraints as follows:
-conv + 0.9 e 0 conv - 1.0 e 0
Tw1 ((3%). In addition, assume the rate constant k0 and the heat-transfer coefficients have been estimated to within (10% of the nominal (or average) values. At the design stage, typically, the nominal values for the uncertain parameters are used to solve the optimization problem for the optimal values V* and A*. Because the actual unknown values of the uncertain parameters are likely to be different from the nominal values, the V* and A* thus obtained are not truly optimal and can lead to wrong inferences on the design. For example, the nominal design can lead to violation of process constraints at the operational stage. Therefore, it is imperative that we account for parameter uncertainties through estimation of the flexibility of the CP. More often than not, the probability distribution function for the uncertain parameters is not known; instead, only the intervals that bound these parameters are known. It is common to assume a uniform probability distribution within such intervals and use stochastic optimization for construction of flexible processes. However, stochastic optimization methods are very computationally intensive. Therefore, it is important to develop deterministic methods for flexible processes. We have developed an efficient technique for estimating the flexibility of the chemical process. In this paper, we will focus our interests on fixed designs. 2. Mathematical Basis Halemane and Grossmann formulated the necessary and sufficient condition of the chemical process (CP) flexibility in the form2
F1(d) e 0
(1)
where the feasibility function, F1(d), has been defined as
F1(d) ) max min max gj(d,z,t), t∈Τ
z∈Z
j∈J
J ) {1, 2, ..., m} (2)
Here, d is a vector of design variables, z is a vector of control variables, Z is a region of admissible values of control variables, and T ) {t: ht e t e t)} is the domain for the uncertain parameter t. If eq 1 is satisfied, then one can find values for the control variables that will ensure the satisfaction of the design specifications. In other words, the process is flexible. The reduced design constraints
gj(d,z,t) e 0,
j ) 1,.., m
-(T2 - Twl) + 11.1 e 0
are obtained from the original mathematical models,
T2 - 389 e 0
φ(d,x,z,t) ) 0
(3)
-T2 + 311 e 0
gj (d,x,z,t) e 0
(4)
T2 - T1 e 0 Assuming we want to design a CP as shown in Figure 1, we need to find optimal sizes of the reactor and heat exchanger. However, we need to consider uncertainty. It is known that, during the operational stage, the input flow rate, F0, the initial feed temperature, T0, and the temperature of the cooling water, Tw1, can vary within certain limits, for example, F0 ((10%), T0 ((2%), and
where x is a vector of state variables. Equation 4 represents design specifications. Equation 3 represents material and energy balances and has the same dimension as x. From (3), assuming one can determine x as an explicit function of d, z, and t (i.e., x ) x(d,z,t)), then (4) becomes
gj(d,z,t) ) gj (d,x(d,z,t),z,t) e 0
2370
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
Let us introduce the function h(d,t):
h(d,t) ) min max gj(d,z,t) z∈Z
j∈J
define the function F1i such that
(5)
When problem (5) is being solved, d and t are parameters. Thus, h is a function of d and t. Now, F1(d) can be represented as
F1(d) ) max h(d,t) t∈T
(6)
Therefore, the evaluation of F1(d) is reduced to the maximization of h(d,t) with respect to t. In general, h(d,t) is multiextremal and nondifferentiable as pointed out by Halemane and Grossmann,2 and therefore very difficult to solve. However, several researchers have made progress in this area. Next, we will discuss some strategies for solving the flexibility function. When gj is jointly convex in z and t, the global solution of the problem will lie at one of the vertices of the parameter set T.2 On the basis of this, we can determine F1(d) by evaluating h(d,t) at every vertex and choose the largest value of h(d,t).2 In the latter approach, the computational effort is, in general, proportional to the number of vertices, 2r (r is the number of uncertain parameters). Swaney and Grossmann introduced the flexibility index, which characterizes the largest uncertainty region that the design can handle for feasible operation.3 With regard to the calculation of the flexibility index, they proposed two algorithms: a heuristic vertex search method and an implicit enumeration scheme that avoids exhaustive enumeration of all the vertices. Nevertheless, these algorithms rely on the assumption that the global solution lies at one of the vertices. To overcome this problem, Grossmann and Floudas proposed an active set strategy,4 which uses the fact that under certain conditions the number of active constraints in (5) is equal to q + 1 (q is the number of control variables). It reduces eq 5 to the identification of all possibly active constraint sets and to the solving of a nonlinear programming problem for each active set. If the function gj(d,z,t) is quasi-concave in z and t and strictly quasi-convex in z for fixed t, then the active constraint strategy gives a global solution of eq 2.4 There is, however, a large computational overhead when the number of active sets is large. In addition, there are hard requirements on the convexity (concavity) of the function gj(d,z,t). Indeed, with respect to the variable z, gj(d,z,t) must be simultaneously quasi-concave and strictly quasi-convex. Kabatek and Swaney5 suggested a modification based on Swaney and Grossmann3 which can find nonvertex solution using heuristics. Ostrovsky and Volin et al. developed the ULB (upper and lower bound) method that can obtain nonvertex solutions for the flexibility function.6 Moreover, it will lead to a global solution if the following conditions are satisfied: (1) The function gj(d,z,t) is quasi-convex in z. (2) The function gj(d,z,t) is quasi-concave in t. The ULB method is based on the branch and bound strategy. It seeks the maximum of h(d,t) in t by partitioning the region T into subregions Ti (in the rest of the paper, we will sometimes use the subscript to designate the subregion). In developing the ULB method, the main issue is to calculate the upper bound of h(d,t) on Ti. Suppose the region T is partitioned into N subregions Ti (i ) 1, ..., Nk). For each subregion Ti,,
F1i ) max min max gj(d,z,t)
(7)
g(d,z,t) ) max gj(d,z,t)
(8)
t∈Ti
z∈Z
j∈J
Also, let j∈J
Now (7) becomes
F1i ) max min g(d,z,t) t∈Ti
z∈Z
Changing the order of the first two operations in the expression, we have
F2i ) max max g(d,z,t) z∈Z
t∈Ti
There exists the relation7
min max f(x,y) g max min f(x,y) x
y
y
x
(9)
where x and y are vectors of continuous or discrete variables. Comparing the expressions for F1i and F2i and using (9), we obtain
F2i(d) g F1i(d)
(10)
Substituting g(d,z,t) from (8) into the expression for F2i, we obtain
F2i ) min max max gj(d,z,t) z∈Z
t∈Ti
j∈J
There also exists the following well-known relation,
max max f(x,y) ) max max f(x,y) x
y
y
x
where x and y are vectors of continuous or discrete variables. Using the above equality, we finally obtain
F2i ) min max max gj(d,z,t) z∈Z
j∈J
t∈Tl
(11)
From (6) and (10) we obtain
∀ t ∈ Ti
F2i(d) g max h(d,t), t∈Ti
(12)
We note that F2i(d) is an upper bound of h(d,t) in Ti . It is easy to show that eq 11 is equivalent to
F2i ) min u z,u
max gj(d,z,t) e u, t∈Ti
j ) 1, ..., m
(13)
where u is an auxiliary variable. Because eq 12 holds for each Ti, the upper bound FU 1 of F1 in the region T is of the form i FU 1 ) max F2 g F1 i
(14)
It follows directly from (11) that the following inequality holds:
F2,p(d) g F2,q(d) if Tq ⊂ Τp
(15)
Now we consider the calculation of a lower bound of the function h(d,t) on T. First, we calculate a lower
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2371
bound in the subregion Ti. The following inequality
F1(d) ) max h(d,t) g h(d,ti), t∈Ti
∀ ti ∈ Ti
(16)
holds for any point in Ti. It is known that eq 5 can be reduced to the following:
h(d,ti) ) max u z,u
gj(d,z,ti) e u,
(17)
j ) 1, ..., m
After solving (17) for any point ti, we obtain a lower bound of h(d,t) on Ti.. However, we want to obtain the tightest possible lower bound for F1i(d). Therefore, the selection of ti is critical. We use the following heuristic: we take the solution of the upper bound problem in (13) as ti. The rational for this is that the upper bound problem gives an approximation to F1i(d), which tends to F1i(d) as the size of Ti tends to zero. Because we want the tightest possible lower bound of F1(d) on T, it is reasonable to select the lower bound FL1 as follows:
FL1 ) max h(d,ti), i
i ) 1, ..., N
(18)
Figure 2. Flowchart of BB algorithm.
The new approach BB_ACTIVE for solving P1 and P2 is based on
F ) min f(d,d0)
(22)
F1(d) e b
(23)
d
We next identify the following flexibility issues: (1) P1: Determination of the flexibility of the chemical process (CP). (2) P2: Estimation of the feasibility function for a fixed design d. We note that the ULB method calculates the feasibility function without distinguishing between P1 and P2. In the case of P1, we only need to determine if the CP is feasible for the design d (i.e., the sign of the feasibility function). Preliminary computational experiments have confirmed our expectation that P1 is less computationally intensive than P2. For P1 we can use the following stopping criteria:
FL1 g 0
(19)
FU 1 ) max F2i(d) e 0
(20)
i
In fact, if eq 19 holds, then the CP design is infeasible because there exists a point at which h(d,t) is positive. On the other hand, if eq 20 holds, then the CP design is feasible according to the relation in (14). For P2, the stopping criterion is L FU 1 - F1 e
(21)
In this paper, we have used a modification of the ULB method to solve P1 and P2. This modification is referred to as the BB method as shown in Figure 2. So far, all the above-mentioned techniques for calculating the feasibility function use an enumeration procedure. Halemane and Grossmann2 and Grossmann and Floudas4 used explicit enumeration. Swaney and Grossman3 and Ostrovsky et al.6 employed an efficient implicit enumeration (branch and bound) strategy. However, the latter is exponential in nature. In the worst case, it can become an exhaustive search. To avoid this possibility, we propose a new approach, namely, BB_ACTIVE. 3. BB_ACTIVE Approach for Calculation of Feasibility Function Let us solve either P1 or P2 for d ) d0. That is, we must determine either the sign or the value of F1(d0).
where f(d,d0) and the auxiliary parameter b must be selected in such a way that the solution of (22) gives the solution of P1 and P2. We choose the following functional form for f(d,d0) n
f(d,d0) )
(di - d0i )2 ∑ i)1
(24)
where n is the dimension of d. Let d* and f* be optimal values of d and f. It is clear that the following relations hold:
f(d*,d0) > 0 if d* * d0
(25)
f(d*,d0) ) 0 if d* ) d0
(26)
First, let us consider b ) 0 in (23). Also, consider the case when F1(d0) e 0. Because there exists the feasible point d ) d0 for which the objective function f(d,d0) is equal to zero, on one hand, the minimum of the function cannot be greater than zero. On the other hand, f(d,d0) is nonnegative (see (24)). It follows from here that the optimal value f(d*,d0) ) 0 and consequently d* ) d0. If eq 23 is active at the optimal point, then F1(d*) ) F1(d0) ) 0. In this case, the solution of (22) gives the solution of P1 and P2. If at the optimal point eq 23 is inactive, then F1(d*) ) F1(d0) < 0. Here, we only obtain the solution of P1 because we only know that F1(d0) < 0 (i.e., the value F1(d0) is unknown). Now, we consider the case when F1(d0) > 0. It is clear that at the optimal point f(d*,d0) > 0. To prove the latter, suppose instead f(d*,d0) ) 0, then d* ) d0 and F1(d*) ) F1(d0) ) 0. However, this contradicts our supposition that F1(d0) > 0. Thus, we have
F1(d*) e 0 if f(d*,d0) ) 0 F1(d*) > 0 if f(d*,d0) > 0 Consequently, the solution of eq 22 gives the solution of problem P1. Consider the method of solving P2 with
2372
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
Consequently, eq 30 can be replaced by the following constraints:
F2,1(d) e b, ...,F2,l(d) e b, ..., F2,Nk(d) e b Now, eq 29 can be transformed to the form
FU,(k) ) min f(d,d0) d
Figure 3. Dependence of F(b) on b.
the help of eq 22 for the case F1(d0) > 0 (the case F1(d0) < 0 can be considered similarly). In (22), b is a parameter. Therefore, in general, the optimal value F of the objective function has to be some function of the parameter b: F )F(b). Let F1(d0) < b. Then, there exists an admissible point d ) d0 such that f(d0,d0) ) 0. Consequently, the minimum of the objective function of eq 22 will be equal to zero as well and we have F(b) ) 0. If F1(d0) > b, then using the same arguments which we used for the case b ) 0, one can show that F(b) > 0. Thus, F(b) is given by
(32)
F2,1(d) e b, ..., F2,l(d) e b, ..., F2,N(d) e b Let d(k) be the solution to this problem. Later, we will need the following theorem:9 Theorem 1: The problem
min f(x) x∈X
(33) i ) 1, ..., m
min φi(x,y) e 0, y∈Y
is equivalent to
{
0 if b g F1(d0) F(b) ) >0 if b < F1(d0)
(27)
min f(x) x,yi
φi(x,yi) e 0,
It is clear that the following equation
F(b) ) 0
(28)
has an infinite number of roots (see also Figure 3). Therefore, the calculation of F1(d0) is equivalent to the determination of the smallest root of (28). Because calculation of F(b) requires solution of eq 22, we present a solution strategy for the latter for fixed b based on an approach developed by Ostrovsky et al.6,8 for the two-step optimization problem. The new strategy takes into account the fact that the objective function in (22) does not depend on t. This method is a two-level iterative procedure that partitions the region T into subregions. The upper level serves to partition T using information obtained from the lower level. More specifically, the lower level is used to calculate the upper FU,(k) and lower FL,(k) bounds of F at the kth iteration of the upper level. We consider briefly all the main algorithms of the method. Calculation of an Upper Bound. At the kth iteration of the upper level, let the region T be partitioned into Nk subregions Tl, l ) 1, ..., Nk. At the lower level, we obtain an upper bound of F(b) as follows:
FU,(k) ) min f(d,d0) d
max F2i(d) e b,
i ) 1, ..., Nk
(29) (30)
i ) 1, ..., m
where yi (i ) 1, ...,) are new search variables. Using these statements and the expression for F2i in (11), we can transform (32) into
FU,(k) ) min f(d,d0) d,zl
max max gj(d,zl,t) e b, j∈J
l ) 1, ..., Nk
t∈Tl
Furthermore, using eq 31, we can transform the above problem into
FU,(k) ) min f(d,d0) d,zl
l
max gj(d,z ,t) e b, t∈Tl
(34)
l ) 1, ..., Nk, j ) 1, ..., m
where zl is the control variable vector on Tl. A method of solving a similar problem was discussed by Ostrovsky et al.8 Calculation of a Lower Bound. Let zl,(k) and d(k) be the solution of eq 34. Consider the following problem,
FL,(k) ) min f(d,d0) d
(35) Because eq 14 holds, any point from an admissible region, DF, of (22) belongs to an admissible region of DU of (29), and consequently (DF ∈ DU). From here the minimum of (29) cannot be less than the minimum of (22). Thus, FU,(k) is an upper bound of F(b) (i.e., FU,(k) g F(b)). The following equivalent relation holds:
max f(x) e 0 S ∀ xf(x) e 0 x
(31)
h(d,ti) e b,
ti ∈ S(k)
where S(k) is a set of points from T. Because eq 16 holds, the admissible region of eq 22 is a subregion of that of eq 35. Therefore, the minimum of (22) cannot be less than the minimum of (35). Consequently, FL,(k) is a lower bound for F(b) as follows:
FL,(k) e F(b)
(36)
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2373
Using Theorem 1 and the expression of h(d,t) in (5), one can reduce eq 35 to the problem
FL,(k) ) min f(d,d0) d,zi
max gj(d,zi,ti) e b j∈J
ti ∈ S(k)
where zi is the control variable vector corresponding to the point ti and S(k) is some set of points in T. Using eq 31, one can reduce the problem into Figure 4. Two-level iterative strategy for solving (16) for fixed b.
FL,(k) ) min f(d,d0) d,zi
gj(d,zi,ti) e b,
∀ ti ∈ S(k) j ) 1, ..., m
(37)
We can use any points from T to construct the set S(k). However, we want to select ti such that eq 37 gives a reasonably good approximation to F(b). We use the following heuristic: to construct S(k), we will employ the set of points ti obtained by solving the upper bound problem in (34). The rational is that the upper bound of (34) gives an approximation to F1i(d). This approximation tends to F1i(d) if the sizes of all subregions Tli tend to zero. Partitioning the Set of Subregions. At the upper level, partitioning of a selected set of subregions is performed at every iteration. The partitioning strategy strongly affects the computational complexity of the procedure. The simplest way is to partition all the subregions. In this case, the dimensions of eqs 34 and 37 can become very large. To alleviate this problem, we employ the following heuristic: at the kth iteration, partition only those subregions which contain at least one active constraint at the solution of eq 34. Thus, the lth subregion is partitioned if
∃ j max gj(d(k),zl,(k),t) ) b t∈Tl
We will use the following technique for partitioning a subregion. At the kth iteration, let the lth subregion have at least one active constraint. Then, the lth subregion will be partitioned into two subregions, Tp and Tq, by the addition of the constraints ts e tls and ts g tls,
Tp ) {t: ts e tls; t ∈ Tl}, Tq ) {t: ts g tls; t ∈ Tl} where tl is a point in Tl and tls is the sth component of the vector tl. For partioning, the branching point tl and the branching variable tls are employed. There are many ways to choose the branching point and the branching variable. We suggest a simple strategy. We take the midpoint of the lth subregion as tl. Next, we select one component of the variable t as a branching variable. We systematically alternate through all r components of t. The stopping criterion for solving eq 22 for fixed b is problem-dependent. For P1, the iterative procedure stops if one of the following conditions is satisfied:
|FU,(k) - FL,(k)| e 1
(40)
On one hand, if the iterative procedure terminates as a result of eq 38, then the CP is not feasible. On the other hand, if the iterative procedure terminates as a result of eq 39, then the CP is feasible. This two-level strategy is shown in Figure 4. Now, we compare the BB method and the BB_ACTIVE method. First, consider the BB method. The region T was partitioned into Nk subregions, Ti(k), at (k) the kth iteration. The upper bounds F2i (i ) 1, ..., Nk) were calculated in previous iterations. The BB method uses the following strategy. It partitions the subregion (k) (k) Tl(k) with the largest upper bound (F2l g F2l , ∀ i) into two or more subregions. Define as (k) - max h(d,t) ) F2l t∈Tl
(41)
If is small, then it is very likely that the global maximum of h(d,t) on T falls on the lth subregion. Therefore, we partition the lth subregion. Upper bounds for the resulting subregions are evaluated and a new iteration is performed. On one hand, if is not small, then the selection of Tl(k) may be wrong. Therefore, at some iteration p (p > k), the BB method will have to return to some of the subregions at previous iterations (k + 1), (k + 2), ..., p. This has clearly resulted in unfruitful calculations. This becomes worse the larger is. On the other hand, a small is expected to result in only a small amount of unnecessary calculations. With regard to the BB_ACTIVE method, suppose at the kth iteration we obtained FU,(k) by solving eq 32. Let only the lth subregion contain an active constraint. Using our heuristic, we must partition the lth subregion into p and q subregions. Therefore, at the (k + 1) iteration the upper bound FU,(k+1)will be obtained by solving the problem
FU,(k+1) ) min f(d,d0) d
(42) F2,l(d) e 0, F2,(l-1)(d) e 0, F2,p(d) e 0, F2,q(d) e 0, ..., F2,N(d) e 0 Because Tp ⊂ Tl and Tq ⊂ Tl (see (15)), the admissible region of eq 42 is contained in the admissible region of eq 32. Consequently, the following inequality holds:
FL,(k) > 1
(38)
FU,(k+1) e FU,(k)
FU,(k) > 2
(39)
Therefore, partitioning the subregion with active constraint improves the upper bound. Now, suppose that only one of the constraints corresponding to the lth
For P2, the stopping criterion is
2374
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
subregion is nonactive, which means that the solution d(k) satisfies the condition
f(d(k),d0) e f(d,d0) ∀ d (k)
(k)
(k)
F2,l(d ) ) 0, ..., F2,l-1(d ) ) 0, F2,l(d ) < 0,
(43)
F2,l+1(d(k)) ) 0, ..., F2,N(d(k)) ) 0 It follows that d(k) is the solution of the problem
min f(d,d0) d
F2,l(d) e 0, ..., F2,l-1(d) e 0, F2,l+1(d) e 0, ..., F2,N(d) e 0 Suppose we partition the lth subregion into p and q subregions, then at the next iteration we must solve the problem:
FU,(k+1) ) min f(d,d0) d
(44) F2,l(d) e 0, ...,F2,l-1(d) e 0, F2,p(d) e 0, F2,q(d) e 0, ..., F2,N(d) e 0
V
A
γ
feasibility function
I
I1 or I2
CPU time (s)
5.0 5.0 5.0 6.5 6.5 1.0 8.0 3.0 3.0
10.0 10.0 10.0 9.7 9.7 3.0 12.0 5.0 5.0
2.0 1.5 1.0 1.0 0.8 1.0 1.0 1.0 0.6
0.05648 0.04195 0.02879 0.001781 -0.001847 0.3242 -0.01528 0.09712 0.0847
156 6 2 2 1 11 12 3 4
1 1 1 1 1 1 2 1 1
128 5 2 2 1 8 7 2 1
a I, total number of iteration. I , the iteration number of the 1 first positive lower bound. I2, the iteration number of the first negative upper bound.
Table 2. Results for Example 1 by BB_ACTIVE Methoda V
A
γ
feasibility function
I
I1 or I2
CPU time (s)
5.0 5.0 5.0 6.5 6.5 1.0 8.0 3.0 3.0
10.0 10.0 10.0 9.7 9.7 3.0 12.0 5.0 5.0
2.0 1.5 1.0 1.0 0.8 1.0 1.0 1.0 0.6
0.0625 0.03125 0.03125 0.001563 0 0.3174 0 0.09375 0.07813
6 7 7 3 1 12 1 7 8
1 1 1 1 1 1 1 1 1
7 7 3 2 1 11 1 3 2
a I, total number of iterations. I , number of the first iteration 1 for FL,(k) > . I2, number of first iteration for FU,(k) < .
Because Tp ⊂ Tl and Tq ⊂ Tl, we have
F2,l(d) g F2,p(d), F2,l(d) g F2,q(d)
Table 1. Results for Example 1 by BB Methoda
(45)
Consequently, the constraints corresponding to the p and q subregions will be nonactive again and we obtain
FU,(k+1) ) FU,(k) Thus, partitioning subregions with a nonactive constraint has no effect. Therefore, we employ the heuristic given earlier, on one hand, to avoid unnecessary partitioning in contrast to the BB method. On the other hand, the BB_ACTIVE method is more computationally intensive at every iteration. Thus, one can expect that the BB_ACTIVE method will be superior to the BB method if the upper bound algorithm used in the BB method gives a poor approximation to the maximum of h(d,t) on each Ti (i ) 1, ..., Nk). The reverse is true if the upper bound calculation for BB gives a good approximation to the maximum of h(d,t). Next, we will present three examples to illustrate the application of the BB_ACTIVE method and to compare it to the BB method. 4. Computational Experiments To compare the BB and BB_ACTIVE methods, we present results from both methods. All the calculations were performed on a DELL Pentium II 400-MHz personal computer. For all the examples, we represent the uncertainty region in the form N T(γ) ) {ti: tN i (1 - γβi) e ti e ti (1 + γβi)}
where tN i is the nominal value of the uncertain parameters, γ is the parameter to determine the size of the uncertain parameter range, and βi is a deviation fraction. Example 1. This is the problem illustrated in section 1 (Motivational Example). We investigated the chemical process (CP) flexibility for different values of the design
variables and for different sizes of the uncertainty region. Here, the vector of βi’s is [0.1, 0.02, 0.03, 0.1, 0.1]. The results from the BB and BB_ACTIVE algorithms are shown in Tables 1 and 2, respectively. From Tables 1 and 2, we can see that the flexibility function increases when the range of uncertain parameters is increased for a fixed design. To compare the BB and BB_ACTIVE methods, let us consider the first row in the tables. For a fixed design, V ) 5.0 m3, A ) 10.0 m2, and when γ ) 2.0, the uncertain parameters vary as follows: 36.288 < F0 < 54.432, 319.68 < T0 < 346.32, 282 < Tw1 < 318, 7.84 < kR < 11.76, and 1308 < U < 1962. The corresponding value of the flexibility function is 0.05648 for the BB method; this required 156 iterations and 128 CPU s (the C++ time function time( ) was used to calculate the CPU time). We note that, at the first iteration, the lower bound of 0.01453 is positive, which means that the process is not feasible for this design. The BB_ACTIVE method obtained similar results in only 6 iterations and 7 CPU s. Tables 1 and 2 also show the results for different designs and ranges of uncertain parameters. Up to the second significant digit after the comma, both methods gave the same values of feasibility function. The results consistently favor the BB_ACTIVE method over BB. It is interesting to note that the method of Halemane and Grossmann2 required 73.9 s (on a DEC-20 workstation) for evaluating the feasibility function (V ) 6.6 m3, A ) 9.7 m2, γ ) 1.0), while the BB and BB_ACTIVE methods required 2 s. It is not clear how much faster the DELL Pentium II 400 MHz is compared to a DEC20. Therefore, a direct comparison with a reported CPU time is not possible. However, if we assume that the Pentium II is an order of magnitude faster than the DEC-20, then we can still conclude that we have developed an efficient algorithm. Example 2. In this example, we will consider the CSTR optimization problem as shown in Figure 5 with the nomenclature presented in Table 3.10 Two first-order
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2375 Table 4. Results for Example 2 by BB and BB_ACTIVE Methodsa BB feasibility function
γ
I
0.1 U ) -27.81; 200 L ) -50.4 0.2 U ) 2.89; 200 L ) 1.89 0.3 U ) 65.42; 200 L)63.67
BB_ACTIVE I1 CPU feasibility or I2 time(s) function
I
I1 CPU or I2 time (s)
1
289
0
1
1
1
2
215
2.25
6
1
4
5
229
64.87
16
1
16
a I, total number of iterations. I , the iteration number of the 1 first positive lower bound. I2, the iteration number of the first negative upper bound.
Figure 5. Flowsheet of Example 2. Table 3. Nomenclature of Example 2 V ) reactor volume, m3 F ) volumetric flow rate, m3/min Q ) rate of heat removal, J/min CA0, CB0 ) inlet concentrations of A and B, mol/m3 T0 ) inlet temperature, K PB ) rate of production of B, mol/min EA, EB ) activation energies, J/mol k0A, k0B ) pre-exponential constants, min-1 ∆H1, ∆H2 ) molar heats of the reactions, J/mol F ) system density, kg/m3 Cp ) system specific heat, J/(kg K) R ) ideal gas constant, J/(mol K)
reactions in series (A f B f C) with Arrhenius-type kinetics are assumed. The volume of the reactor is taken as a design variable. We analyze the process flexibility for V ) 0.391 m3. Volumetric flow rate F, reaction temperature T, feed temperature T0, initial concentrations of A and B, CA0 and CB0, are the designated control variables. The bounds for these control variables are 0.9036 e F e 1.1044, 280.44 e T e 342.76, 298.3 e T0 e 329.7, 2962.1 e CA0 e 3273.9, 325.375 e CB0 e 359.625, respectively. Uncertain parameters t ) [EA, EB, k0A, k0B, ∆H1, ∆H2, F, Cp], and the nominal values tN ) [3.64 × 104, 3.46 × 104, 8.4 × 105, 7.6 × 104, 2.12 × 104, 6.36 × 104, 1180, 3.2 × 103]. The deviation for each uncertain parameter is equal to [0.1, 0.01, 0.1, 0.1, 0.05, 0.05, 0.05, 0.05]. To eliminate state variables CA, CB, and Q, we use the following: (1) Energy balance:
FCpF(T0 - T) + k0A exp(-EA/RT)CA(-∆H1)V + k0B exp(-EB/RT)CB(-∆H2)V - Q ) 0 (2) Material balance equation:
F(CA0 - CA) - k0A exp(-EA/RT)CAV ) 0 F(CB0 - CB) + k0A exp(-EA/RT)CAV k0B exp(-EB/RT)CBV ) 0 We impose constraints on the production of B and the heat duty Q as
g1 ) -PB + 549.6 e 0 g2 ) PB - 650.4 e 0 g3 ) Q - 2.79e7 e 0 g4 ) -Q + 2.286e7 e 0
Figure 6. Flowsheet of Example 3.
where
PB ) FCB ) FCB )
FCB0 + k0A exp(-EA/RT)VCA F + k0B exp(-EB/RT)V
and
Q ) FCpF(T0 - T) + k0A exp(-EA/RT)CA(-∆H1)V + k0B exp(-EB/RT)CB(-∆H2)V The results from the BB and BB_ACTIVE methods are shown in Table 4. Considering the first row in Table 4, we see that for γ ) 0.1, the process is feasible. When γ ) 0.2 and 0.3, the BB method could not converge on the solution after 200 iterations. Therefore, we present the upper and lower bounds of F1(d). In contrast, we obtained accurate results from BB_ACTIVE in much fewer iterations and smaller CPU time. Although the BB method did not converge after 200 iterations for some γ, we know the process is not flexible within that range of uncertain parameters after 2 or 5 iterations. For example, at γ ) 0.2, after 2 iterations, the lower bound of the flexibility function is equal to 1.8946, and at γ ) 0.3, the lower bound of the flexibility function is 63.674 after 5 iterations. Example 3. In this example, we investigate the flexibility of the heat exchanger network as shown in Figure 6. The optimized network for processing two hot streams (FS2 and FS3) and three cold streams (FS1, FS4, FS5) were obtained by Masso and Rudd.11 We define the
2376
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
heat exchanger area, cooler and heater as design variables, and study the fixed design at [A1, A2, A3, A4, AC, AH] ) [213, 764, 317, 60, 143, 12] (ft2). The flow rates for cooling water and hot steam are taken as control variables with bounds of 20 000 e Fwc e 40 000 and 100 e Fh e 1000. The inlet temperature and heat capacity of each stream are considered as uncertain parameters with nominal values tN ) [TS1, TS2, TS3, TS4, TS5, CpS1, CpS2, CpS3, CpS4, CpS6] ) [100, 440, 520, 200, 350, 0.8, 0.7, 0.68, 0.85, 0.62] (lb/h). The deviation for each uncertain parameter is as follows: [0.01, 0.01, 0.01, 0.01, 0.01, 0.001, 0.001, 0.001, 0.001, 0.001]. Other data for the process are given in Table 5. The outlet temperature for each stream is allowed to vary by (10 °F, and the minimum allowable stream temperature difference for each heat exchanger, cooler, and heater is applied. Therefore, we have the following constraints (gj e 0):
g1 ) -Ts5,2 + 395 g2 ) Ts5,2 - 405 g3 ) Twc2,2 - 180 g4 ) -Ts2,2 + 155
Table 5. Specific Data for Example 3 stream no.
flow rate (lb/h)
input temp. (°F)
output temp. (°F)
heat capacity (Btu/(lb °F))
1 2 3 4 5
20 000 40 000 35 000 31 000 32 000
100 440 520 200 350
430 150 300 400 410
0.80 0.70 0.68 0.85 0.62
g7 ) Ts6,2 - 415 g8 ) -Ts6,2 + 405 g9 ) Ts4,2 - 355 g10 ) -Ts4,2 + 345 g11 ) Ts7,2 - 155 g12 ) -Ts7,2 + 145 g13 ) Ts3,2 - 305 g14 ) -Ts3,2 + 295 g15 ) Ts1,2 - 435 g16 ) -Ts1,2 - 425 g17 ) Ts4,2 - Ts7,1 +20 g18 ) Twc1,2 - T1 + 20 g19 ) T2 - T6 + 20 g20 ) Ts5,2 - Ts2,1 + 20 g21 ) Twc2,2 - T3 + 20 g22 ) Ts1,2 - T4 + 20
tw ) 180 °F
Minimum Allowable Approach ∆T’s heat exchanger ∆THE ) 20 °F steam heater ∆T H ) 25 °F water cooler ∆T C ) 20 °F Overall Heat-Transfer Coefficients heat exchanger uHE ) 150 Btu/(h (sq ft) °F) steam heater uH ) 200 Btu/(h (sq ft) °F) water cooler uC ) 150 Btu/(h (sq ft) °F) steam pressure Ps ) 450 lb/sq in abs., saturated (Hv ) 1206.1 Btu/lb, HL ) 438.7 Btu/lb, latent heat ) 767.5 Btu/lb, T ) 456 °F) cooling water temperature twI ) 100 °F Table 6. Results for Example 3a
g5 ) -Ts2,2 + 145 g6 ) Tw,c12 - 180
Cpw ) 1 Btu/(lb °F)
heat capacity of cooling water maximum water output temperature
BB
BB_ACTIVE
γ
feasibility function
I
I1 or I2
0.5 0.8 1.0
0.2857 0.4696 0.5916
98 95 87
20 29 29
CPU feasibility CPU time(s) function I I1 or I2 time (s) 900 840 1140
0.2969 0.4766 0.6016
9 9 9
1 1 1
6 8 7
a I, total number of iterations. I , the iteration number of the 1 first positive lower bound. I2, the iteration number of the first negative upper bound.
when F1 was evaluated, and at the 29th iteration, the lower bound becomes positive for the first time. On the other hand, the BB_ACTIVE method required only 9 iterations. At the first iteration, the lower bound was positive, so we concluded that the network was infeasible for the given design d. From analysis of the results of the flexibility analysis of the above three examples, one can draw the following conclusions: (1) Both methods required very few iterations for solving P1 (in most cases, only 2 iterations). (2) In some cases, the BB method either required a large number of iterations for evaluating the flexibility function or did not converge to the solution after 200 iterations. The larger the number of iterations, the larger the storage of information about the subregions. The BB_ACTIVE method does not appear to have this problem because it requires a few iterations for each b. (3) In most cases, BB_ACTIVE shows a large advantage over the BB method when the flexibility function is evaluated.
g23 ) T5 - Ts3,1 +20 g24 ) Ts6,2 - Tsteam + 25 We now analyze the flexibility of the heat exchanger network for different sizes of the uncertainty region and present the results as shown in Table 6. From Table 6, we observe that, for all the values of γ ) (0.5, 0.8, 1.0), the feasibility function is positive and consequently the network is not feasible for the given design and uncertainty regions of these sizes. For γ ) 1.0, the BB method, on one hand, required 87 iterations
5. Conclusions We have considered issues regarding estimation of the flexibility of chemical processes through the feasibility function. We considered two problems, namely, the determination of whether a process is feasible (i.e., the sign of the feasibility function) (P1) and the evaluation of the feasibility function (P2). We have developed a very effective method, BB_ACTIVE, that does not employ any enumeration procedure for the solving of the flexibility problem.
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2377
In most cases, BB_ACTIVE shows a very significant advantage over the BB method with regard to the process flexibility function calculation. Acknowledgment We gratefully acknowledge funding from the NSF under Grant CTS-9726135. Literature Cited (1) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall: Englewood Cliffs, NJ, 1997. (2) Halemane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AIChE J. 1983, 29 (3), 425. (3) Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. AIChE J. 1985, 31 (4), 621. (4) Grossmann, I. E.; Floudas, C. A. Active Constraint Strategy for Flexibility Analysis in Chemical Processes. Comput. Chem. Eng. 1987, 11 (6), 675. (5) Kabatek, U.; Swaney, R. E. Worst-Case Identification in Structured Process Systems. Comput. Chem. Eng. 1992, 18, 11063.
(6) Ostrovsky, G. M.; Volin, Yu. M.; Barit, E. I.; Senyavin, M. M. Flexibility Analysis and Optimization of Chemical Plants with Uncertain Parameters. Comput. Chem. Eng. 1994, 18 (8), 755. (7) McKinsey, J. C. C. Introduction to the Theory of Games; McGraw-Hill: New York, 1952/ (8) Ostrovsky, G. M.; Volin, Yu. M.; Senyavin, M. M. An Approach to Solving the Optimization Problem under Uncertainty. Int. J. Syst. Sci. 1997, 28 (4), 379. (9) Ostrovsky, G. M.; Volin, Yu. M.; Golovasuhin, D. V. Optimization Problem of Computer System under Uncertainty. Comput. Chem. Eng. 1998, 22, 1007. (10) Bernardo, F. P.; Saraiva, P. M. Robust Optimization Framework for Process Parameter and Tolerance Design. AIChE J. 1998, 44 (9), 2007. (11) Masso, A. H.; Rudd, D. F. The Synthesis of System Designs. AIChE J. 1969, 15 (1), 10.
Received for review July 14, 1999 Revised manuscript received March 15, 2000 Accepted March 20, 2000 IE9905207
