A Multiparametric Programming Approach for ... - ACS Publications

Mar 3, 1997 - Reactive Scheduling by a Multiparametric Programming Rolling Horizon Framework: A Case of a Network of Combined Heat and Power Units...
0 downloads 0 Views 228KB Size
Download PDF
Ind. Eng. Chem. Res. 1997, 36, 717-728

717

A Multiparametric Programming Approach for Linear Process Engineering Problems under Uncertainty Joaquı´n Acevedo† and Efstratios N. Pistikopoulos* Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College, London SW7 2BY, U.K.

In this paper, a parametric programming approach is proposed for the analysis of linear process engineering problems under uncertainty. A novel branch and bound algorithm is presented for the solution of mixed-integer linear programming problems where n right-hand-side parameters are allowed to vary independently. The procedure, based on the solution of multiparametric linear programs at each node of the tree search and special bounding procedures, identifies the different optimal integer solutions and their corresponding optimal value functions as the uncertain parameters vary within given ranges. Three examples are presented to illustrate the basic steps of the algorithm and its applicability to process engineering problems. Introduction Applications of mathematical programming approaches for the solution of process engineering problems have increased considerably in the last 2 decades. In particular, mixed-integer formulations have been used extensively to take into account the discrete (e.g., existence of process units in a synthesis problem) and continuous decisions (e.g., operating conditions such as flowrates and temperatures) involved in these problems. In many cases, these formulations can be posed as mixed-integer linear programming (MILP) problems as follows:

min cx + dy y,x

s.t.

(1)

Ax + Ey e b xg0 y ∈ {0, 1}

where y is the vector of 0-1 integer variables and x is the vector of continuous variables. All matrices and vectors (A, E, b, c, d) are constant with conforming dimensions. The objective function is typically given by an economic criterion, while the mixed-integer constraints represent mass and energy balances, design equations, and logical constraints. MILP formulations as in (1) have been applied, for example, to the synthesis of heat-exchanger networks (Papoulias and Grossmann, 1983b), synthesis of distillation sequences (Andrecovich and Westerberg, 1985; Raman and Grossmann, 1993), synthesis of utility systems (Papoulias and Grossmann, 1983a), batch plant design (Voudouris and Grossmann, 1992), scheduling (Kondili et al., 1993), and planning problems (Shah and Pantelides, 1991). In practice, process systems typically involve significant uncertainty, arising internally from the process (e.g., flowrate and temperature variations) or from the model (e.g., physical properties, transfer coefficients) or externally to the process (e.g., product demands, envi* Author to whom correspondence should be addressed. Telephone: (44)171-594 6620. Fax: (44)171-594 6606. Email: [email protected]. † Present address: Department of Chemical Engineering, I.T.E.S.M., Monterrey, Mexico. Telephone: (52)(8)3582000 × 5436. E-mail: [email protected]. S0888-5885(96)00451-4 CCC: $14.00

ronmental conditions). Accounting for these types of uncertainty transforms the deterministic problem in (1) into a parametric MILP of the following form:

min c(θ1)x + d(θ2)y y,x

s.t.

(2)

A(θ3)x + E(θ4)y e b + Fθ5 x g 0;

y ∈ {0, 1} θ∈Ξ

where θ ) {θ1, θ2, θ3, θ4, θ5} represents the vector of variations of the parameters of model (1) in a given space Ξ ) {θ|Gθ e g}. The solution of problem (2) has been approached mainly through “multiperiod” or scenario-based formulations (Sahinidis et al., 1989; Shah and Pantelides, 1992), where the uncertain parameters are discretized into a number of deterministic realizations, and stochastic formulations (Reinhart and Rippin, 1986; Ierapetritou and Pistikopoulos, 1994; Liu and Sahinidis, 1996), where the uncertain parameters are described by probabilistic distribution functions. Flexibility considerations have also been included in these models for the design and retrofit of continuous (Pistikopoulos and Grossmann, 1988) and batch processes (Straub and Grossmann, 1990). A different way to incorporate uncertainty into mathematical programs is by using sensitivity analysis and parametric programming techniques. The theory of parametric programming provides a systematic way to analyze the effect of parameter changes on the optimal solution of a mathematical programming model, to quantify the robustness of this solution, and to compare it to other optimal or near-optimal solutions that may arise as the parameters of the model move away from a “nominal” point. For process engineering problems, sensitivity analysis has been mainly used for the analysis of linear and nonlinear continuous models (Takamatsu et al., 1970; Ganesh and Biegler, 1987; Seferlis and Hrymak, 1996). The application of these techniques to mixed-integer problems, however, has been rather limited, mainly due to the computational complexity of the problem and the theoretical limitations of the algorithms presented (e.g., monotonicity, single-parameter changes, etc.). © 1997 American Chemical Society

718 Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997

For mixed-integer linear programs, mainly righthand-side (RHS) parameterizations have been studied. Approaches that have been proposed in the open literature for their solution range from generalizations of the branch and bound algorithm (Roodman, 1974; Marsten and Morin, 1977; Ohtake and Nishida, 1985), use of cutting plane methods (Holm and Klein, 1984; Jenkins and Peters, 1987), and bounding procedures based on the convexity/concavity properties of the parametric LP solution (Geoffrion and Nauss, 1977; Jenkins, 1982); a thorough review of the theory and applications of pMILP can be found in Jenkins (1990). More recently, Pertsinidis (1992) presented two algorithms for the solution of parametric MILP problems and their application to process synthesis problems under uncertainty. All these algorithms, however, share a common limitation: they can only be applied to problems with a single uncertain parameter or several uncertain parameters varying in a single direction at the same time (scalar parametrization). This paper presents a novel algorithm for the multiparametric analysis of process engineering problems under uncertainty, where n parameters in the mixedinteger linear model can vary independently. The algorithm is based on a branch and bound procedure where multiparametric linear programs are solved at each node of the enumeration tree and rigorous bounding procedures are established as a function of the uncertain parameters. The final solution is given as a map of the optimal integer configurations (e.g., optimal structures) and their corresponding optimal value functions, projected onto the uncertain parameter space. The paper is organized as follows. First, an outline of the proposed procedure is presented; an algorithm for the solution of multiparametric linear programs as described by Gal (1995) is briefly explained, and a special procedure for comparing different optimal value functions is proposed. Next, we describe in detail the proposed algorithm for the solution of multiparametric mixed-integer linear problems; theoretical and implementation issues are also discussed. Finally, three examples are presented to illustrate the analytical steps of the proposed procedure as well as its application to process engineering problems. Outline of the Proposed ProceduresTheoretical Developments We consider a special case of the multiparametric mixed-integer linear model for process systems in (2), where only the right-hand-side parameters are allowed to vary, leading to the following problem:

z(θ) ) min cx + dy y,x

s.t.

(3)

Ax + Ey e b + Fθ x g 0;

