- Email: [email protected]
Worst-case identification in structured process systems
0098-l 354/92 SS.00+ O.&l copyright Q 1992 Pergsmon F%as Ltd
Compufers hem. &gng, Vol. 16, No. 12, pp. 1063-1071,1992 Printed in Great Britain. AU rights mservcd
WORST-CASE
IDENTIFICATION IN STRUCTURED PROCESS SYSTEMS U. I(ABATBKt and R. E. SWANEY$
Department
of Chemical Engineering, University of Wisconsin, Madison, Wl 53706, U.S.A.
(Received 20 September 199o;jnal reuiGm received 14 May 1992; receivedfor publication 28 May 1992) A--An improved method for identifying the worat-caae set of parameter combinations in controllable systems is presented that is based on branch and bound search in the space of the uncertain parameters. A characteristic of controllable structured systems frequently appears wherein the limiting condition is local&d in a part of the system that has only limited coupling to the other parts. While there may be many uncertain parameters varying in parta upstream of the limiting constraints, often there will remain suf&ient control dcgoxs of freedom in those other parts to compensate. The branch and bound method in ita simplest form cannot recognize this situation, and will usually engage in searching among combinations of such parameters, when in the end they have no effect on the limiting condition. Modifications to the basic method are presented that allow these parameters to be nxagnizcd apriori. The complexity of the branch and bound search is reduced accordingly, while the rigor of the original method is maintained. As shown in three heat exchanger network examples, a substantial reduction in effort is often achieved.
INTRODUCTION
It is a reality of practice that chemical processing systems cannot simply be designed for some fixed nominal set of conditions, but must be designed to successfully accommodate a range of conditions. The origins of these variations fall into two categories. In the first are expected variations in operating conditions such as feed and product specifications, catalyst activity, utility supplies and equipment fouling. In the second are the uncertainties inherent in the design models themselves and in the physical data employed. While the former are known to vary, and the latter do not actually vary but are unknown, both types of uncertainty impose requirements for the design to meet. The task of design under uncertainty is to ensure that the worst case within the range of expected variations (or realizations) can be successfully handled. This requires that the worst-case combination be identitied, and then designed for. Worst-case identification in a complex process system often is not a trivial task. The existence of control adjustments that can be made in response to variations contributes to the complexity of the analysis. A formal framework for addressing these concerns was presented by Swaney and Grossmann (1985a) in the form of a “flexibility index”. In that framework,
the problem variables are classiiied as: design variables d, which am selected during design but are fixed during operation; control variables E, which may be varied to suit circumstances during operation; uncertain parameters 0, which exhibit a range of values during operation that must be accommodated, and state variables x, which are dependent variables in the process model. The ranges of variation in the uncertain parameters %,are described by deviations A@;, Aii,? below and above a nominal operating point 07. The index value F measures the fraction of the target ranges of variation that the design can accommodate for simultaneous independent variation of the parameters. Process operating constraints and product speci&tions are specifM by the set of inqualities g(d, z, 0, x) C 0 while the process model is represented by the set of equations h(d, 2, 8, X) = 0. The latter are suflicient to determine the state variables x. While generally the variables x and the relations h( . , . , *, x) = 0 are r&ah& explicitly in implementation, for discussion it is conveniat to treat the state variables as implicit functions x = h -‘(d, z, 0) and to eliminated them. The constraint functions may then be expressed in the reduced form g(d, z, 0, x) =f(z, 6), with the dependence on the design variables implied but suppressed. The analysis problem may be formulated as:
(1)
s.t. .f(z, e) < 0,
pmellt addreux lnstitut mrsystcmdynamik und r&geltecbnik, ulkniitit stmtgart, PfafFenwaldkg Stuttgart go, Germany. ~To whom all correspondence should be w.
F=mjnyS
9. D7ooo
8 =P+88,
(2) (3) (4)
1063
1064
U. KABATEKand R. E. SWANEY
to represent whether or not a constraint is active. For a given combination of these, the maximization part of (1) becomes implicit, so the problem becomes one of minimizing over 6 and the binary variables. The solution is then found either by solving the resulting mixed integer program (an MILP for the linear case), or by enumerating the feasible combinations of active constraints and solving the NLP in 8 for each combination. Since this method searches among L=max 6, (5) the constraints, it may be described as a “dual” d* method. For further applications of the active cons.t. f;(z, el) 6 0 (6) straint approach, see Pistikopoulos and Grossmann i=l m, (1989). @=@N+@ ‘*“’ (7) > Regarding relative merits, the two approaches offer somewhat different capabilities. Although both By evaluating each constraint at its individual worstprocedures are rigorous and effective for purely linear case, 6’= 8**‘, the solution to this NLP will provide models, the active set method usually requires fewer a lower bound to (1). binary decisions and fewer underlying combinations. The basis of the implicit enumeration algorithm is The dual approach is thus well-suited for linear to employ a branch and bound procedure to search problems, wherein its strengths can be fully exploited. the domain ee(-AO;,AOT}, j= 1,. . . , no, to For complex nonlinear problems, the primal approach locate the worst-case &*. At a node in the search presents its own advantages. First, the primal method tree, each parameters 4 is either “unassigned” or is easily implemented for complex system models. “assigned”. Unassigned parameters j e J take on the by a small set of values in the individual constraints given by 8; = &;CJ. The subproblems are representable NLPs loosely coupled in a multipcriod manner, each Assigned parameters j E J, on the other hand, are of which is a parametric optimization of the confixed either at 8,= -AO; or 8/ = AO:t and take ventional process model. It is therefore not difficult on the same value in every constraint. The tree is to adapt existing flowsheet optimizers to this task. developed by taking a node with the lowest lower Equivalent simplicity of implementation is obtained bound, selecting one of its unassigned parameters to for the dual approach only if the method of enumerbecome assigned, and replacing the node by separate ating active sets is employed. Moreover, a feasible descendant nodes, one for each of the possible starting point is available for every NLP calculation assignment values. This amounts to partitioning the required by the primal method. Secondly, since the search domain according to the value of the assigned resulting solution is based on a series of primal-feasiile parameter and computing individual bounds for calculations, any impact of general nonlinearities on the partitions. The validity of the collective bound is the global validity of the result is likely to be conmaintained, since for every assigned parameter, each servative, in the form of suboptimal control settings. extreme value is represented by a node in the tree. Global rigor rests on the ability to locate parameter The procedure is continued until enough parameters values that are worst-case for each constraint indihave been assigned to close the gap between the lower vidually, and more intensive methods may be applied bound and the solution. The worst-case parameter to this subtask if needed. The impact of nonlinearities combination 8* is then readily identied within in the dual formulation is less clear. Finally, it should the partition having the lowest bound. Since this be noted that both approaches usually are capable procedure executes the worst-case search directly in of locating nonvertex solutions (i.e. 8: interior to the variables 8, it may be. described as a “primal” [ - AO; , A@; ] for some j) when they occur in typical method. problems. In the primal method these are usually The second approach is the “active constraint detectable and accessible by descent from a vertex strategy” outlined by Floudas and Grossmann using a local minimization. (1987). The idea is to search among the possible The purpose of the present paper is to present an combinations of constraints that will be active at improved branch and bound strategy for the primal the worst-case condition. Binary variables are defined method. The modified method permits a reduction in the number of branches in the search tree by recog?A generalization of this exists for solution points that arc nizing the presence of “compensatable parameters” not at the parameter bounds; a branch mightconceivably within the problem. The resulting reduction in solugenerate an assigned value 8, interior to [ - AO;-, AfI,+1, possibly with 4 = 4(z, 0). tion effort can be substantial in typical problems. The solution O* = ON + 6*8* represents the worstcase operation at which limiting conditions are reached. Two main approaches have been proposed to solve (1). The first is the implicit enumeration method presented by Swaney and Grossmann (1985b). This method is based on a lower bound program of the following form:
1065
Worst-ease identiticatiou in structured process systems F5 T14 F3
T15
I T7
I-T10 T5
+ Tll F4
F,n+$
T2
($3
I
T9
i Fig.
1.
T13
Example I network.
COMPENSATABLE
PMMGXERS
This term is introduced to describe a behavior frequently exhibited by structured process models. O&en there are a number of parameters whose values have no impact on the feasibility of the worst-case limiting condition. This situation is principally a commquence of system spar&y. Limiting conditions usually arc encountered in only a part of the system (e.g. where the “bottleneck” is located), while other parts operate inside of their limits and still have degrees of freedom for compensation. An illustration of this behavior is provided by Example 1, shown in Fig. 1. The data describing this small heat exchanger network am given in Tables l-3. Heat transfer coe!Bcients vary with flowrate and
Table 1.
Exampla1 mwrtain
e,
i :
e/”
r{ (m% kW-‘) rr (m’K kW-‘1
r) &K kW-'j
0.5 0.5
r: (mzK kW_‘) F, &t~a-‘) Fa (ks a-‘) JS &I-:) ? (s)
ix 10’ ii 5:
T: (IO Tl(K) T,, IKl
723 388 313
parametera -AO;
ae; +0.5 +0.5
+0.5
+0.5 +1.0 +2.0 +2.0 +3.0 +5.0 +5.0 +5.0 +s.o
-0.5 -0.5 -0.5 -0.5 -1.0 -2.0 -2.0 -8.0 -5.0 -5.0 -5.0 -5.0
Table 3. ii’ (kWm”K-I)
k
constant heat capacities of I kJ kg-’ K-’ apply for every stream; details are given in the Examples section. In this example, flowratc and inlet temperatures for streams Fl-F4 are treated as uncertain parameters, along with the fouling resistancea in each exchanger. Each exchanger has a variable bypass; together with the split fraction controlling the flow into E4 they comprise the five control variables. Outlet temperatures must meet inequality con&mints. The limiting condition in this example occurs in the lower part of the network around El and E3, with El operating at zero bypass. Of the 12 uncertain parameters, only Fl, F3, F4, Tl, Tl2 and the fouling in El and E3 impact the limiting condition. The original branch and bound procedure generates 57 nodes in solving this problem, with branches on six of the variables. That method is inherently able to avoid branching on the four fouling resistances or on F4 or T12, since in this problem they have uniform worst-case values in each of the limiting constraints. However, branches are made on the uncertain parameters F2, T4 and T7. Variations in these parameters in fact do not effect the result. This is a consequence of the structure of the network and of the fact that control degrees of freedom remain in the upper part of the network. Temperature TS is the only connection between variations in these parameters and the limiting constraints, The upper part of the network is able to adjust T8 to the optimum value required by the constrained lower part, in spite of variations in F2, T4 and T7. These parameters are thus termed compensatable. In the original method, branching paranmters are selected by computing the sensitivity of the node lower bound program to their assignment. Although the solution node is insensitive to changes in a compensatable parameter, intermediate parent nodes
Example1 cxdunspr e (kgr-‘)
data
fi: (%Wm-*K-l)
e (kgs-‘)
(m’K$W-‘)
:
ii
1.5
20 10
1.5
iti
0
3 4
120 12
1.5 1.5
ii
1.5 6.0
30 10
: 0
U. lCm~mx
1066
and R. E.SWANEY
may respond significantly to aaaignment of the parameter. This behavior results when the individual constraint worst-case values 87.’ are not uniform for the parameter. The effect of the compensatable parameters on the node lower bounds does not disappear until their assignment makes them uniform in the active constraints. The presence of such parameters artificially complicates the worst-case identification task for the original method. Through proper recognition of compensatability, however, the added complexity can be essentially eliminated. In the case of Example 1, the number of nodes in the branch and bound tree reduces from 57 down to 9. MODIFIED
METHOD
The modified method seeks to avoid unnecessary branching by recognixing the presence of compensatable parameters a priori. Since there is no general a priori way to identify such parameters with certainty, the method works with a projected (postulated) set of compensatable parameters, which is revised as needed during the course of the algorithm. The projected set is then either confirmed or refuted when branching appears complete. This approach permits efficient solution with nearly a minimum of branching, yet retains the rigor of the earlier method. In the following, frequent use is made of linearized perturbation analysis of the solution of the node lower bound program (5). Of particular use are the predicted changes in the bound accompanying changes in the vertex values C. Let 12 0 be the vector of multipliers for the constraints (6) at the solution of the lower bound program (5) for the node. For the changes in parameter 4 from the values 8: at the no& to the uniform assignments A6: and -A@;, the respective bound change predictions are given by:
The projected set of compensatable parameters is generated as follows. At a node k early in the tree, parameters Rj, j lJL are assigned to fixed vertex values g):, while the unassigned parameters 0;, j E Jk taken on the various individual worst-case values @’ for each constraint i. Since parameter compensatability is not usually evident until all intluential parameters are assigned values that are uniform for every constraint, a complete assignment is selected. A likely candidate for the worst-case vertex within the node’s subtree may be predicted by selecting
assignments that give the best predicted bound increase according to (8,9). A complete assigmnent is thus obtained by setting the unassigned parameters according to:
y=fJ;,
jeJL
(10)
(5) is then solved for CT’= 6k for this projected vertex. At this vertex solution, parameter sensitivities are again computed from (8,9). Compensatability is readily evident, and the projected set of compensatable parameters C” for node k is taken as:
C*={jly,Z=y;=O,jeJk}.
(11)
The compensatability
property of a parameter derives from and depends on the set of constraints that are active at the worst-case vertices. The active set associated with the projected set C’ is also indicated at the vertex solution by {iI&> 0}, and will be denoted A(Ck). With a projected set identified, the branch and bound procedure continues as in the original method, except that parameters in C* are reserved from selection as branching parameters. Branching parameters are selected from the complement set Cc. Each new subnode k’ generated originally inherits the set Cc’ = C* of its parent. However, the eventual applicability of a set C” to vertices below node k’ is likely to be contingent upon having the same active set that generated Cy also be active at those vertices. For this reason, if the active set I$’ exhibited at the solution to the node’s lower bound program differs from A(C’), a new set of compensatable parameters is projected for node k’ using the procedure given above. C’ is then taken as the new projected set and branching continues. Note that the use of projected sets CK in these steps implies no assumptions, but merely serves to guide the branching strategy. When a node z is produced by the last branch indicated within the set Cc of its parent node k (i.e. no further branching would be indicated in c), the compensatable parameter set Cr = Ck is usserred. In this step, the parameters in C” are assigned to uniform values 8, in the active constraints i E A (CE). while retaining their individual worst-case values 8:’ in the inactive constraints. This action seeks to merge the combinations of the parameters in Cr into a single node. The bound computation (5) is performed for E with Cr asserted. At the solution y,+ , y,- are determined from (8,9). The asserted set is then con8rmed if: yj’=y;=O,
VjeC4
(12)
Worst-ease identificationin structuredproasa systems If the asserted set is continned, the bound computation for z represents a valid lower bound, and no further branching at the node is required. If the set is refuted, however, the projected set is abandoned. A new projected set is determined for the node as above. The bound computation is recomputed with this second set asserted. If the second set is also refuted, the lower bound is again recomputed with the set. not asserted, and node E remains open for further branching. Confirmation criterion The purpose of the confirmation test is to verify that variations of the compensatable parameters within the limits -AB; < 4 6 AtI,? do not a&ct the value of the lower bound computation (5). SuIIiciency of the criterion (12) derives from the following assumption: Regularity condition-At the solution to the lower bound (5) for node E with C” asserted, let Rz, 0) be the vector of ti active constraint functions [i.e. x(z, 0) =&(z, 0) for i = 1, . . . , ff; i' EZ:],with corresponding multipliers 1, > 0. Then the regularity condition holds if:
rank
(13)
over the feasible values of 8 and z. Stationarity conditions at the solution require that:
to (5)
which means the rank of g/azr cannot exceed R - 1, The regularity condition requires the constraint gradients to have maximum independence in the control variables within this limit. The condition will hold at the solution point except in cases with nonunique multipliers, which could arise either from redundancy in the constraints or unlikely coincidence. A simple remedy is to stipulate a solution with the minimal number of positive multipliers, with A mdurxd accordingly. The extension to neighboring values of 0 and z is a modest condition that will nearly always hold in practice. Under the regularity assumption, the null space of the active constraint gradients in the controls is of rank one, spanned by x(z, 0) that is defined implicitly by x’(aJlaz’) = 0 over the feasible values of z and 8. Consider alternate values P(z), s e [O, 11, where B’(0) = 6r is the value at the solution to (5). Under the requirement that: Ix(z, t9]f$(z,
8.)]%
= 0,
(15)
1067
over e’(s) and z(s), the imphcit function theorem provides the existence of z(z) such that nz(z), e'(z)]= 0. Variations in the compensatabk parameters correspond to various 0; (s), j E C’, while holding 0;) je CI, constant. The existence of the condition r; =r; = 0 in (12) implies: XTg=O,
VjJj@,
I
at the solution to (5), and so provides that condition (15) will hold at that point. There is nothing inherent in the formulation of (5) itself to induce condition (16). Thus, excepting extreme coincidence, the appearance of condition (12) indicates a structural characteristic of the system. This characteristic exhibited byxz, 0) will not disappear with variations in e,, j E d, so requirement (15) may be expected to hold over these variations. The net result of this argument is that the control variables will possess sufficient degrees of freedom to maintain: l;(z,e) = 0, ie z: (17) for the active constraints over variations in the compensatable parameters. No change is implied in the objective function variable S. The remain@ requirement for validity of the lower bound with the compensatable parameters asserted is to ensure feasibility of the inactive constraints. This is immediately provided, since they are evaluated at their individual worst-case values @’ for / E Cxc Jr. An asserted parameter set may thus be confirmed or refuted by checking condition (12).
Detine
Jx = index set of unassigned parameters at node k, Jk = index set of assigned parameters at node k, K = index set of open nodes, Bk = index set of branch candidate parameters, #Bk=cardinality of set Bk. Branch and bound 0. set
P=0, co=0, A(CD)=0, +&.ceat,,==nb, K = (0). Evaluate node 0.
1. Detelmine A ==argr${L,}.
2. Branch on parameter I = i(k). (a) if #B,t-0, Stop; F=t,.
U. KAMTEK and R. E. SWANEY
1068
(b) form subnodes
k’ and k”:
6:=85=8;, &‘=
1. Project vertex 9, j e Jk, using (10) with y+. yfrom node evaluation. 2. Solve (5) with 8’= dk and obtain active set A(CI). 3. Compute y:, y; at projected vertex using (8,9). 4. Determine Ck from (11).
je:JL,
-A@,-,
&=A0:,
~=‘=JQ{Z}, A (CK) = A(P)
Parameter set projection @ode k)
CK=CY=Ck, = R (C”).
3. Evaluate node k’. 4. Evaluate node k”. 5. Replace node k by nodes k’ and k” in K, and go to step 1.
Node evaluation (node k’, parent k) 0. If # Bh = 1 perform Step 1. Otherwise perform Step 2 with ilag not set. 1. (Assert C”): (a) select #yet-Ae-,Ae+]
forjoCY;
for (j d P) or [j E Ck’, i E A (F)], @*’ for (j e Jv n c’) or [j e CK, i 4 A(P)];
VT (b) set fl=
(c) solve (5) to obtain LK; (d) compute Y_T , yl from (8,9) for j e JK; = 0 for j E CK and min{Y,+, y;} = 0 (e) ifY/‘=Y; for j E JK, them C’ is continned and LK is valid; node evaluation is complete; (f) otherwise, project a new parameter set and repeat Steps (a-e) with the revised Cy; (8) if the revised Cy is not confirmed, perform Step (2) with flag set. 2. (CK not asserted):
solve (5) -to obtain & and Z$; if A (Ck’) $ Zs,project a new parameter set and revise c”; compute y,+, y,- from (8,9) for j E J”; determine branch candidates within set M = JKnE’:
W3 = am~IY,}, ?I=
-4 {
(18)
IY,+ r,-l, a) -
&={j]Y~)O,
ifmin{v;,7;) r otherwise,
0
jeM}
0, (19) (20)
(f) if # BK > 0, node evaluation is complete; (g) otherwise, repeat (e) for M = Jy; (h) if redetermined #BE = 0 or flag is set, node evaluation is complete; (i) otherwise go to Step 1.
In Step 2a of the branch and bound procedure, one may wish to project to a vertex to obtain a complete assignment, and then evaluate the vertex in order to confirm the solution. This will deal with the rare chance in which a gap remains when the branch selection indicates no candidates. Further branching would be performed if required. At the solution vertex it is of value to test for local descent in the function 6 *(0) to check for a nonvertex solution accessible by local descent from the vertex. The test is performed by checking for negative values of A’tTf/atl,(-8,) at the vertex solution. A local minimization of S*(B) may then be performed if so indicated. The simplest way to obtain individual constraint worst-case values 8.’ is to assume monotonicity in each parameter and to select - A0,- or Af?: based on the signs of the gradients @J&9,. In this case y,+, y,should be monitored for the appearance of negative values, which would indicate a failure of the monotonicity assumption. In response, duplicate constraints with multiple 8**’at both extremes can be employed. Alternatively, a local search may be performed to locate a possible nonvertex &I*‘, with duplicate constraints then employed as indicated. EXAMPLES
Comparative results for the mod&d method are presented here for three heat exchanger network examples. In these examples, heat transfer coefficients are functions of flowrate, and enthalpies are quadratic functions at temperature. In each problem stream flowrates and inlet temperatures, and exchanger fouling, are uncertain parameters. Each exchanger has a variable bypass, and stream outlet temperatures are specified by inequalities. For this formulation the reduced constraints are monotonic in the uncertain parameters, giving certain determination of the P. The modeling equations may be found in the Appendix. For Example 1 described earlier, the modified method requires nine nodes, reduced from 57 with the original procedure. The set of compensatable parameters is correctly identified by projecting at tbe top node. The active set remains unchanged throughout the tree, so the projected set remains in use throughout.
identihtion in strucked processsystnns
Worst-case
1069
Tabk~4.Eumpk12unsrtlin~ e,
/ 1
2 3 4 5 6
F7
F4~l~& (&
t: 13 14 15 16 17
Tll
Jw
-A$
1 t::
+0.75 +1.0 +0.9 +1.0 +1.0 +0.7 +1.0 +1.0
-0.75 -1.0 -0.9 -1.0 -1.0 -0.7 -1.0 -1.0
ii 30 30 50 553 373 323 423 373 333 353
+1.0 +I.0 +1.0 + 1.0 +5.0 +5.0 + 5.0 + 5.0 f5.0 fS.0 + 5.0
-1.0 -1.0 -1.0 -1.0 -5.0 - 5.0 -5.0 -5.0 -5.0 -5.0 -5.0
0.75 1.0 0.9
r{ r5(m (m;K KkW-) kW-:) I, (m* K kW-‘)
ri (III*K kW-1.
5 9 10
I
e/w
T19
I (m*KkW. -I rrI (m’K kW, -‘) F, Or6s-‘) 4 Orrs-‘) F2 (Ltl a-‘) F. (Ls 8-l) 4 (kg s-v &ckes-:)
9 FrT: (K) T, (K) T,ci(K) T,, (K)
0.7 m 8
)
Tn. (K) T;; *)
20
Example 3
F3
Fig. 2. Example 2 network. Exampl% 2 The network for Example 2 is shown in Fig. 2. Unartain parameters, constraints, exchanger data and heat capacities are listed, respectively, in Tables 4-7. The “bottleneck” occurs in 332,with temperatures T3, TS and T6 at their limits. T7 is the only communication between the uncertain parameters in the lower part of the network (aside from F3) and the active constraints, and there are enough control degrees of freedom to compensate for their variations. Of the 24I uncertain parameters, the six fouling resistances do not give rise to branching with either method. However, the original method develops branches on each of the other 14 parameters. Bounding restricts the total number of nodes to 1067. This problem has nine compensatable parameters (F4, FS, F6, F7, T9, TlO, T13, T16, T18). With the mod&d method this set is properly identified at the in etfect throughout iirstnodeinthetree,andremains the tree. The problem is solved after 31 nodes. Tabk 6. lhmpk k
6: QW m-‘K-l)
The network for Example 3 is shown in Fig. 3. Uncertain parametem, constraints, exchanger data and heat capacities are listed in Tables 8-11, mspectively. The “bottleneck” appears in E3, with temperatures T8 and Tll at their limits. This example illustrates a case wherein the set of compensated parameters changes during tree development. There are two active sets that appear at various nodes in the tree: one in which T14 is at its limit, and one in which T16 is at its limit. At the top node, the projected compensatable parameter set is (Fl, F2, F6, Tl, T4, T6, T12, TlS). The second active set is encountered at the eleventh node evaluation, with a new, smaller projected parameter set that no longer contains (F6,T12,TlS). The active set changes at several Table 5.
i
&le 2 constraints Lower limit
Variable
2 exchanger data
c (kg 8-l)
fi: (kW m-‘K-l)
b
e
I-‘)
(m’K?W-I)
upper limit
1070
U. KAIWIBK and R.
E. SW=
Table 7. Exampk 2 heat atpacltiea
I
C’ (k.l kg’! K- ‘)
(K)
::I: 1.0 4.0 1.9 2.0 2.5
273 273 293 298 273 273 273
: 3 4
7
Table 8. Ihmpk
E’ (kl k$K-‘1
f,j
2)
1.0 2.8
373 373 373 473 373 373 373
:3 2.0 2.0 2.3
K
0.5
+0.5
-0.5
: 4 5 6 7 8
I, r5(m pi (m’K (rn: KkW-) kW-‘) kW-:) rt im’K kW-‘j I, (m’K kW-‘) 4 (Lg p-‘) F2 (kg 8-l) F, OrIx-‘) F4 Or.ss-‘) & ha-‘) Fs Orgs-‘) T’ (K) T4 (K) Ts (K) T. (K)
0.5 0.5 0.5 12 10 10 14 8
+0.5 +os +0.5 +0.5 +1.0 +1.0 +1.0 +1.0 +1.0 +1.0 +5.0 f5.0 +5.0 +5.0 +5.0 +5.0
-0.5 -0.5 -0.5 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0
1 2 3 4 5 f
;:
I
F5
T12
F4 Fig. 3.
Example 3 network.
places in the tree, and at one node the anti-cycling Sag is invoked. The node involved gets pruned after a singIe branch. Other than the five fouling resistances which never require branching, the original method develops branches on each of the remaining 12 parameters. Rounding restricts the total number of nodes to 175. The modi&d method required only 45 nodes to solve this problem. Table 12 summarims the problem sizes and the relative performance of the original and modified
: 4 5
100 100 100
1 t 1
t: 20 20
$
x
ii
:
(R=@)-' (R't?')-' (R'@)-' (R'@-'
C (W
3900
545
Ts (K) G (K) T” (K) TM (K) TM (K)
0 450 400 0
4% al 4%
CONCLUSIONS An improved branch and bound strategy has been presented for the primal method of worst-case identification. The modified method permits a reduction in
data
i
&%K-1) (kg:-‘) O;wm:‘K-‘)
k
limit
method in terms of node counts. Also shown in the table is the average number of equivalent single-state NLPs per node evaluation. This number reSects the need for distinct values of 6**’ in the node NLPs. Since a separate set of state variables must be determined for each distinct &‘, this number may be used when comparing the computational requirements of the node NLPs to the work required to solve a basic model optimization in a single state. For the modified method, this average number includes the requirements for the few repeated node evaluations and the single state solutions performed to project the parameter sets. NLP solutions were computed using a successive quadratic programming strategy starting from the feasible nominal point (6 = 0).
Table 10. Examok 3 achauer
.
Uppa (R’u’)-’
0
8 9
F6
3 oonstmints Lower limit
Variable
I
ae:
z 450 530 330 470 370
Table 9. Examde
T10
-A@,-
1
t: 15 16 17
F2
paramercrs
e/”
@J
1: 11 12
Fl
3 unartain
i
1 I 1 1
4
($-‘)
(m* K kW-‘)
10 z 20
0 0
8
Worst-case
identification
in stmchued
1071
process sys-
Tabk 11. Example 3 heat capacities
a I
(W kg-l:
I 2 3 4 :
fi
K-‘)
(K)
2.0 2.0 2.0 1.0 1.0 1.0
o;l$K-‘)
280 280 280 280 280 280
g\
2.0 2.0 2.0 1.0 1.0 1.0
380 380 380 :: 380
Table 12.Summaryof results Node evaluations
F!quivakot NLPS
Variables
:
20 14 6 5 46 62
32,767 511
1067 57
31 9
1697 119
3
17 5 52
8191
175
45
365
the number of branches in the search tree by recognizing the presence of “compensatable parameters” within the problem. For the examples studied, the modified method was very effective at identifying the compensatable parameters and avoiding unnecessary branching. The branch and bound effort is reduced accordingly; in some problems a reduction of more than an order of magnitude may be realized. Since compcnsatable variations are commonly encountered, the method’s advantages should be of value in many typical problems.
Akzl*zk~k-Qk=O,
x:iik-
@AAD).
Floudas C. A. and I. E. Grossmann, Active constraint strategy for flexibility analysis in chemical processes. Computers them. hgineering 11, 675-693 (1987).
Pistikopoulos E. N. and I. E. Grossmann, Optimal retrofit design for improving process flexibility in nonlinear eystem~1. Fixed degree of flexibility. Computers them. Bgineering 13, 1003-1016 (1989). Swaney R. E. and I. E. Grossmann, An index for operational ilexibility in chemical process design. Part I: formulation and theory. AXhE 3I 31, 621-530 (1985a). Swaney R. E. and I. E. Grossmann, An index for opcrational flexibility in chemical process design. Part II: computationalalgorithms. AIChE Jl31,631-641(1985b).
(A4) (A5) (A6)
3(x:, x:, = 0,
where
x:, x: > 0, Otherwise
and
x:e-4 - x:e-4 x:-x: ’
I
Z,I_’ e-q-4
’ Q=m-_
r:(Ff)-“” + r:(Ff””
REFERENCES
(A3)
T&$+ T&=0,
x:a;L-T*,+Pti=O,
&= Acknowledgement-U. Kabatek gratefully acknowledges partial support by the German Academic Exchange !3ervice
169
+ r; + rt - R* = 0,
(A7) (A@
zk*O, Rf’(7’
1
- 1 Q 0.
(A9)
Equations for splitter k T:,-T;=O,
(AlO)
TL-Tt=O, F$, - z:e
(All) = 0,
F”, - (1 - z:,zq = 0,
o
I.
(Al2) (A13) (A14)
Equattons spec~ying constants APPENDIX
(Al5)
Modeling Equations for Heat Exchanger k Fj$f
+ b;(T:,
+ Th)/2](T:, - Th) + Qk = 0,
Fk[nr ,-i-i)-Qk=O, 5c + b’(Tk EOD + T”)/2J(T* n
(Al) W)
(AW
Compufers hem. &gng, Vol. 16, No. 12, pp. 1063-1071,1992 Printed in Great Britain. AU rights mservcd
WORST-CASE
IDENTIFICATION IN STRUCTURED PROCESS SYSTEMS U. I(ABATBKt and R. E. SWANEY$
Department
of Chemical Engineering, University of Wisconsin, Madison, Wl 53706, U.S.A.
(Received 20 September 199o;jnal reuiGm received 14 May 1992; receivedfor publication 28 May 1992) A--An improved method for identifying the worat-caae set of parameter combinations in controllable systems is presented that is based on branch and bound search in the space of the uncertain parameters. A characteristic of controllable structured systems frequently appears wherein the limiting condition is local&d in a part of the system that has only limited coupling to the other parts. While there may be many uncertain parameters varying in parta upstream of the limiting constraints, often there will remain suf&ient control dcgoxs of freedom in those other parts to compensate. The branch and bound method in ita simplest form cannot recognize this situation, and will usually engage in searching among combinations of such parameters, when in the end they have no effect on the limiting condition. Modifications to the basic method are presented that allow these parameters to be nxagnizcd apriori. The complexity of the branch and bound search is reduced accordingly, while the rigor of the original method is maintained. As shown in three heat exchanger network examples, a substantial reduction in effort is often achieved.
INTRODUCTION
It is a reality of practice that chemical processing systems cannot simply be designed for some fixed nominal set of conditions, but must be designed to successfully accommodate a range of conditions. The origins of these variations fall into two categories. In the first are expected variations in operating conditions such as feed and product specifications, catalyst activity, utility supplies and equipment fouling. In the second are the uncertainties inherent in the design models themselves and in the physical data employed. While the former are known to vary, and the latter do not actually vary but are unknown, both types of uncertainty impose requirements for the design to meet. The task of design under uncertainty is to ensure that the worst case within the range of expected variations (or realizations) can be successfully handled. This requires that the worst-case combination be identitied, and then designed for. Worst-case identification in a complex process system often is not a trivial task. The existence of control adjustments that can be made in response to variations contributes to the complexity of the analysis. A formal framework for addressing these concerns was presented by Swaney and Grossmann (1985a) in the form of a “flexibility index”. In that framework,
the problem variables are classiiied as: design variables d, which am selected during design but are fixed during operation; control variables E, which may be varied to suit circumstances during operation; uncertain parameters 0, which exhibit a range of values during operation that must be accommodated, and state variables x, which are dependent variables in the process model. The ranges of variation in the uncertain parameters %,are described by deviations A@;, Aii,? below and above a nominal operating point 07. The index value F measures the fraction of the target ranges of variation that the design can accommodate for simultaneous independent variation of the parameters. Process operating constraints and product speci&tions are specifM by the set of inqualities g(d, z, 0, x) C 0 while the process model is represented by the set of equations h(d, 2, 8, X) = 0. The latter are suflicient to determine the state variables x. While generally the variables x and the relations h( . , . , *, x) = 0 are r&ah& explicitly in implementation, for discussion it is conveniat to treat the state variables as implicit functions x = h -‘(d, z, 0) and to eliminated them. The constraint functions may then be expressed in the reduced form g(d, z, 0, x) =f(z, 6), with the dependence on the design variables implied but suppressed. The analysis problem may be formulated as:
(1)
s.t. .f(z, e) < 0,
pmellt addreux lnstitut mrsystcmdynamik und r&geltecbnik, ulkniitit stmtgart, PfafFenwaldkg Stuttgart go, Germany. ~To whom all correspondence should be w.
F=mjnyS
9. D7ooo
8 =P+88,
(2) (3) (4)
1063
1064
U. KABATEKand R. E. SWANEY
to represent whether or not a constraint is active. For a given combination of these, the maximization part of (1) becomes implicit, so the problem becomes one of minimizing over 6 and the binary variables. The solution is then found either by solving the resulting mixed integer program (an MILP for the linear case), or by enumerating the feasible combinations of active constraints and solving the NLP in 8 for each combination. Since this method searches among L=max 6, (5) the constraints, it may be described as a “dual” d* method. For further applications of the active cons.t. f;(z, el) 6 0 (6) straint approach, see Pistikopoulos and Grossmann i=l m, (1989). @=@N+@ ‘*“’ (7) > Regarding relative merits, the two approaches offer somewhat different capabilities. Although both By evaluating each constraint at its individual worstprocedures are rigorous and effective for purely linear case, 6’= 8**‘, the solution to this NLP will provide models, the active set method usually requires fewer a lower bound to (1). binary decisions and fewer underlying combinations. The basis of the implicit enumeration algorithm is The dual approach is thus well-suited for linear to employ a branch and bound procedure to search problems, wherein its strengths can be fully exploited. the domain ee(-AO;,AOT}, j= 1,. . . , no, to For complex nonlinear problems, the primal approach locate the worst-case &*. At a node in the search presents its own advantages. First, the primal method tree, each parameters 4 is either “unassigned” or is easily implemented for complex system models. “assigned”. Unassigned parameters j e J take on the by a small set of values in the individual constraints given by 8; = &;CJ. The subproblems are representable NLPs loosely coupled in a multipcriod manner, each Assigned parameters j E J, on the other hand, are of which is a parametric optimization of the confixed either at 8,= -AO; or 8/ = AO:t and take ventional process model. It is therefore not difficult on the same value in every constraint. The tree is to adapt existing flowsheet optimizers to this task. developed by taking a node with the lowest lower Equivalent simplicity of implementation is obtained bound, selecting one of its unassigned parameters to for the dual approach only if the method of enumerbecome assigned, and replacing the node by separate ating active sets is employed. Moreover, a feasible descendant nodes, one for each of the possible starting point is available for every NLP calculation assignment values. This amounts to partitioning the required by the primal method. Secondly, since the search domain according to the value of the assigned resulting solution is based on a series of primal-feasiile parameter and computing individual bounds for calculations, any impact of general nonlinearities on the partitions. The validity of the collective bound is the global validity of the result is likely to be conmaintained, since for every assigned parameter, each servative, in the form of suboptimal control settings. extreme value is represented by a node in the tree. Global rigor rests on the ability to locate parameter The procedure is continued until enough parameters values that are worst-case for each constraint indihave been assigned to close the gap between the lower vidually, and more intensive methods may be applied bound and the solution. The worst-case parameter to this subtask if needed. The impact of nonlinearities combination 8* is then readily identied within in the dual formulation is less clear. Finally, it should the partition having the lowest bound. Since this be noted that both approaches usually are capable procedure executes the worst-case search directly in of locating nonvertex solutions (i.e. 8: interior to the variables 8, it may be. described as a “primal” [ - AO; , A@; ] for some j) when they occur in typical method. problems. In the primal method these are usually The second approach is the “active constraint detectable and accessible by descent from a vertex strategy” outlined by Floudas and Grossmann using a local minimization. (1987). The idea is to search among the possible The purpose of the present paper is to present an combinations of constraints that will be active at improved branch and bound strategy for the primal the worst-case condition. Binary variables are defined method. The modified method permits a reduction in the number of branches in the search tree by recog?A generalization of this exists for solution points that arc nizing the presence of “compensatable parameters” not at the parameter bounds; a branch mightconceivably within the problem. The resulting reduction in solugenerate an assigned value 8, interior to [ - AO;-, AfI,+1, possibly with 4 = 4(z, 0). tion effort can be substantial in typical problems. The solution O* = ON + 6*8* represents the worstcase operation at which limiting conditions are reached. Two main approaches have been proposed to solve (1). The first is the implicit enumeration method presented by Swaney and Grossmann (1985b). This method is based on a lower bound program of the following form:
1065
Worst-ease identiticatiou in structured process systems F5 T14 F3
T15
I T7
I-T10 T5
+ Tll F4
F,n+$
T2
($3
I
T9
i Fig.
1.
T13
Example I network.
COMPENSATABLE
PMMGXERS
This term is introduced to describe a behavior frequently exhibited by structured process models. O&en there are a number of parameters whose values have no impact on the feasibility of the worst-case limiting condition. This situation is principally a commquence of system spar&y. Limiting conditions usually arc encountered in only a part of the system (e.g. where the “bottleneck” is located), while other parts operate inside of their limits and still have degrees of freedom for compensation. An illustration of this behavior is provided by Example 1, shown in Fig. 1. The data describing this small heat exchanger network am given in Tables l-3. Heat transfer coe!Bcients vary with flowrate and
Table 1.
Exampla1 mwrtain
e,
i :
e/”
r{ (m% kW-‘) rr (m’K kW-‘1
r) &K kW-'j
0.5 0.5
r: (mzK kW_‘) F, &t~a-‘) Fa (ks a-‘) JS &I-:) ? (s)
ix 10’ ii 5:
T: (IO Tl(K) T,, IKl
723 388 313
parametera -AO;
ae; +0.5 +0.5
+0.5
+0.5 +1.0 +2.0 +2.0 +3.0 +5.0 +5.0 +5.0 +s.o
-0.5 -0.5 -0.5 -0.5 -1.0 -2.0 -2.0 -8.0 -5.0 -5.0 -5.0 -5.0
Table 3. ii’ (kWm”K-I)
k
constant heat capacities of I kJ kg-’ K-’ apply for every stream; details are given in the Examples section. In this example, flowratc and inlet temperatures for streams Fl-F4 are treated as uncertain parameters, along with the fouling resistancea in each exchanger. Each exchanger has a variable bypass; together with the split fraction controlling the flow into E4 they comprise the five control variables. Outlet temperatures must meet inequality con&mints. The limiting condition in this example occurs in the lower part of the network around El and E3, with El operating at zero bypass. Of the 12 uncertain parameters, only Fl, F3, F4, Tl, Tl2 and the fouling in El and E3 impact the limiting condition. The original branch and bound procedure generates 57 nodes in solving this problem, with branches on six of the variables. That method is inherently able to avoid branching on the four fouling resistances or on F4 or T12, since in this problem they have uniform worst-case values in each of the limiting constraints. However, branches are made on the uncertain parameters F2, T4 and T7. Variations in these parameters in fact do not effect the result. This is a consequence of the structure of the network and of the fact that control degrees of freedom remain in the upper part of the network. Temperature TS is the only connection between variations in these parameters and the limiting constraints, The upper part of the network is able to adjust T8 to the optimum value required by the constrained lower part, in spite of variations in F2, T4 and T7. These parameters are thus termed compensatable. In the original method, branching paranmters are selected by computing the sensitivity of the node lower bound program to their assignment. Although the solution node is insensitive to changes in a compensatable parameter, intermediate parent nodes
Example1 cxdunspr e (kgr-‘)
data
fi: (%Wm-*K-l)
e (kgs-‘)
(m’K$W-‘)
:
ii
1.5
20 10
1.5
iti
0
3 4
120 12
1.5 1.5
ii
1.5 6.0
30 10
: 0
U. lCm~mx
1066
and R. E.SWANEY
may respond significantly to aaaignment of the parameter. This behavior results when the individual constraint worst-case values 87.’ are not uniform for the parameter. The effect of the compensatable parameters on the node lower bounds does not disappear until their assignment makes them uniform in the active constraints. The presence of such parameters artificially complicates the worst-case identification task for the original method. Through proper recognition of compensatability, however, the added complexity can be essentially eliminated. In the case of Example 1, the number of nodes in the branch and bound tree reduces from 57 down to 9. MODIFIED
METHOD
The modified method seeks to avoid unnecessary branching by recognixing the presence of compensatable parameters a priori. Since there is no general a priori way to identify such parameters with certainty, the method works with a projected (postulated) set of compensatable parameters, which is revised as needed during the course of the algorithm. The projected set is then either confirmed or refuted when branching appears complete. This approach permits efficient solution with nearly a minimum of branching, yet retains the rigor of the earlier method. In the following, frequent use is made of linearized perturbation analysis of the solution of the node lower bound program (5). Of particular use are the predicted changes in the bound accompanying changes in the vertex values C. Let 12 0 be the vector of multipliers for the constraints (6) at the solution of the lower bound program (5) for the node. For the changes in parameter 4 from the values 8: at the no& to the uniform assignments A6: and -A@;, the respective bound change predictions are given by:
The projected set of compensatable parameters is generated as follows. At a node k early in the tree, parameters Rj, j lJL are assigned to fixed vertex values g):, while the unassigned parameters 0;, j E Jk taken on the various individual worst-case values @’ for each constraint i. Since parameter compensatability is not usually evident until all intluential parameters are assigned values that are uniform for every constraint, a complete assignment is selected. A likely candidate for the worst-case vertex within the node’s subtree may be predicted by selecting
assignments that give the best predicted bound increase according to (8,9). A complete assigmnent is thus obtained by setting the unassigned parameters according to:
y=fJ;,
jeJL
(10)
(5) is then solved for CT’= 6k for this projected vertex. At this vertex solution, parameter sensitivities are again computed from (8,9). Compensatability is readily evident, and the projected set of compensatable parameters C” for node k is taken as:
C*={jly,Z=y;=O,jeJk}.
(11)
The compensatability
property of a parameter derives from and depends on the set of constraints that are active at the worst-case vertices. The active set associated with the projected set C’ is also indicated at the vertex solution by {iI&> 0}, and will be denoted A(Ck). With a projected set identified, the branch and bound procedure continues as in the original method, except that parameters in C* are reserved from selection as branching parameters. Branching parameters are selected from the complement set Cc. Each new subnode k’ generated originally inherits the set Cc’ = C* of its parent. However, the eventual applicability of a set C” to vertices below node k’ is likely to be contingent upon having the same active set that generated Cy also be active at those vertices. For this reason, if the active set I$’ exhibited at the solution to the node’s lower bound program differs from A(C’), a new set of compensatable parameters is projected for node k’ using the procedure given above. C’ is then taken as the new projected set and branching continues. Note that the use of projected sets CK in these steps implies no assumptions, but merely serves to guide the branching strategy. When a node z is produced by the last branch indicated within the set Cc of its parent node k (i.e. no further branching would be indicated in c), the compensatable parameter set Cr = Ck is usserred. In this step, the parameters in C” are assigned to uniform values 8, in the active constraints i E A (CE). while retaining their individual worst-case values 8:’ in the inactive constraints. This action seeks to merge the combinations of the parameters in Cr into a single node. The bound computation (5) is performed for E with Cr asserted. At the solution y,+ , y,- are determined from (8,9). The asserted set is then con8rmed if: yj’=y;=O,
VjeC4
(12)
Worst-ease identificationin structuredproasa systems If the asserted set is continned, the bound computation for z represents a valid lower bound, and no further branching at the node is required. If the set is refuted, however, the projected set is abandoned. A new projected set is determined for the node as above. The bound computation is recomputed with this second set asserted. If the second set is also refuted, the lower bound is again recomputed with the set. not asserted, and node E remains open for further branching. Confirmation criterion The purpose of the confirmation test is to verify that variations of the compensatable parameters within the limits -AB; < 4 6 AtI,? do not a&ct the value of the lower bound computation (5). SuIIiciency of the criterion (12) derives from the following assumption: Regularity condition-At the solution to the lower bound (5) for node E with C” asserted, let Rz, 0) be the vector of ti active constraint functions [i.e. x(z, 0) =&(z, 0) for i = 1, . . . , ff; i' EZ:],with corresponding multipliers 1, > 0. Then the regularity condition holds if:
rank
(13)
over the feasible values of 8 and z. Stationarity conditions at the solution require that:
to (5)
which means the rank of g/azr cannot exceed R - 1, The regularity condition requires the constraint gradients to have maximum independence in the control variables within this limit. The condition will hold at the solution point except in cases with nonunique multipliers, which could arise either from redundancy in the constraints or unlikely coincidence. A simple remedy is to stipulate a solution with the minimal number of positive multipliers, with A mdurxd accordingly. The extension to neighboring values of 0 and z is a modest condition that will nearly always hold in practice. Under the regularity assumption, the null space of the active constraint gradients in the controls is of rank one, spanned by x(z, 0) that is defined implicitly by x’(aJlaz’) = 0 over the feasible values of z and 8. Consider alternate values P(z), s e [O, 11, where B’(0) = 6r is the value at the solution to (5). Under the requirement that: Ix(z, t9]f$(z,
8.)]%
= 0,
(15)
1067
over e’(s) and z(s), the imphcit function theorem provides the existence of z(z) such that nz(z), e'(z)]= 0. Variations in the compensatabk parameters correspond to various 0; (s), j E C’, while holding 0;) je CI, constant. The existence of the condition r; =r; = 0 in (12) implies: XTg=O,
VjJj@,
I
at the solution to (5), and so provides that condition (15) will hold at that point. There is nothing inherent in the formulation of (5) itself to induce condition (16). Thus, excepting extreme coincidence, the appearance of condition (12) indicates a structural characteristic of the system. This characteristic exhibited byxz, 0) will not disappear with variations in e,, j E d, so requirement (15) may be expected to hold over these variations. The net result of this argument is that the control variables will possess sufficient degrees of freedom to maintain: l;(z,e) = 0, ie z: (17) for the active constraints over variations in the compensatable parameters. No change is implied in the objective function variable S. The remain@ requirement for validity of the lower bound with the compensatable parameters asserted is to ensure feasibility of the inactive constraints. This is immediately provided, since they are evaluated at their individual worst-case values @’ for / E Cxc Jr. An asserted parameter set may thus be confirmed or refuted by checking condition (12).
Detine
Jx = index set of unassigned parameters at node k, Jk = index set of assigned parameters at node k, K = index set of open nodes, Bk = index set of branch candidate parameters, #Bk=cardinality of set Bk. Branch and bound 0. set
P=0, co=0, A(CD)=0, +&.ceat,,==nb, K = (0). Evaluate node 0.
1. Detelmine A ==argr${L,}.
2. Branch on parameter I = i(k). (a) if #B,t-0, Stop; F=t,.
U. KAMTEK and R. E. SWANEY
1068
(b) form subnodes
k’ and k”:
6:=85=8;, &‘=
1. Project vertex 9, j e Jk, using (10) with y+. yfrom node evaluation. 2. Solve (5) with 8’= dk and obtain active set A(CI). 3. Compute y:, y; at projected vertex using (8,9). 4. Determine Ck from (11).
je:JL,
-A@,-,
&=A0:,
~=‘=JQ{Z}, A (CK) = A(P)
Parameter set projection @ode k)
CK=CY=Ck, = R (C”).
3. Evaluate node k’. 4. Evaluate node k”. 5. Replace node k by nodes k’ and k” in K, and go to step 1.
Node evaluation (node k’, parent k) 0. If # Bh = 1 perform Step 1. Otherwise perform Step 2 with ilag not set. 1. (Assert C”): (a) select #yet-Ae-,Ae+]
forjoCY;
for (j d P) or [j E Ck’, i E A (F)], @*’ for (j e Jv n c’) or [j e CK, i 4 A(P)];
VT (b) set fl=
(c) solve (5) to obtain LK; (d) compute Y_T , yl from (8,9) for j e JK; = 0 for j E CK and min{Y,+, y;} = 0 (e) ifY/‘=Y; for j E JK, them C’ is continned and LK is valid; node evaluation is complete; (f) otherwise, project a new parameter set and repeat Steps (a-e) with the revised Cy; (8) if the revised Cy is not confirmed, perform Step (2) with flag set. 2. (CK not asserted):
solve (5) -to obtain & and Z$; if A (Ck’) $ Zs,project a new parameter set and revise c”; compute y,+, y,- from (8,9) for j E J”; determine branch candidates within set M = JKnE’:
W3 = am~IY,}, ?I=
-4 {
(18)
IY,+ r,-l, a) -
&={j]Y~)O,
ifmin{v;,7;) r otherwise,
0
jeM}
0, (19) (20)
(f) if # BK > 0, node evaluation is complete; (g) otherwise, repeat (e) for M = Jy; (h) if redetermined #BE = 0 or flag is set, node evaluation is complete; (i) otherwise go to Step 1.
In Step 2a of the branch and bound procedure, one may wish to project to a vertex to obtain a complete assignment, and then evaluate the vertex in order to confirm the solution. This will deal with the rare chance in which a gap remains when the branch selection indicates no candidates. Further branching would be performed if required. At the solution vertex it is of value to test for local descent in the function 6 *(0) to check for a nonvertex solution accessible by local descent from the vertex. The test is performed by checking for negative values of A’tTf/atl,(-8,) at the vertex solution. A local minimization of S*(B) may then be performed if so indicated. The simplest way to obtain individual constraint worst-case values 8.’ is to assume monotonicity in each parameter and to select - A0,- or Af?: based on the signs of the gradients @J&9,. In this case y,+, y,should be monitored for the appearance of negative values, which would indicate a failure of the monotonicity assumption. In response, duplicate constraints with multiple 8**’at both extremes can be employed. Alternatively, a local search may be performed to locate a possible nonvertex &I*‘, with duplicate constraints then employed as indicated. EXAMPLES
Comparative results for the mod&d method are presented here for three heat exchanger network examples. In these examples, heat transfer coefficients are functions of flowrate, and enthalpies are quadratic functions at temperature. In each problem stream flowrates and inlet temperatures, and exchanger fouling, are uncertain parameters. Each exchanger has a variable bypass, and stream outlet temperatures are specified by inequalities. For this formulation the reduced constraints are monotonic in the uncertain parameters, giving certain determination of the P. The modeling equations may be found in the Appendix. For Example 1 described earlier, the modified method requires nine nodes, reduced from 57 with the original procedure. The set of compensatable parameters is correctly identified by projecting at tbe top node. The active set remains unchanged throughout the tree, so the projected set remains in use throughout.
identihtion in strucked processsystnns
Worst-case
1069
Tabk~4.Eumpk12unsrtlin~ e,
/ 1
2 3 4 5 6
F7
F4~l~& (&
t: 13 14 15 16 17
Tll
Jw
-A$
1 t::
+0.75 +1.0 +0.9 +1.0 +1.0 +0.7 +1.0 +1.0
-0.75 -1.0 -0.9 -1.0 -1.0 -0.7 -1.0 -1.0
ii 30 30 50 553 373 323 423 373 333 353
+1.0 +I.0 +1.0 + 1.0 +5.0 +5.0 + 5.0 + 5.0 f5.0 fS.0 + 5.0
-1.0 -1.0 -1.0 -1.0 -5.0 - 5.0 -5.0 -5.0 -5.0 -5.0 -5.0
0.75 1.0 0.9
r{ r5(m (m;K KkW-) kW-:) I, (m* K kW-‘)
ri (III*K kW-1.
5 9 10
I
e/w
T19
I (m*KkW. -I rrI (m’K kW, -‘) F, Or6s-‘) 4 Orrs-‘) F2 (Ltl a-‘) F. (Ls 8-l) 4 (kg s-v &ckes-:)
9 FrT: (K) T, (K) T,ci(K) T,, (K)
0.7 m 8
)
Tn. (K) T;; *)
20
Example 3
F3
Fig. 2. Example 2 network. Exampl% 2 The network for Example 2 is shown in Fig. 2. Unartain parameters, constraints, exchanger data and heat capacities are listed, respectively, in Tables 4-7. The “bottleneck” occurs in 332,with temperatures T3, TS and T6 at their limits. T7 is the only communication between the uncertain parameters in the lower part of the network (aside from F3) and the active constraints, and there are enough control degrees of freedom to compensate for their variations. Of the 24I uncertain parameters, the six fouling resistances do not give rise to branching with either method. However, the original method develops branches on each of the other 14 parameters. Bounding restricts the total number of nodes to 1067. This problem has nine compensatable parameters (F4, FS, F6, F7, T9, TlO, T13, T16, T18). With the mod&d method this set is properly identified at the in etfect throughout iirstnodeinthetree,andremains the tree. The problem is solved after 31 nodes. Tabk 6. lhmpk k
6: QW m-‘K-l)
The network for Example 3 is shown in Fig. 3. Uncertain parametem, constraints, exchanger data and heat capacities are listed in Tables 8-11, mspectively. The “bottleneck” appears in E3, with temperatures T8 and Tll at their limits. This example illustrates a case wherein the set of compensated parameters changes during tree development. There are two active sets that appear at various nodes in the tree: one in which T14 is at its limit, and one in which T16 is at its limit. At the top node, the projected compensatable parameter set is (Fl, F2, F6, Tl, T4, T6, T12, TlS). The second active set is encountered at the eleventh node evaluation, with a new, smaller projected parameter set that no longer contains (F6,T12,TlS). The active set changes at several Table 5.
i
&le 2 constraints Lower limit
Variable
2 exchanger data
c (kg 8-l)
fi: (kW m-‘K-l)
b
e
I-‘)
(m’K?W-I)
upper limit
1070
U. KAIWIBK and R.
E. SW=
Table 7. Exampk 2 heat atpacltiea
I
C’ (k.l kg’! K- ‘)
(K)
::I: 1.0 4.0 1.9 2.0 2.5
273 273 293 298 273 273 273
: 3 4
7
Table 8. Ihmpk
E’ (kl k$K-‘1
f,j
2)
1.0 2.8
373 373 373 473 373 373 373
:3 2.0 2.0 2.3
K
0.5
+0.5
-0.5
: 4 5 6 7 8
I, r5(m pi (m’K (rn: KkW-) kW-‘) kW-:) rt im’K kW-‘j I, (m’K kW-‘) 4 (Lg p-‘) F2 (kg 8-l) F, OrIx-‘) F4 Or.ss-‘) & ha-‘) Fs Orgs-‘) T’ (K) T4 (K) Ts (K) T. (K)
0.5 0.5 0.5 12 10 10 14 8
+0.5 +os +0.5 +0.5 +1.0 +1.0 +1.0 +1.0 +1.0 +1.0 +5.0 f5.0 +5.0 +5.0 +5.0 +5.0
-0.5 -0.5 -0.5 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0
1 2 3 4 5 f
;:
I
F5
T12
F4 Fig. 3.
Example 3 network.
places in the tree, and at one node the anti-cycling Sag is invoked. The node involved gets pruned after a singIe branch. Other than the five fouling resistances which never require branching, the original method develops branches on each of the remaining 12 parameters. Rounding restricts the total number of nodes to 175. The modi&d method required only 45 nodes to solve this problem. Table 12 summarims the problem sizes and the relative performance of the original and modified
: 4 5
100 100 100
1 t 1
t: 20 20
$
x
ii
:
(R=@)-' (R't?')-' (R'@)-' (R'@-'
C (W
3900
545
Ts (K) G (K) T” (K) TM (K) TM (K)
0 450 400 0
4% al 4%
CONCLUSIONS An improved branch and bound strategy has been presented for the primal method of worst-case identification. The modified method permits a reduction in
data
i
&%K-1) (kg:-‘) O;wm:‘K-‘)
k
limit
method in terms of node counts. Also shown in the table is the average number of equivalent single-state NLPs per node evaluation. This number reSects the need for distinct values of 6**’ in the node NLPs. Since a separate set of state variables must be determined for each distinct &‘, this number may be used when comparing the computational requirements of the node NLPs to the work required to solve a basic model optimization in a single state. For the modified method, this average number includes the requirements for the few repeated node evaluations and the single state solutions performed to project the parameter sets. NLP solutions were computed using a successive quadratic programming strategy starting from the feasible nominal point (6 = 0).
Table 10. Examok 3 achauer
.
Uppa (R’u’)-’
0
8 9
F6
3 oonstmints Lower limit
Variable
I
ae:
z 450 530 330 470 370
Table 9. Examde
T10
-A@,-
1
t: 15 16 17
F2
paramercrs
e/”
@J
1: 11 12
Fl
3 unartain
i
1 I 1 1
4
($-‘)
(m* K kW-‘)
10 z 20
0 0
8
Worst-case
identification
in stmchued
1071
process sys-
Tabk 11. Example 3 heat capacities
a I
(W kg-l:
I 2 3 4 :
fi
K-‘)
(K)
2.0 2.0 2.0 1.0 1.0 1.0
o;l$K-‘)
280 280 280 280 280 280
g\
2.0 2.0 2.0 1.0 1.0 1.0
380 380 380 :: 380
Table 12.Summaryof results Node evaluations
F!quivakot NLPS
Variables
:
20 14 6 5 46 62
32,767 511
1067 57
31 9
1697 119
3
17 5 52
8191
175
45
365
the number of branches in the search tree by recognizing the presence of “compensatable parameters” within the problem. For the examples studied, the modified method was very effective at identifying the compensatable parameters and avoiding unnecessary branching. The branch and bound effort is reduced accordingly; in some problems a reduction of more than an order of magnitude may be realized. Since compcnsatable variations are commonly encountered, the method’s advantages should be of value in many typical problems.
Akzl*zk~k-Qk=O,
x:iik-
@AAD).
Floudas C. A. and I. E. Grossmann, Active constraint strategy for flexibility analysis in chemical processes. Computers them. hgineering 11, 675-693 (1987).
Pistikopoulos E. N. and I. E. Grossmann, Optimal retrofit design for improving process flexibility in nonlinear eystem~1. Fixed degree of flexibility. Computers them. Bgineering 13, 1003-1016 (1989). Swaney R. E. and I. E. Grossmann, An index for operational ilexibility in chemical process design. Part I: formulation and theory. AXhE 3I 31, 621-530 (1985a). Swaney R. E. and I. E. Grossmann, An index for opcrational flexibility in chemical process design. Part II: computationalalgorithms. AIChE Jl31,631-641(1985b).
(A4) (A5) (A6)
3(x:, x:, = 0,
where
x:, x: > 0, Otherwise
and
x:e-4 - x:e-4 x:-x: ’
I
Z,I_’ e-q-4
’ Q=m-_
r:(Ff)-“” + r:(Ff””
REFERENCES
(A3)
T&$+ T&=0,
x:a;L-T*,+Pti=O,
&= Acknowledgement-U. Kabatek gratefully acknowledges partial support by the German Academic Exchange !3ervice
169
+ r; + rt - R* = 0,
(A7) (A@
zk*O, Rf’(7’
1
- 1 Q 0.
(A9)
Equations for splitter k T:,-T;=O,
(AlO)
TL-Tt=O, F$, - z:e
(All) = 0,
F”, - (1 - z:,zq = 0,
o
I.
(Al2) (A13) (A14)
Equattons spec~ying constants APPENDIX
(Al5)
Modeling Equations for Heat Exchanger k Fj$f
+ b;(T:,
+ Th)/2](T:, - Th) + Qk = 0,
Fk[nr ,-i-i)-Qk=O, 5c + b’(Tk EOD + T”)/2J(T* n
(Al) W)
(AW