y ∈ {0, 1} θ∈Ξ

The most common solution technique for deterministic MILPs is the branch and bound algorithm based on linear programming relaxations (see, e.g., Nemhauser and Wolsey, 1988). Branch and bound is a tree enumeration method which relies on the solution of linear

Figure 1. Definition of critical regions.

problems, obtained by relaxing the integrality condition of y at each node of the tree. At each level of the tree, two new nodes are generated by fixing the value of one of the (fractional) integer variables to the two possible solution values (0, 1). To avoid a complete enumeration of the tree, the algorithm alternates between two key steps: (1) solving relaxed problems at each node of the tree to evaluate upper and lower bounds on the final solution and (2) comparing these bounds in order to fathom nodes from which it is proved that the optimal solution will not be found. Following this general strategy, a multiparametric branch and bound (mBB) procedure can be devised as follows. A relaxation of problem (3) can be obtained by considering y as a continuous variable in [0, 1], giving rise to the following multiparametric linear programming formulation:

zj(θ) ) min cx + dyj y,x

s.t.

(4)

Ax + Eyj e b + Fθ x g 0;

0 e yj e 1 θ∈Ξ

At the root node, the solution of problem (4), zj(θ), provides a parametric lower bound on the solution of the multiparametric MILP problem in (3), z(θ). At successive nodes, a subset of 0-1 variables is fixed at the discrete values, providing parametric lower bounds on any integer solution that could be generated from these nodes. On the other hand, the optimal solution of a terminal node in the tree, where all integer variables have been fixed, represents a parametric upper bound on the final solution, denoted by zˆ (θ). These bounds can then be used to eliminate parts of the uncertain space defined for the remaining nodes or to completely fathom a given node. In mBB, the comparison of these bounds must be made considering the different optimal value functions that can be obtained from the solution of (4) as the uncertain parameters θ vary in Ξ. Thus, associated to a node i in the tree, a set of k subspaces or critical regions, CRi,k ∈ Ξi, can be defined such that, in each critical region, a single optimal solution of the linear model exists. This concept is illustrated in Figure 1, where (z* A, x* A) and (z* B, x* B) represent two different sets of optimal value and solution functions. CRA and CRB represent the space in θ where each solution is optimal,

Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997 719

i.e., the critical regions, such that

∀θ ∈ CRA

z* A e z* B CRA )

one can pass from B1 to B2 in one dual (pivot) step. By extension, their corresponding critical regions CR1 and CR2 are said to be neighbors, if and only if B1 and B2 are neighbors, from which it follows that CR1 ∩ CR2 ) 0e . In this way, CRA and CRB in Figure 1 represent two neighboring regions, such that x* A and x* B are given by neighboring optimal bases. The objective is then to find all possible optimal solutions of (4), their optimal value functions zjk(θ), and their corresponding critical regions, CRk, such that

{

C1 g 0 θL e θ e θU

and

∀θ ∈ CRB

z*A g z* B CRB ) θL

{

∪CRk ) K;

C1 e 0 L

θ eθe

θU

θU

where and are the vectors of lower and upper bounds of the uncertain parameters, respectively. A given space k in any node i can then be discarded if one of the following (fathoming) criteria holds: (1) Problem (4) is infeasible for all values of θ ∈ CRi,k (infeasibility criterion). (2) An integer solution is found for all θ ∈ CRi,k (integrality criterion). (3) The solution of the node is greater than the current best upper bound defined in the same space (dominance criterion). Subsequent nodes derived from node i will inherit only those regions which cannot be discarded. Node i will be fathomed only if all its k critical regions can be discarded. Solution of the Multiparametric Linear Problems. An algorithm for the solution of RHS multiparametric linear problems (mpLP) proposed by Gal and Nedoma (1972) is used here for the solution of the relaxed problems in (4) at each node of the mBB tree. The procedure, based on the Simplex algorithm for deterministic LPs, incorporates the matrix of coefficients (F) of the uncertain parameters into the Simplex tableau, extending the optimality conditions to become a function of vector θ. The three major steps in the algorithm are (i) find an initial vector of parameter values, θ0, for which a feasible solution exists, (ii) find an initial optimal solution of the deterministic LP obtained by fixing the value of θ ) θ0, and (iii) by pivoting over the initial optimal tableau, find the set of optimal solutions that define zj(θ), ∀θ ∈ Ξ. These steps are based on the definition of the optimality conditions of the parametric problem as follows. Let B be an optimal basis and F the related index set of basic variables. The primal feasibility condition in terms of θ is then given by: -1 B-1 F b(θ) ) BF (b + Fθ) g 0

or equivalently

-FFθ e bF FF ) B-1 F F;

(5)

bF ) B-1 F b

Since the dual solution does not depend on θ for RHS parametrizations, a critical region, i.e., the space in Ξ where B remains optimal, can be uniquely defined by the conditions in (5) and the set of constraints Gθ e g. Two such optimal bases, B1 and B2, are said to be neighboring optimal bases or neighbors if and only if (i) there exists θ* ∈ Ξ such that B1 and B2 are both optimal bases to the multiparametric LP for θ* and (ii)

K∈Ξ

where k is the index of critical regions and K is the set of all feasible vector parameters, that is, all possible values of θ for which problem (4) has a finite optimal solution, provided that the various CRk do not overlap. A detailed algorithm is presented in Appendix I; more details on the mathematical background and the algorithm itself are given in the book by Gal (1995). Comparison Procedure. In order to apply the dominance criterion, the solution of problem (4) at any given node i in the tree, zji(θ), must be compared to the current upper bound, zˆ (θ), so as to find the space of θ values where zji(θ) e zˆ (θ). However, these two objective value functions can only be compared in the space of θ where both functions are defined, that is, the intersection of their corresponding critical regions. The first step in the procedure is then to find CRint ) CRi ∩ CRUB. Following the algorithm of Gal and Nedoma (1972) for multiparametric LPs, the solution zji(θ) will be defined as a piecewise linear function, formed by k optimal value linear functions, zji,k(θ), and their corresponding critical regions, CRi,k, describing the optimality space of each function. Equivalently, the upper bound will also be defined by a set of linear functions, zˆ k′(θ), and their corresponding critical regions of optimality, UB . CRk′ By definition (see Appendix I), each of these critical spaces is given by closed nonempty polyhedra, that is, a set of linear inequalities in θ; it then follows that CRint can be defined by the set of constraints from both CRi,k UB and CRk′ . The definition of this new region CRint is graphically depicted in Figure 2, where the constraints defining CRi,k are shown in solid lines and the conUB are given by dashed lines. straints defining CRk′ Mathematically, this set of constraints can be expressed as the equivalent set of nonredundant constraints, which, in turn, can be identified through a redundancy test for linear inequalities (see Appendix II). In general, four cases can result from such a redundancy test: Case 1 (Figure 2a): All the constraints from CRi,k UB are redundant. This implies that CRi,k ⊃ CRk′ , and UB int therefore CR ) CRk′ . UB Case 2 (Figure 2b): All the constraints from CRk′ UB are redundant. This implies that CRi,k ⊂ CRk′ , and therefore CRint ) CRi,k. Case 3 (Figure 2c): Constraints from both regions are nonredundant. The two spaces intersect with each other, and CRint is given by the space delimited by the nonredundant constraints. Case 4 (Figure 2d): The problem is infeasible. The two spaces are apart from each other; CRint) {0e }. If there exists an intersection space between the two regions, CRint * {0e }, the second step of the procedure

720 Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997

Figure 2. Comparing two critical regions.

is to check whether zji,k(θ) e zˆ k′(θ) in CRint. This can be achieved by defining a new constraint given by

zdiff(θ) ) zji,k(θ) - zˆ k′(θ) e 0

(6)

and solving a new redundancy test for this constraint in CRint. The solution of this problem will now belong to one of three cases as depicted in Figure 3: Case 1 (Figure 3a): The new constraint is redundant; this would mean that the solution of the node in CRint is less than the upper bound, and therefore the space must be kept for further analysis. Case 2 (Figure 3b): The new constraint generates an empty region; this would mean that the solution of the node in CRint is greater than the upper bound, and therefore the space can be discarded from the analysis. Case 3 (Figure 3c): The new constraint is nonredundant; this would mean that the solution of the node is less than the upper bound in only part of CRint, and therefore the rest of the space can be discarded from the analysis. Based on the previous considerations, a formal procedure to compare the solution of node i to a current upper bound would be as follows: 1. Select a critical region k from the solution of the node and a space k′ from the current upper bound. UB 2. Set CRint ) CRi,k ∩ CRk′ . If the intersection space of the two regions is empty, go back to step 1. 3. Delete the redundant constraints from CRint. 4. Set zdiff(θ) ) zji,k(θ) - zˆ k′(θ). Add the constraint in (6) to the set of constraints defining CRint and check if the new constraint is redundant. One of the following cases must occur: (a) zji,k(θ) e zˆ k′(θ), ∀θ ∈ CRint: keep

CRint in the solution of the node. (b) zji,k(θ) g zˆ k′(θ), ∀θ ∈ CRint: CRint can be removed from further analysis if the rest of CRi,k can be defined, i.e., continue for CRi,k UB - (CRi,k ∩ CRk′ ). (c) zji,k(θ) e zˆ k′(θ) for only some θ ∈ CRint: partition CRi,k in two (smaller) spaces: (i) the space for which zji,k(θ) e zˆ k′(θ), i.e., CRint with the additional constraint in (6), and (ii) the part of CRi,k UB that does not intersect with CRk′ , i.e., CRi,k - (CRi,k ∩ UB CRk′ ). Update the critical regions according to the result. 5. Repeat the procedure for all necessary combinations of k and k′. It should be noted that a region CRint can only be deleted from further analysis if the rest of the original space CRi,k, i.e., (CRi,k - CRint), can be sufficiently defined as a new polyhedron or a set of polyhedra. In general, it can be proved that this can be guaranteed if UB CRi,k ⊂ CRk′ or if the common predecessor space (a common predecessor space corresponds to a critical region CRp, defined in a node at a higher level of the UB tree, such that CRi,k ⊂ CRp and CRk′ ⊂ CRp) of CRi,k UB and CRk′ is completely covered by neighboring regions in CRUB. This condition ensures that the different spaces defined by the intersection of CRi,k and the UB various CRk′ can be exactly defined for each k′ as new polyhedra and therefore (CRi,k - CRint) can also be defined. A Branch and Bound Algorithm for Multiparametric MILPs Based on the theoretical developments presented in the previous section, an mBB algorithm for the solution

Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997 721

Figure 3. Possible outcomes from the bounding procedure.

Figure 4. Algorithm for multiparametric MILPs.

of RHS parametric MILP problems as in (3) is formally stated in this section. As shown in Figure 4, it involves the following basic steps: 1. Initialize the list of nodes to solve, N, as an empty set and its counter n ) 0. Set the upper bound zˆ (θ) ) ∞.

2. Solve the fully relaxed problem in (4). If no solution is found, stop the procedure; the problem has no integer solution. If any critical region presents an integer solution, update the upper bound and discard the region from further analysis. 3. Select a 0-1 variable from yj and generate two new nodes fixing the selected variable to the discrete values 0 and 1. Set the uncertain parameter space of these nodes, Ξn+1 and Ξn+2, equal to the remaining critical regions. Add the new nodes to N, and update n ) n + 2. 4. If N is empty, stop the procedure with the solution as given by zˆ (θ). Otherwise, take node i from N and remove it from the list; n ) n - 1. 5. Solve problem (4) with the subset of fixed integer variables and uncertain parameter space Ξi as defined for node i. If the problem is infeasible for all θ ∈ Ξi, go back to step 4. 6. With the procedure proposed in the previous section, compare the solution of node i to the current upper bound. 7. If the node is terminal (all values of y are fixed), update the upper bound function according to the solution of the comparison procedure. 8. If every critical region from the solution of node i is removed or the node represents an integer solution, return to step 4; otherwise, return to step 3. Note that, since the definition of critical regions is inherited from node to node, some of these regions, created by the relaxation of the integer variables in the upper nodes of the tree, may result in the same optimal solution when these variables are fixed; if all these

722 Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997

regions are neighbors covering completely a single critical region in a predecessor node, they can be aggregated in a single larger space by removing the constraints defining the hyperplanes separating them. The resulting space will be equal to the original critical region. For the case when the same optimal solution is found in two or more critical regions, but the regions cannot be aggregated, it may be advantageous when defining a new upper bound to re-solve the relaxation problem at this node considering the entire initial uncertainty parameter space. This will yield the minimum number of critical regions covering the widest possible uncertain space. Furthermore, the initial optimal solution (tableau) obtained at this node would still be valid, and therefore only the definition of critical spaces and neighbor regions need be reevaluated. If a new integer solution is found but neither CRi,k UB nor CRk′ can be discarded or redefined, then the new upper bound will consist of three distinct areas which must be evaluated in the following order: (1) CRint including the constraint in (6) if nonredundant, (2) UB CRk′ , and (3) CRi,k. This, however, does not restrict the order used in the comparison procedure, although, from a computational viewpoint, starting from the largest space is expected to be more efficient. Further refinements can also be implemented in order to improve the performance of the algorithm. In particular, several heuristic rules already available for the selection of the next candidate subproblem (node selection) and the new variable upon which we may branch at a given level of the enumeration tree (branching variable selection) could, in principle, be applied. The performance of such rules on parametric problems, however, should be further studied. Also note that, at the limiting case when the number of uncertain parameters is exactly 1, this algorithm would reduce to that of Ohtake and Nishida (1985). In the case of a single parameter, however, there are two major simplifications as compared to the multiparametric case: first, the critical regions are given by “intervals” in θ which can be evaluated directly from the dual solution of the equivalent deterministic problem (see e.g. Gal, 1995), and, second, since the optimal solution at each node is given as a function of a single parameter, the comparison procedure is reduced to finding the intersection between two linear functions in a onedimensional space. As the number of uncertain parameters increases to at least two uncertain parameters, these characteristics disappear and thus the identification of the critical regions and the comparison of two different optimal solutions become the core of the new algorithm. Finally, note that, as in any standard branch and bound algorithm, the convergence of the proposed procedure in a finite time is ensured by exhaustiveness, since the solution of each node will yield a finite number of critical regions and the maximum number of nodes is bounded by the finite number of integer variables. Numerical Examples Three examples are presented in this section to illustrate the basic steps and the potential of the proposed algorithm. The first two are modifications of the numerical examples presented in Gal (1995), to which integer variables have been included, together with some logical constraints added to connect the

integer variables to continuous decisions. The third example corresponds to a process synthesis under uncertainty problem for the selection of a utility plant at minimum cost. In the case of the first two examples, the branching variables were selected following their order in the model in a first-in-first-out (FIFO) searching strategy, whereas in the last example, a simple last-in-first-out (LIFO) strategy was adopted for the tree search, while for the branching variable, the integer variable with a relaxed value closer to 0.5 was selected. Example 1. Consider the following numerical example:

min - 3x1 - 8x2 + 4y1 + 2y2 x

s.t.

x1 + x2 e 13 + θ1 5x1 - 4x2 e 20 -8x1 + 22x2 e 121 + θ2 4x1 + x2 g 8 x1 - 10y1 e 0 x2 - 15y2 e 0 x g 0;

y ∈ {0, 1}

0 e (θ1, θ2) e 10 Following the proposed mBB algorithm, the solution of this problem is given by the following: Node 0 (Fully Relaxed Problem). 1. Set Ξ ) {θ|0 e θ e 10}, zˆ ) ∞. 2. A feasible solution is found for θ0 ) (0, 0). 3. Solving the problem at θ0, we obtain one optimal function given by

zj0,1 ) -73.3 - 4.00444θ1 - 0.17556θ2 which is valid in the region CR0,1 given by the following constraints:

0.07333θ1 - 0.00333θ2 e 0.45 θ2 e 10 The optimal solution related to this optimal function is given by

x1(θ) ) 5.5 + 0.73θ1 - 0.03θ2 x2(θ) ) 7.5 + 0.27θ1 + 0.03θ2 y1(θ) ) 0.55 + 0.07θ1 - 0.004θ2 y2(θ) ) 0.5 + 0.02θ1 + 0.002θ2 Note that, since CR0,1 corresponds only to a part of the original uncertain space and the solution has no neighboring bases, the space left out of the analysis represents an infeasible space.

Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997 723

Node 1 (Fix y1 ) 0). 1. Set Ξ1 ) CR0,1. 2. The problem is found infeasible in Ξ1; therefore, this node is fathomed. Node 2 (Fix y1 ) 1). 1. Set Ξ2 ) CR0,1. 2. A feasible solution is found for θ0 ) (0, 0). 3. Solving the problem at θ0, we obtain one optimal function given by

variables, two continuous variables, and three uncertain parameters varying in the range [0, 5].

zj2,1 ) -71.5 - 4.29778θ1 - 0.16222θ2

x1 + 2x2 e 12 + θ1 - θ3

which is valid in the region CR2,1, which coincides with the original space CR0,1. Node 3 (Fix y1 ) 1; y2 ) 0). 1. Set Ξ3 ) CR2,1. 2. Solving the problem at θ0 ) (0, 0), we obtain an optimal function given by z3,1) -8 for all θ in Ξ3. Its critical region CR3,1 is given by the same constraints. 3. Since this node represents an integer solution, the upper bound is updated; i.e., y*) (1, 0), zˆ ) zj3,1, and CRUB ) CR3,1. Node 4 (Fix y1 ) 1; y2 ) 1). 1. Set Ξ4 ) CR2,1. Solving the problem at θ0 ) (0, 0), we obtain one optimal function given by

zj4,1 ) -70.5 - 4.33333θ1 - 0.16667θ2 which is valid in the region CR4,1 given by the same constraints covering the entire space Ξ4. No neighboring optimal bases are found. 2. Comparing this solution with the upper bound, CRint ) CR4,1 ) CRUB, and the new constraint is given by

min - 3x1 - 2x2 + 10y1 + 5y2 x

s.t.

x1 e 10 + θ1 + 2θ2 x2 e 10 - θ1 + θ2 x1 + x2 e 20 - θ2

x1 - 20y1 e 0 x2 - 20y2 e 0 -x1 + x2 g 4 - θ3 y1 + y2 g 1 x g 0;

y ∈ {0, 1}

0eθe5 1. The fully relaxed problem with Ξ ) {θ|0 e θ e 5} is first considered. A feasible solution is found for θ ) (0, 0, 0), where the optimal solution of the problem is given by

zj0,1 ) -9.33333 - 1.58333θ1 + 0.41667θ3

{

4θ1 - 3θ2 - 2θ3 e 14 θe5 ∀θ

CR0,1 )

By applying one dual step, one neighbor basis is identified with an optimal value function and corresponding space of optimality given by

zj0,2 ) -31.5 + 4.75θ1 - 4.75θ2 - 2.75θ3 zdiff ) zj4,1 - zˆ ) -62.5 - 4.33333θ1 - 0.16667θ2 e 0 This constraint is obviously redundant for any positive value of θ, and, therefore, z4,1 represents the new upper bound, with CR4,1 defining its optimality space. Since no more nodes can be created, this upper bound represents the final optimal solution of the problem given by

{

4θ1 - 3θ2 - 2θ3 g 14 θ1 e 5

CR0,2 )

Note that this solution covers the entire initial uncertainty space, with the two regions separated by one plane. 2. A new node is created by fixing y1 ) 0, Ξ1,1 ) CR0,1, and Ξ1,2 ) CR0,2. The solution of this node is given by

zj1,1 ) -7 - θ1 + θ3 y* ) (1, 1) CR1,1 ) z*(θ) ) -70.5 - 4.33333θ1 - 0.16667θ2

{

3θ1 - 2θ2 - θ3 e 8 θe5 ∀θ

zj1,2 ) -15 + 2θ1 - 2θ2 CR )

{

0.07333θ1 - 0.00333θ2 e 0.45 θ2 e 10

Note that, in this numerical example, only one optimal solution is found at each node and the final solution also corresponds to a single 0-1 integer vector combination. Example 2. In the following example, the comparison procedure is illustrated in detail, and its several cases are discussed. The problem consists of two integer

{

-3θ1 + 2θ2 + θ3 e -8 CR1,2 ) 4θ1 - 3θ2 - 2θ3 e 14 θe5 θ ∈ {θ1, θ3} zj1,3 ) -15 + 2θ1 - 2θ2 CR1,3 )

{

-4θ1 + 3θ2 + 2θ3 e -14 θ1 e 5

Since the node presents an integer solution, i.e., y2 ) 1 for all values of θ in CR1,1, CR1,2, and CR1,3, this solution

724 Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997

represents an upper bound of the final optimal solution. Moreover, since there are two regions (CR1,2 and CR1,3) with the same optimal value function (zj1,2 ) zj1,3), a more compact representation of the regions of optimality to define the new upper bound can be obtained by resolving the node for the entire initial uncertain space. This then yields the following upper bound:

zˆ 1 ) -7 - θ1 + θ3 CRUB 1

{

1.5θ1 - θ2 - 0.5θ3 e 4 ) θe5 ∀θ

zˆ 2 ) -15 + 2θ1 - 2θ2 CRUB 2 )

{

-1.5θ1 + θ2 + 0.5θ3 e -4 θe5 θ ∈ {θ1, θ3}

3. A new node is created by fixing y1 ) 1, Ξ2,1 ) CR0,1, and Ξ2,2 ) CR0,2. The solution of this node yields:

Adding the new constraint to the intersection set results again in a nonredundant constraint, yielding a solution given by

zj2,1b ) -3.33333 - 1.58333θ1 + 0.16667θ3

{

CR2,1b )

1.33333θ1 - θ2 - 0.66667θ3 e 4.66667 -3θ1 + 2θ2 + θ3 e -8 -3.58333θ1 + 2θ2 + 0.166667θ3 e 11.66667 θe5 ∀θ

5. Taking CR2,2, it is obvious from the previous analysis that it does not intersect with CRUB 1 , so this option is not considered. Comparing it then to CRUB 2 , we find that CR2,2 ⊂ CR2UB and therefore CRint ) CR2,2. The difference between the two optimal objective values results in a new constraint:

zdiff ) -10.5 + 2.75θ1 - 2.75θ2 - 3θ3 e 0 which is found nonredundant, yielding a new, smaller, space given by

zj2,1 ) -3.33333 - 1.58333θ1 + 0.16667θ3

CR2,1 )

{

1.33333θ1 - 1θ2 - 0.66667θ3 e 4.66667 θe5 ∀θ

zj2,2 ) -25.5 + 4.75θ1 - 4.75θ2 - 3θ3

CR2,2 )

{

-4θ1 + 3θ2 + 2θ3 e -14 θ1 e 5

The lower bound given by this solution is now compared to the current upper bound to define the space in Ξ that can be removed by applying the dominance criterion. 4. Comparing CR2,1 and CRUB 1 , it is found that int ) CRUB. Comparing CRUB ⊂ CR and therefore CR 2,1 1 1 the two optimal objective functions:

zdiff ) zj2,1 - zˆ 1 ) 3.66667 - 0.58333θ1 - 0.83333θ3 Adding the new constraint {zdiff e 0}, it is found nonredundant. CR2,1 is then divided into two spaces. The first one is given by CRUB 1 , with the new constraint included as

zj2,2a ) -25.5 + 4.75θ1 - 4.75θ2 - 3θ3

{

-4θ1 + 3θ2 + 2θ3 e -14 CR2,2a ) 2.75θ1 - 2.75θ2 - 3θ3 e 10.5 θ1 e 5 The mBB procedure continues with CR2,1a, CR2,1a, and CR2,2a as parameter spaces of further interest in the analysis. 6. A new node is created by fixing y1 ) 1, y2 ) 0, and Ξ3 ) CR2. Its solution results in one infeasible space, CR2,2a, and a new solution for the two regions, CR2,1a and CR2,1b. Comparing this solution to the current upper bound, it is found that z3 g zˆ for all possible values of θ; therefore, this node is fathomed. 7. A new node is created by fixing y1 ) 1, y2 ) 1, and Ξ4 ) CR2. Its solution for the three parameter spaces in yields

zj4,1 ) 0.33333 - 1.66667θ1 + 0.33333θ3

{

1.5θ1 - θ2 - 0.5θ3 e 4 CR4,1 ) -0.58333θ1 - 0.83333θ3 e 3.66667 θe5 ∀θ zj4,2 ) 0.33333 - 1.66667θ1 + 0.33333θ3

{

CR4,2 ) zj2,1a ) -3.33333 - 1.58333θ1 + 0.16667θ3

{

1.5θ1 - θ2 - 0.5θ3 e 4 CR2,1a ) -0.58333θ1 - 0.83333θ3 e 3.66667 θe5 ∀θ The second space is given by (CR2,1 ∩ CRUB 2 ). In this region, however, z2,1 must be lower than zˆ 2, which generates a new constraint: diff

z

) 11.66667 - 3.58333θ1 + 2θ2 + 0.166667θ3 e 0

1.33333θ1 - θ2 - 0.66667θ3 e 4.66667 -3θ1 + 2θ2 + θ3 e -8 -3.58333θ1 + 2θ2 + 0.166667θ3 e 11.66667 θe5 θ ∈ {θ1, θ3} zj4,3 ) -23 + 5θ1 - 5θ2 - 3θ3

{

-4θ1 + 3θ2 + 2θ3 e -14 CR4,3 ) 2.75θ1 - 2.75θ2 - 3θ3 e 10.5 θ1 e 5 8. Comparing this solution to the current upper bound CRUB, we identify two regions in CR2UB where

Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997 725

Figure 6. Enumeration tree for example 3. Table 1. Mathematical Model for Example 3 cost functions

Figure 5. Superstructure for the utility plant synthesis.

the new integer solution has a better optimal value than the previous one. The upper bound is then updated, resulting in the following five independent regions:

mass balances

y* ) (0, 1) zˆ 1 ) -7 - θ1 + θ3 CRUB 1 )

{

{

energy balances

1.5θ1 - θ2 - 0.5θ3 e 4 θe5 ∀θ

zˆ 2 ) -15 + 2θ1 - 2θ2

-3θ1 + 2θ2 + θ3 e -8 UB CR2a ) 3.66667θ1 - 2θ2 - 0.33333θ3 e 15.33333 θe5 θ ∈ {θ1, θ3}

{

demands

logical constraints

-3θ1 + 2θ2 + θ3 e -8 UB CR2b ) -3θ1 + 3θ2 + 3θ3 e -8 θe5 θ ∈ {θ1, θ3} y* ) (1, 1) zˆ 3 ) 0.33333 - 1.66667θ1 + 0.33333θ3 CRUB 3 )

{

1.33333θ1 - θ2 - 0.66667θ3 e 4.66667 -3.66667θ1 + 2θ2 + 0.33333θ3 e -15.33333 θe5 θ ∈ {θ1, θ3} zˆ 4 ) -23 + 5θ1 - 5θ2 - 3θ3

{

-4θ1 + 3θ2 + 2θ3 e -14 CRUB 4 ) 3θ1 - 3θ2 - 3θ3 e 8 θ1 e 5 which represents the final solution to the parametric problem. Example 3. Synthesis of a Utility Plant. The last example (taken from Biegler et al., 1996) corresponds to the problem of synthesizing a utility plant for uncertain demands of electricity (E) and steam at medium and low pressures (MP, LP). The superstructure of the problem, shown in Figure 5, contains the following alternatives. Steam at different pressures (HP ) 4.83 MPa, MP ) 2.07 MPa, LP ) 0.34 MPa) can be raised with high- and medium-pressure boilers; let-down valves can be used. For the turbines, it is assumed that high-pressure turbines can be expanded down to low pressure and can have extractions to medium pressure. Medium to low turbines are backpressure turbines.

Cost ) CHPB + CMPB + CHPT + CMPT + CExt CHPB ) 90000YHPB + 9600HP CMPB ) 40000YMPB + 8500MP CHPT ) 45000(YT1 + YT2) + 25EHPT CHPT ) 25000(YT3 + YT4) + 14.5EMPT CExt ) 20000(YT1′ + YT2′) HP ) HPMP1 + HPLP1 + HPMP2 + HPLP2 + HPLD MPProd ) MP + HPMP1 + HPMP2 + HPLD - MPLP1 - MPLP2 - MPLD LPProd ) HPLP1 + HPLP2 + MPLP1 + MPLP2 + MPLD EHPT ) EHP1 + EHP2 EMPT ) EMP1 + EMP2 EHP1 ) 42.6HPMP1 + 121HPLP1 EHP2 ) 42.6HPMP2 + 121HPLP2 EMP1 ) 78.4MPLP1 EMP2 ) 78.4MPLP2 MPProd g 20 + 10θ1 LPProd g 80 + 10θ4 EHP1 + EMP1 g 7000 + 1000θ2 EHP2 + EMP2 g 4000 + 1000θ3 HPLP1 - MaxFlowYT1′ e 0 HPLP2 - MaxFlowYT2′ e 0 HPMP1 + HPLP1 - MaxFlowYT1 e 0 HPMP2 + HPLP2 - MaxFlowYT2 e 0 MPLP1 - MaxFlowYT3 e 0 MPLP2 - MaxFlowYT4 e 0 YT1 - YHPB e 0 YT2 - YHPB e 0 YT1′ - YT1 e 0 YT2′ - YT2 e 0 YT1 + YT3 ) 1 YT2 + YT4 ) 1

Table 2. Uncertain Demands for the Utility Plant demand

nominal value

variation

MP steam (θ1) power 1 (θ2) power 2 (θ3) LP steam (θ4)

25 ton/h 7500 kW 4500 kW 85 ton/h

(5 (500 (500 (5

Power demands can be satisfied with any of these turbines, but only one turbine can be assigned to each demand. The model, based on the MILP formulation of Papoulias and Grossmann (1983a), is given in Table 1. The multiparametric MILP consists of 30 constraints, 8 integer variables, 24 continuous variables, and 4 uncertain parameters, defining the range of values that the electricity and steam demands may take. These uncertain demands and their variations are given in Table 2; the nominal values of these parameters are taken from the original example. The enumeration tree for the solution of this problem following the proposed mBB algorithm is shown in Figure 6. Each node has been enumerated according to the order in the solution procedure; the subscript indicates the number of critical regions defined in the

726 Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997 Table 3. Optimal Structures and Their Critical Regionsa (1, 0, 0, 1, 1, 0, 0, 0)

Cost1(θ) ) 130.06 + 9.60θ1 + 13.69θ2 + 2.50θ3 CR1 ) -12.75θ2 + 10θ4 e 9.28 Cost2(θ) ) 121.15 + 9.60θ1 + 1.45θ2 + 2.50θ3 + 9.60θ4

{

-10θ1 + 23.47θ3 - 10θ4 e 6.10 12.75θ2 - 10θ4 e -9.28 Cost3(θ) ) 115.29 + 1.45θ2 + 25.03θ3

CR2 )

(1, 0, 1, 0, 0, 1, 1, 0)

CR3 )

{

CR4 )

{

10θ1 - 23.47θ3 + 10θ4 e -6.10 12.75θ2 - 10θ4 e -9.28 Cost4(θ) ) 141.01 + 62.20θ1 + 10.43θ2 + 9.38θ3 -3.38θ1 - 3.26θ2 + 6.88θ3 e 0.05 3.52θ1 - 8.26θ2 - 8.26θ3 + 10θ4 e 3.86

Cost5(θ) ) 137.30 + 9.60θ1 + 2.50θ2 + 1.45θ3 + 9.60θ4 CR5 ) (1, 0, 0, 1, 1, 0, 0, 1)

-11.19θ2 - 1.05θ3 + 9.60θ4 e 3.76 -5.43θ1 + 12.75θ2 + 12.75θ3 - 15.43θ4 e -5.97

Cost6(θ) ) 119.86 + 6.22θ1 + 9.38θ2 + 10.43θ3 CR6 )

a

{ {

6.22θ1 + 7.93θ2 - 14.60θ3 e 11.57 12.75θ2 - 10θ4 e -9.28

Cost ) $/yr 104.

node, while the superscript indicates the number of these regions that yield an integer solution. Nodes where the upper bound was updated are indicated by shadowed circles; crossed circles denote infeasible nodes for any value of θ in the corresponding parameter space. The order of the vector of integer variables y starts with the boilers (HPB, MPB), then the high-pressure turbines (T1, T2), the low-pressure turbines (T3, T4), and the extensions (T1′, T2′). The final solution is given in Table 3, where the three optimal integer configurations are tabulated with their optimal value functions and corresponding critical regions (initial bounds on θ have been omitted). For the first structure, y*1 ) (1, 0, 0, 1, 1, 0, 0, 0), the critical regions presented cover the entire uncertain space, while the other two structures are given in terms of regions embedded in the previous solution; both regions of optimality for y* 2 ) (1, 0, 1, 0, 0, 1, 1, 0) are in the first optimal region of y*1, while the region of optimality of y* 3 ) (1, 0, 0, 1, 1, 0, 0, 1) is embedded in the third critical region of y* 1. The information provided by the parametric solution as given in Table 3 can then be used in several ways: (i) a complete analysis of the sensitivity of the “optimal” structures (integer solutions) obtained can be easily performed with respect to each uncertain parameter, which can then be used to make a decision on the final implementation/design (for example, select the most robust structure/design); (ii) “important” uncertain parameters can be determined, thus providing guidelines for effective control so as to minimize, if possible, variations in the plant; (iii) Table 3 essentially represents a map of all scenarios of the optimal operation of the plant for any possible realization of values of the uncertain parameterssin this respect, it could be used to define (on-line) optimal operating conditions (for example, scheduling) given any value of the uncertainty; (iv) it can be used as a vehicle for the computation of an expected objective value in the context of stochastic programming (see Acevedo and Pistikopoulos, 1996b). For comparison, a deterministic MILP problem was obtained by fixing the value of the uncertain parameters at θ ) 0. The solution of this problem following the

same branching sequence requires the evaluation of 19 nodes; the extra nodes solved to obtain the multiparametric solution were nodes 6 and 7 (both infeasible) and nodes 16 and 17. In general, it can be expected that the number of nodes needed to obtain a multiparametric solution will increase with respect to that of the equivalent deterministic problem. Furthermore, the parametric solution of each node, in addition, requires the evaluation of the critical regions and the redundancy test (employed for comparing two different solutions), thus increasing the solution times. The size of these latter problems, however, will be much smaller in most cases, involving as variables only the uncertain parameters and, in the case of the comparison procedure, typically very few constraints. Concluding Remarks A novel algorithm has been presented for the parametric solution of mixed-integer linear problems, enabling the analysis of the effect of independent variations of several parameters on the optimal solution. The solution procedure, based on the branch and bound algorithm for deterministic MILPs, involves two main calculation steps: (i) the solution of multiparametric LP problems at each node of the tree and (ii) the evaluation of the uncertain parameters space for which a node must be solved. A procedure for comparing optimal value functions in terms of the vector of uncertain parameters was also proposed to guide the search and eliminate parts of the tree. The approach was illustrated through some example problems. Although computational experience with larger models is still underway, a good performance of the algorithm is expected compared to the alternative of solving a (possibly large) number of deterministic MILPs for different values of the uncertain parameter vector; moreover, the parametric procedure will render complete information regarding the sensitivity of the optimal solution(s) with respect to changes on the uncertain parameters. The extension of this solution approach to mixedinteger nonlinear formulations can also be envisaged,

Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997 727

coupling, for example, this algorithm with a (scalar) parametric outer-approximation procedure as proposed by Acevedo and Pistikopoulos (1996a). Acknowledgment J.A. gratefully acknowledges financial support from the Mexican Government through the Consejo Nacional de Ciencia y Tecnologı´a, CONACyT. Financial support from the Commission of European Communities under Grants ERB CHBI CT93 0484 and JOE3 CT95 0017 is also acknowledged. Nomenclature Subindex i ) index of the nodes in the enumeration tree k ) index of critical regions Sets and Constants A ) matrix of constraint coefficients of the continuous variables B ) optimal basis of an LP E ) matrix of constraint coefficients of the integer variables F ) matrix of constraint coefficients of the uncertain parameters in the RHS G ) matrix of constraint coefficients of the uncertain parameters in the LHS of the purely parametric constraints b ) vector of RHS values of constraints c ) vector of cost coefficients of the continuous variables d ) vector of cost coefficients of the integer variables g ) vector of RHS values of purely parametric constraints Variables x ) continuous variables y ) 0-1 integer variables y* ) optimal value of the 0-1 integer variables yj ) relaxed integer variables z(θ) ) objective function of the mpMILP problem zi(θ) ) objective function of the mpMILP problem at node i zji(θ) ) objective function of the relaxed mpMILP problem at node i zˆ i(θ) ) best upper bound for the mpMILP problem CR ) critical region θ ) vector of uncertain parameter Ξ ) initial uncertain parameter space Ξi ) space in Ξ to consider in node i

Appendix I. Algorithm for mpLP Problems (Gal, 1995) An algorithm for solving parametric LP problems is given by Gal and Nedoma (1972). The basic principles and mathematical foundations are given here without repeating proofs that can be found in their original work. Further details can also be found in the excellent book by Gal (1995). Definition 1. If there exists a vector θ such that the multiparametric LP problem has a finite optimal solution with respect to this θ, then this θ is said to be a feasible vector parameter. Let K denote the set of all feasible vector parameters. Definition 2. The basis BF is said to be an optimal basis if and only if there exists θ ∈ K such that BF is primal and dual feasible with respect to θ. Theorem 1. Given a fixed θ0 ∈ Rs, if there exists a finite optimum to the multiparametric LP problem, then, for every θ ∈ ReS, the problem has a finite optimal solution or has no feasible solution.

Theorem 2. K is a closed convex polyhedral set. Theorem 2 implies that it is possible to express K as a union of a finite number of closed convex polyhedral sets CRk. By the proof of Theorem 1, it can also be stated that the dual solution of the multiparametric linear problem does not depend on θ, and therefore the optimality of BF with respect to θ is given by the primal feasibility conditions, thus uniquely defining CRF. Theorem 3. The value function zj(θ) is linear for all θ ∈ CRF and, thus, continuous over K. Definition 3. Two bases B1 and B2 are said to be neighboring bases (or neighbors) if and only if (1) there exists θ* ∈ K such that both B1 and B2 are optimal bases to (4) for θ* and (2) it is possible to pass from B1 to B2 by one dual step (and vice versa). Definition 4. Consider two optimal bases B1 and B2 and CR1 and CR2, the corresponding feasible regions of θ as defined by (5). The two regions CR1 and CR2 are said to be neighbors if and only if B1 and B2 are neighbors. Based on these definitions and theorems, a procedure for the solution of problem (4) is now stated. The procedure consists of two parts. First, an initial optimal solution is found for fixed (feasible) values of the uncertain parameters. In the second part, all neighbor bases of the current optimal tableau are found in a recursive procedure, determining a set of k nonoverlapping regions CRk such that

∪CRk ) K;

K∈Ξ

with K as given in Definition 1. The steps of the procedure for solving parametric LPs can then be summarized as follows: Phase 1: Finding an Initial Solution. 1. Find a feasible point (x0, θ0) by considering θ as a variable in problem (4). If no feasible point can be found, then stop; problem is infeasible for any θ. 2. With the feasible point, fix θ ) θ0 and solve the resulting LP. Phase 2: Finding All Optimal Solutions. 1. Form two lists of optimal bases: list V keeps the optimal basis for which all possible neighbor bases have been found, and list W keeps the optimal bases which have been defined (as neighbors of those bases in V) but which have not been examined for possible neighbors. V is initially empty, while W is initialized with the optimal bases found in phase 1. 2. Take an element of W (basis i) and find all possible neighbor bases. Insert in W those bases which are not in V or W. 3. By one dual step, move from basis i to each neighbor to determine its optimal solution. 4. Repeat the procedure until W is empty. V then represents the set of optimal bases that defines the solution of problem (4). Although the complete procedure has not been fully automated in this work, Gal (1995) presents a detailed study on the implementation and application of this algorithm for finding the different optimal value functions and the corresponding (critical) regions of optimality. The implementation of these ideas will be essential for the evaluation of the proposed multiparametric MILP algorithm for the solution of larger models.

728 Ind. Eng. Chem. Res., Vol. 36, No. 3, 1997

Appendix II. Redundancy in Linear Programming Consider the system of linear inequalities: N

ai,jxj e bi ∑ j)1

i ) 1, ..., m

A constraint k from the set i N

ak,jxj e bk ∑ j)1 is redundant with respect to the set of feasible solutions if there exists an optimal solution to the following problem:

min k N

s.t.

ai,jxj + i ) bi ∑ j)1

i ) 1, ..., m

such that k is a basic variable and the coefficients of the row where k occurs in the optimal basis are nonpositive (except the coefficient of k). If min k > 0, then the constraint is said to be strongly redundant; if min k ) 0, then it is weakly redundant. Note that, in our case, the rigorous minimization of each k is not necessary, since any feasible solution where k ) 0 is a sufficient condition to discard the constraint as redundant. An efficient algorithm for determining redundant constraints, following the same procedure as that for finding neighboring regions, is given in Chapter 2 of Gal’s book (1995). Literature Cited Acevedo, J.; Pistikopoulos, E. N. A Parametric MINLP Algorithm for Process Synthesis Problems under Uncertainty. Ind. Eng. Chem. Res. 1996a, 35, 147. Acevedo, J.; Pistikopoulos, E. N. A Hybrid Parametric/Stochastic Programming Approach for Mixed-Integer Linear Problems under Uncertainty. Ind. Eng. Chem. Res. 1996b, accepted for publication. Andrecovich, M. J.; Westerberg, A. W. A Simple Synthesis Method Based on Utility Bounding for Heat-Integrated Distillation Sequences. AIChE J. 1985, 31, 363. Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods for Chemical Process Design; book in preparation, 1996. Gal, T. Post-optimal Analyses, Parametric Programming, and Related Topics, 2nd ed.; Walter de Gruyter: Berlin, Germany, 1995. Gal, T.; Nedoma, J. Multiparametric Linear Programming. Mang. Sci. 1972, 18, 406. Ganesh, N.; Biegler, L. T. A Reduced Hessian Strategy for Sensitivity Analysis of Optimal Flowsheets. AIChE J. 1987, 33, 282. Geoffrion, A. M.; Nauss, R. Parametric and Post-optimality Analysis in Integer Linear Programming. Mang. Sci. 1977, 23, 453. Holm, S.; Klein, D. Three Methods for Post-optimal Analysis in Integer Linear Programming. Math. Prog. Study 1984, 21, 97.

Ierapetritou, M. G.; Pistikopoulos, E. N. Novel Optimization Approach of Stochastic Planning Models. Ind. Eng. Chem. Res. 1994, 33, 1930. Jenkins, L. Parametric Mixed-Integer Programming: An Application to Solid Waste Management. Mang. Sci. 1982, 28, 1270. Jenkins, L. Parametric Methods in Integer Linear Programming. Ann. Oper. Res. 1990, 27, 77. Jenkins, L.; Peters, D. A Computational Comparison of Gomory and Knapsack Cuts. Comput. Oper. Res. 1987, 14, 449. Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A General Algorithm for Short Term Scheduling of Batch OperationssI. MILP Formulation. Comput. Chem. Eng. 1993, 17, 211. Liu, M. L.; Sahinidis, N. V. Optimization in Process Planning under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 4154-4165. Marsten, R. E.; Morin, T. L. Parametric Integer Programming: The Right-Hand Side Case. Ann. Discr. Math. 1977, 1, 375. Nemhauser, G. L.; Wolsey, L. A. Integer and Combinatorial Optimization; Wiley-Interscience: New York, 1988. Ohtake, Y.; Nishida, N. A Branch and Bound Algorithm for 0-1 Parametric Mixed-Integer Programming. Oper. Res. Lett. 1985, 4, 41. Papoulias, S.; Grossmann, I. E. A Structural Optimization Approach in Process Synthesis. Part I: Utility Systems. Comput. Chem. Eng. 1983a, 7, 695. Papoulias, S.; Grossmann, I. E. A Structural Optimization Approach in Process Synthesis. Part II: Heat Recovery Networks. Comput. Chem. Eng. 1983b, 7, 707. Pertsinidis, A. On the parametric optimization of mathematical programs with binary variables and its application in the chemical engineering process synthesis. Ph.D. Dissertation, Carnegie Mellon University, Pittsburgh, PA, 1992. Pistikopoulos, E. N.; Grossmann, I. E. Optimal Retrofit Design for Improving Process Flexibility in Linear Systems. Comput. Chem. Eng. 1988, 7, 719. Raman, R.; Grossmann, I. E. Symbolic Integration of Logic in Mixed-Integer Linear Programming. Comput. Chem. Eng. 1993, 17, 909. Reinhart, H. J.; Rippin, D. W. T. The Design of Flexible Batch Chemical Plants. AIChE Annual Meeting, New Orleans, Nov 1986; Paper 50e. Roodman, G. M. Post-optimality Analysis in zero-one Programming by Implicit Enumeration: The Mixed-Integer Case. Nav. Res. Log. Quart. 1974, 21, 595. Sahinidis, N. V.; Grossmann, I. E.; Fornari, R. E.; Chathrathi, M. Optimization Model for Long-Range Planning in Chemical Industry. Comput. Chem. Eng. 1989, 13, 1049. Seferlis, P.; Hrymak, A. N. Sensitivity Analysis for Chemical Process Optimization. Comput. Chem. Eng. 1996, 20, 1117. Shah, N.; Pantelides, C. C. Optimal Long-Term Campaign Planning and Design of Batch Operations. Ind. Eng. Chem. Res. 1991, 30, 2308. Shah, N.; Pantelides, C. C. Design of Multipurpose Batch Plants with Uncertain Production Requirements. Ind. Eng. Chem. Res. 1992, 31, 1325. Straub, D. A.; Grossmann, I. E. Integrated Stochastic Metric of Flexibility for Systems with Discrete State and Continuous Parameter Uncertainties. Comput. Chem. Eng. 1990, 14, 967. Takamatsu, T.; Hashimoto, I.; Ohno, H. Optimal Design of a Large Complex System from the Viewpoint of Sensitivity Analysis. Ind. Eng. Chem. Process Des. Dev. 1970, 9, 368. Voudouris, V.; Grossmann, I. E. Mixed-Integer Linear Programming Reformulations for Batch Process Design with Discrete Equipment Sizes. Ind. Eng. Chem. Res. 1992, 31, 1314.

Received for review July 29, 1996 Revised manuscript received December 16, 1996 Accepted December 16, 1996X IE960451L X Abstract published in Advance ACS Abstracts, February 1, 1997.