3976
Ind. Eng. Chem. Res. 1999, 38, 3976-3987
Algorithms for the Solution of Multiparametric Mixed-Integer Nonlinear Optimization Problems Vivek Dua and Efstratios N. Pistikopoulos* Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College, London SW7 2BY, U.K.
In this paper we present novel theoretical and algorithmic developments for the solution of mixedinteger optimization problems involving uncertainty, which can be posed as multiparametric mixed-integer optimization models, where uncertainty is described by a set of parameters bounded between lower and upper bounds. In particular, we address convex nonlinear formulations involving (i) 0-1 integer variables and (ii) uncertain parameters appearing linearly and separately and present on the right-hand side of the constraints. The developments reported in this work are based upon decomposition principles where the problem is decomposed into two iteratively converging subproblems: (i) a primal and (ii) a master subproblem, representing valid parametric upper and lower bounds on the final solution, respectively. The primal subproblem is formulated by fixing the integer variables which results in a multiparametric nonlinear programming (mp-NLP) problem, which is solved by outer-approximating the nonlinear functions at a number of points in the space of uncertain parameters to derive linear profiles. The aim of the master subproblem is then to propose another set of integer solutions which are better than the integer solutions that have already been analyzed in the primal subproblem. This can be achieved by (i) introducing cuts, (ii) employing outer-approximation (OA) principles, or (iii) using generalized benders decomposition (GBD) fundamentals. The algorithm terminates when the master problem cannot identify a better integer solution. 1. Introduction Process engineering problems usually involve discrete and continuous decisions. Discrete decisions correspond to, for example, selection and arrangement of process units and continuous decisions may correspond to operating conditions such as flow rate. In a mathematical framework, the former decisions can be incorporated by introducing 0-1 integer variables and the latter can be grouped as continuous variables, which results in a mixed-integer nonlinear programming (MINLP) problem of the following form:1
z ) min dTy + f(x) y,x
s.t. Ey + g(x) e b
(1)
x ∈ R ; y ∈ {0, 1} n
m
where y is a vector of 0-1 binary variables, x is a vector of continuous variables, f is a scalar, continuously differentiable and convex function of x, g is a vector of continuously differentiable and convex functions of x, b, and d are constant vectors, and E is a constant matrix. Mixed-integer problems of the form similar to the one outlined in ref 1 find applications2,3 in many process synthesis problems such as heat-exchanger networks,4-9 distillation sequences,10-16 reactor networks,17-21 and total process flow sheets,22,23 in material design problems such as refrigerant design,24-26 solvent design,27,28 and polymer design29,30 and in scheduling and design of batch plants.31-37 * To whom correspondence should be addressed. Tel.: (44) 2075946620. Fax: (44) 2075946606. E-mail: e.pistikopoulos@ ic.ac.uk.
However, uncertainty is inevitably present in any process modelstypically associated with technical parameters, such as heat-transfer coefficients, and commercial parameters, such as cost of raw materials.38 Solution approaches for process synthesis problems under uncertainty can be broadly classified as stochastic and parametric programming techniques; while the stochastic programming approach provides an optimal solution based upon an expected performance criteria for a given probability distribution of the uncertain parameters, the parametric programming approach gives a complete profile of optimal solutions as a function of uncertain parameters.39 In this work, we will focus on developing solution techniques which are based upon parametric programming fundamentals. Parametric programming has been studied for many classes of mathematical programming problems. For linear models without integer variables, the solution is given by linear parametric profiles which are obtained by extending the simplex algorithm to incorporate uncertain parameters,40,41 and for nonlinear programs without integer variables, the solution is obtained by bounding the nonlinear function by linear profiles by using continuity and convexity properties of the objective function.42-48 The presence of integer variables in the model results is an extra complication due to the discrete nature of the solution space. For linear models involving integer variables, the solution techniques are based upon either relaxing the integrality condition or introducing constraints which are obtained from the solution for fixed-integer variables.49-59 For the case when the model is nonlinear and involves integer variables, the solution approaches have been limited to the case when only a single uncertain parameter is present in the model.57,59,60-64 In this work, we propose
10.1021/ie980792u CCC: $18.00 © 1999 American Chemical Society Published on Web 09/16/1999
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3977
algorithms for the solution of multiparametric mixedinteger nonlinear programming problems which are based upon (i) the systematic characterization of the space of uncertain parameters and (ii) decomposition principles. The outline of the rest of the paper is follows. In section 2, a theoretical framework and solution algorithms for addressing multiparametric MINLP problems are proposed. The solution steps of the algorithm are illustrated with two process synthesis examples in section 3 and the concluding remarks are given in section 4. 2. Multiparametric MINLP ProblemssMathematical Developments
Figure 1. Evaluation of zdiff.
In the presence of uncertain parameters in the model, process engineering problems1 can be recasted into multiparametric mixed-integer nonlinear programs (mpMINLP) of the following form:
g(x) at x* to obtain the following multiparametric LP (mp-LP) problem:
zˇ (θ) ) min dTyj + f(x*) + ∇xf(x*)(x - x*) x
z(θ) ) min dTy + f(x)
s.t. Eyj + g(x*) + ∇xg(x*)(x - x*) e b + Fθ
s.t. Ey + g(x) e b + Fθ
θmin e θ e θmax
y,x
θmin e θ e θmax
(2)
x ∈ Rn; θ ∈ Rs
x ∈ Rn; y ∈ {0, 1}m; θ ∈ Rs where θ is a vector of uncertain parameters and F is a constant matrixsnote that the formulation in (2) is similar to (1) except that it involves a vector of uncertain parameters, θ. The solution of problem (2) is approached by decomposing it into (i) a primal and (ii) a master subproblem, representing valid parametric upper and lower bounds, respectively. While the primal subproblem is obtained by fixing the integer variables, the master subproblem identifies new integer variables. 2.1. Initialization. Define an initial upper bound zˆ (θ) ) ∞ in the entire space of θ. An initial feasible point is obtained by solving the following deterministic MINLP:
z* ) min dTy + f(x) y,x,θ
s.t. Ey + g(x) e b + Fθ θmin e θ e θmax
(3)
x ∈ Rn; y ∈ {0, 1}m where θ is treated as a free variable. If (3) has no feasible solution, (2) has no solution; otherwise, the solution of (3) is given by yj, x*, and θ*. 2.2. Primal SubproblemsMultiparametric NLP (mp-NLP). The integer variables in (2) are fixed at y ) yj, which results in a multiparametric NLP (mp-NLP) problem of the following form:
zˆ (θ) ) min dTyj + f(x) x
s.t. Eyj + g(x) e b + Fθ θmin e θ e θmax
(5)
(4)
x ∈ Rn; θ ∈ Rs It has been shown that zˆ (θ) is a convex function of θ43s this convexity property is used here to approximate zˆ (θ) by linear profiles as follows. An outer-approximation of (4) is created by linearizing the nonlinear terms f(x) and
The solution of (5), zˇ (θ), is given by a linear function of θ valid in a critical region, CR, of θ, where a critical region is defined as the feasible region of θ associated with an optimal basis40,58ssince (5) is a relaxation of (4), zˇ (θ) represents a lower bound to zˆ (θ). A number of such relaxations at different fixed points in θ are then created to obtain tighter linear approximations (within ) to zˆ (θ). The points for approximation are systematically identified by using the convexity of zˆ (θ): the difference zdiff ) zˆ (θ) - zˇ (θ) will be maximum at one of the vertices of the CR, and hence checking at the vertices is a sufficient condition for identifying a point of maximum discrepancy. A graphical interpretation for the single parametric case is depicted in Figure 1 (see Appendix A for details on finding the vertices of a CR and Appendix B for the case when an infeasibility is found at a vertex). Thus, θ is fixed at the vertices of a CR at which zdiff is evaluated and then the point at which zdiff is maximum is identified. At this point another outer-approximation is performed and a new mp-LP problem (5) is formulated and solved to obtain a new lower bound zˇ (θ)2. These two lower bounds are then compared to retain a tighter lower bound, by formulating a constraint of the form zˇ (θ)1 - zˇ (θ)2 e 0 (zˇ (θ)1 and zˇ (θ)2 are the two lower bounds) and by using a comparison procedure as described in Acevedo and Pistikopoulos58 (see Appendix C for an outline of the comparison procedure). This comparison effectively subdivides the CR, as graphically shown in Figure 2 (for the case of a single parameter) where the region CR ) [θmin, θmax] is divided into subregions CR1 ) [θmin, θi] and CR2 ) [θi, θmax]. For each subregion such approximations are continued until the nonlinear solution is approximated by linear profiles within a certain tolerance (i.e., zdiff e ). In the regions, CRi, where zdiff e , the upper bound, zˆ (θ)i, is updated as zˆ (θ)i ) zˇ (θ)i. Note that zˇ (θ)i represents underestimators to the nonlinear solution zˆ (θ)i and can be converted to overestimators by adding the tolerance to zˇ (θ)i. The resulting overestimator, zˇ (θ)i + , will approximate the nonlinear solution within the
3978
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999
solutions that have already been analyzed in CRi. Note that the inequality, dTy + f(x) e zˆ (θ)i, constrains the objective function, zj, to take a value which is lower than the current upper bound zˆ (θ)i and the inequality, ∑j∈Jik ik ik ik yik j - ∑j∈L yj e |J | - 1, corresponds to integer cuts which eliminate the integer solutions which have already been explored in the primal subproblem. Note that (6) is solved for each region, CRi, obtained from the solution of the primal subproblem. The integer solution obtained from the solution of (6) is then returned back to the primal subproblem (4). The algorithm terminates in a region, CRi, when there is no feasible solution to the deterministic master subproblem (6) in that region. 2.3.2. Master SubproblemsOuter-Approximation Formulation. A new vector of integer variables can also be obtained by formulating an alternative master subproblem based upon outer-approximations (OA)65 principles, by constructing outer-approximations at the vertices of the CRi, which results in an mp-MILP master subproblem of the following form:
Figure 2. Comparison of lower bounds.
zj(θ) ) min dTy + µ y,x
v
s.t. µ g f(x ) + ∇xf(xv)(x - xv) Ey + g(xv) + ∇xg(xv)(x - xv) e b + Fθ
∑j∈J yikj - ∑j∈L yikj e |Jik| - 1, ik
ik
(7)
k ) 1, ..., Ki
θ ∈ CRi; x ∈ Rn; y ∈ {0, 1}m Figure 3. Conversion of an underestimator to an overestimator.
given tolerancessee Figure 3 for a graphical interpretation. However, irrespective of whether an underestimator or an overestimator is used as the solution of mpNLP, there might be some error in the solutionssee remark (ix) in section 2.5. The final solution of the primal subproblem is given by (i) a vector of linear parametric upper bounds, zˆ (θ)i, valid in the critical regions, CRi, and the corresponding integer solution, and (ii) a set of infeasible critical regions where zˆ (θ)i ) ∞. 2.3. Master Subproblem. In this section three formulations are proposed for obtaining a new vector of integer variables. While the first formulation corresponds to an MINLP formulation, the second and third result in a multiparametric mixed-integer linear programming (mp-MILP) problem. 2.3.1. Master SubproblemsDeterministic Formulation. The following MINLP master subproblem is formulated and solved to identify the next integer solution:
zj ) min z ) dTy + f(x) y,x,θ
s.t. Ey + g(x) e b + Fθ dTy + f(x) e zˆ (θ)i
∑
ik j∈Jikyj
-
∑
ik j∈Likyj
ik
e |J | - 1,
(6) i
k ) 1, ..., K
where µ is a free variable; the superscript v corresponds to the solutions obtained at the vertices of the CR. Note that problem (7) is solved for each region, CRi, to obtain a vector of integer variables, which are returned back to the primal subproblem. However, if the same integer solution is identified in all the regions, then the primal subproblem can be solved in the complete space of θ instead of solving a number of mp-NLPs in all the regions separately. Also note that the linearizations obtained from the primal subproblem for a fixed integer vector can be included for all the corresponding regions. The solution of the mp-MILP, (7), zj(θ), which represents a lower bound, is given by the parametric profiles and the corresponding integer solutions which are valid in corresponding regions of optimality.58,66,67 If in a region, CRi, (7) has no feasible solution or the lower bound, zj(θ), exceeds the upper bound, zˆ (θ), the algorithm stops in that region; otherwise the new vector y is returned to the primal subproblem. 2.3.3. Master SubproblemsGeneralized Benders Decomposition Formulation. The master subproblem can also be reformulated on the basis of generalized benders decomposition (GBD) principles by using Lagrangian relaxations64 at the vertices of the CRi resulting in an mp-MILP of the following form:
z˜ (θ) ) min µ y,µ
s.t. µ g dTy + f(xv) + λT[Ey + g(xv) - b - Fθ]
θ ∈ CRi; x ∈ Rn; y ∈ {0, 1}m
0 g λT[Ey + g(xvinf) - b - Fθ]
where θ is treated as a variable and θ ∈ CRi indicates that θ is bounded by the set of inequalities which define ik ik ik CRi; Jik ) (j|yik j ) 1) and L ) (j|yj ) 0), and |J | is the cardinality of Jik and Ki is the number of integer
θ ∈ CRi; y ∈ {0, 1}m
(8)
where λ represents the Lagrange multipliers, and the subscript “inf” represents the vertices for which the
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3979
primal subproblem is infeasible where a minimization of infeasibilities subproblem is solved.64 Note that (8) is an mp-MILP, the solution of which is analyzed similar to the solution of (7); however, the formulation in (8) has µ as the only continuous variable. 2.4. Multiparametric MINLP Algorithm. On the basis of the above theoretical developments, the main steps of the algorithm can be summarized as follows. Step 0 (Initialization). Define an initial region of θCR, with the best upper bound zˆ (θ) ) ∞, an initial integer structure yj, a feasible point for the linearization, and a tolerance for the approximation of mp-NLP by linear profiles. Step 1 (Primal Subproblem). For each region with a new integer structure, yj: (a) For a given region of θ create outer-approximation of the mp-NLP problem (4) and solve the resulting mp-LP problem (5) to obtain a set of parametric profiles, zˇ (θ), and the corresponding critical regions CR. (b) Identify the point where the difference between the solution of mp-NLP, zˆ (θ), and the solution of the mp-LP, zˇ (θ) is maximum. If zˆ (θ) - zˇ (θ) g in some region CR, return to step 1a with this CR and the point where the difference is maximum as the new point for outer-approximation. (c) If zˆ (θ) - zˇ (θ) e for some region of θ, update the best upper bound function, zˆ (θ), and the corresponding integer solution. Deterministic Master Subproblem (6) Step 2 (Master Subproblem). For each region CR, formulate and solve a deterministic MINLP master subproblem (6) with (i) θ as a variable bounded in the region CR, (ii) an integer cut, and (iii) a parametric cut z e zˆ (θ). Return to Step 1 with the new integer solutions and the corresponding CRs. Step 3 (Convergence). The algorithm terminates when the solution of the deterministic MINLP is infeasible for each region CR. The final solution is given by the current upper bounds zˆ (θ) in the corresponding CRs. Outer-Approximation Master Subproblem (7) Step 2 (Master Subproblem). For each region CR, formulate and solve mp-MILP master subproblem (7) to obtain parametric lower bounds zj(θ) and integer solutions. Step 3 (Convergence). For each region CR, compare the lower bounds zj(θ) to the upper bounds zˆ (θ) (see Appendix C for a comparison procedure). For the regions where zj(θ) g zˆ (θ) or master subproblem (7) is infeasible, the solution is given by zˆ (θ). For the rest of the CRs with corresponding vectors of new 0-1 variables, y, return to Step 1. Generalized Benders Decomposition Master Subproblem (8) Step 2 (Master Subproblem). For each region CR, formulate and solve mp-MILP master subproblem (8) to obtain parametric lower bounds, zj(θ), and integer solutions. Step 3 (Convergence). For each region CR, compare the lower bounds zj(θ) to the upper bounds zˆ (θ) (see Appendix C for a comparison procedure). For the regions where zj(θ) g zˆ (θ) or master subproblem (8) is infeasible,
the solution is given by zˆ (θ). For the rest of the CRs with corresponding vectors of new 0-1 variables, y, return to Step 1. 2.5. Remarks. (i) The formulation given in (2) also includes linear equalities, as they can be rewritten as two opposing inequalities. For nonlinear equalities, however, we assume that they can be relaxed as convex inequalitiess this relaxation can be achieved in many process engineering problems.63,64 This limitation comes from the fact that for fixed integer variables for the parametric solution to be continuous and convex and hence linear approximations to be valid underestimators, the equalities should be affine.42,43 (ii) The formulation given in (2) also incorporates the case when the uncertain parameters are correlated linearly by Pθ e p, where P is a constant matrix and p is a constant vector, as some of the continuous and binary parts of the inequalities in (2) can be made zero to write such correlations. In such a case the correlations, Pθ e p, are imposed on the CR that is normally obtained and a compact representation of the CR is obtained by removing redundant constraints. (iii) For linear models with a single uncertain parameter, the formulation in (6) reduces to the formulation proposed by Pertsinidis et al.57,59 (iv) For nonlinear models with a single uncertain parameter, the formulation in (7) reduces to the formulation given in Acevedo and Pistikopoulos,63 and Pertsinidis et al.57,59 (v) For linear models with more than one uncertain parameter, the formulation in (6) reduces to the formulation given in Dua and Pistikopoulos.66 (vi) At the begining of each iteration, we have a number of starting regions. As the iterations between the primal and master subproblems proceed, the space of θ gets divided into smaller subspaces due to the solution of mp-NLPs and mp-MILPs at each iteration. This may result in a large number of regions in the final solution which may actually not be required; i.e., we could do with a less number of regions and hence a less number of subproblems in each region. This can be achieved by ignoring the subdivision of the region of θ that results from solving an mp-NLP for a y which is greater than the current upper bound. (vii) The deterministic master subproblem (6), which is an MINLP, converges when no feasible solution is found. This MINLP can be solved either by using a branch-and-bound method or by enforcing a (heuristic) limit on the maximum number of iterations for decomposition-based algorithms to avoid the evaluation of all the integer solution alternatives. (viii) Note that the generalized benders decomposition master subproblem (8) is similar to the one proposed by Papalexandri and Dimkou.64 As the authors mentioned, Lagrangian constraints exclude the already analyzed integer solutions, but integer cuts may be included due to numerical errors in the calculation of Lagrange multipliers. Also note that GBD master subproblem is smaller in size than the OA master subproblem, but it is a relaxation of the OA master subproblem, and hence gives a loose lower bound. For the solution of mp-NLPs, the approach outlined in this work is based upon characterizing the space of θ by critical regions, whereas the approach of Papalexandri and Dimkou64 is based upon discretizing the space of θ. (ix) The solution of the primal subproblem (4) is
3980
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999
Figure 4. Flow sheet for example 1.
obtained as a piecewise linear approximation of the nonlinear solution. This approximation may also result in an approximation of the parametric regions corresponding to two different integer solutions and hence may introduce some error in the identification of the exact hyperplane which separates parametric solutions corresponding to two different integer solutionssthis error is however bounded by ( (see Appendix D for details). (x) In the worst case, the computational requirements for the solution of (5) grow exponentially or faster with the size of the problem and are not bounded above by a polynomial in the size of the problem.68 Since the solution of mp-NLP and mp-MILP subproblems invloves a repetitive solution of such mp-LP problems, the computational requirements are further escalated and increase with the nonlinearity and size of the problem. In this work, strategies for alleviating some of these problems have been presentedsthe problem nevertheless still remains computationally intensive. For mpNLPs, linear approximation seems to be the natural strategy. We have suggested efficient identification of points in the space of θ to carry out linear approximations by using convexity properties of the objective function. The number of linear approximations required is a function of the nonlinearity of the model and hence is specific to the problem under consideration. This number of approximations can however be reduced by increasing the tolerance (), albeit at the cost of loss of accuracy. For the solution of mp-MILPs, we used integer and parametric cuts to avoid a complete enumeration of the space of binary variables.66,67 3. Examples In this section, we present two process synthesis example problems to illustrate the solution steps and the typical form of the solution that is obtained from solving such problems parametrically. 3.1. Example 1. Consider a slight modification of a process synthesis under uncertainty example taken from ref 22 for the manufacture of a product C from a chemical B (see Figure 4). The chemical B can either be purchased from the market or manufactured from a chemical A by two different alternatives. The supply of A and the demand of C, given by θ1 and θ2, respectively, are the uncertain parameters. The mathematical model is given in Table 1 and the solution steps are summarized below. 3.1.1. Initialization. (i) An initial feasible point is obtained by solving the formulation in Table 1 as an MINLP with θ as a variable bounded between its lower and upper bounds. The solution is given by y1 ) 1, y2 ) 0, and y3 ) 1. Initialize an upper bound of zˆ (θ) ) ∞ for the full space of θ. Define a tolerance, ) 0.3%, for approximating the nonlinear
Figure 5. Division of critical regionssstep 1. Table 1. Mathematical Model of Example 1 z(θ) ) minx,y 3.5y1 + y2 + 1.5y3 + B2 + 1.2B3 + 1.8(A2 + A3) + 7BP - 11C1 s.t. -0.9B1 + C21/15 e 0 -A2 + exp(B2) e 1 -A3 + exp(B3/1.2) e 1 -A + A2 + A3 ) 0 -B2 - B3 - BP + B1 ) 0 BP e 1.95 y1 + y2 + y3 e 2 C1 e 20y1, B2 e 20y2, B3 e 20y3, A2 e 20y2, A3 e 20y3 A ) θ1; C1 ) θ2 0.50 e θ1 e 0.75 5.50 e θ2 e 6.00
solution by linear profiles. The problem is now divided into two subproblems, (i) a primal subproblem and (ii) a master subproblemsthe solution of the former subproblem, zˆ (θ), represents an upper bound and that of the latter, zj(θ), represents a lower bound. The algorithm converges by iterating between these two subproblems. 3.1.2. Primal Subproblem. (ii) Fix y1 ) 1, y2 ) 0, and y3 ) 1 in Table 1 and solve the resulting mp-NLP with θ as a variable. At the solution thus obtained, linearize the nonlinear constraints. The linearized constraints are given by the following: -0.9B1 + 0.793C1 e 2.359, -A2 + B2 e 0, and -A3 + 1.458B3 e 0.229. (iii) Replacing the nonlinear constraints in Table 1 by these linear constraints, for fixed y, we obtain a multiparametric LP (mp-LP), the solution of which is given by the following linear parametric solution: zˇ (θ)1 ) -2.178θ1 - 4.832θ2 - 14.26; CR1, -0.686θ1 + 0.881θ2 e 4.728, 0.5 e θ1 e 0.75, 5.5 e θ2. Note that CR1, which represents the critical region for zˇ 1(θ), covers only a partial space of θ; the space left out, CR2, is an infeasible region for the current integer configuration. CR1 and CR2 are graphically shown in Figure 5. Since CR1 is an outer-approximation of the nonlinear feasible critical region, this outer-approximation is converted into an inner-approximation by subtracting a small number, 0.005, from the right-hand side of the less than inequality which separates CR1 and CR2. The linear parametric solution, zˇ (θ)1, represents a lower bound to the solution of the mp-NLP subproblem, zˆ (θ). (iv) A point is identified where the difference, zdiff between the solution to mp-NLP, zˆ (θ), and the solution to mp-LP, zˇ (θ), is maximum. This can be achieved by evaluating zˆ (θ) and zˇ (θ) at the vertices of CR1 as shown in Table 2.
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3981
Figure 6. Division of critical regionssstep 2.
Figure 7. Division of critical regionssstep 3.
Table 2. Evaluation at VerticessStep 1 vertex (θ1,θ2)
zˆ (θ)
zˇ (θ)1
(0.50, 5.5) (0.50, 5.75) (0.75, 5.5) (0.75, 5.949)
-41.74 infeasible -42.36 -44.63
-41.92 ** -42.47 -44.64
Table 3. Evaluation at VerticessStep 2 zdiff
(%)
vertex (θ1,θ2)
zˆ (θ)
zˇ (θ)1
zdiff
(%)
0.19 ** 0.11 0.01
0.455 ** 0.26 0.02
(0.75, 5.575) (0.75, 5.944) (0.653, 5.87)
-42.75 -44.61 -44.01
-42.83 -44.62 -44.03
0.08 0.01 0.02
0.189 0.022 0.045
Table 4. Evaluation at VerticessStep 3
(v) Since the subproblem is infeasible at θ ) (0.5, 5.75), a tight feasible point, θf, is obtained from the solution of an infeasibility subproblem (9). The solution is given by θf ) (0.51, 5.742). Nonlinear constraints are linearized at θf ) (0.51, 5.742) to obtain the following linear constraints: -0.9B1 + 0.766C1 e 2.2, -A2 + B2 e 0, and -A3 + 1.258B3 e 0.112. (vi) By again replacing the nonlinear constraints with linear constraints and solving the resulting mp-LP, we obtain the following parametric solution: zˇ (θ)3 ) -2.8178θ1 - 5.0422θ2 - 12.60; CR3, -0.686θ1 + 0.881θ2 e 4.723, -0.796θ1 + 0.851θ2 e 4.48, 0.5 e θ1 e 0.75, 5.5 e θ2. The critical regions are shown in Figure 6. (vii) The two parametric solutions zˇ 1(θ) and zˇ 3(θ) are compared to subdivide the space of θ. The resulting parametric solutions and their corresponding CRs are given by the following: zˇ 1(θ) ) -2.178θ1 - 4.8322θ2 14.26; CR1, -0.686θ1 + 0.881θ2 e 4.723, 0.63θ1 + 0.21θ2 g 1.65, θ1 e 0.75, 5.5 e θ2. zˇ 3(θ) ) -2.8178θ1 - 5.0422θ2 - 12.60; CR3, -0.796θ1 + 0.851θ2 e 4.48, 0.63θ1 + 0.21θ2 e 1.65, 0.5 e θ1 e 0.75, 5.5 e θ2. The graphical interpretation of CRs is shown in Figure 7. Again, 0.005 is subtracted from the right-hand side of the inequality which separates CR3 and CR4. (viii) We again identify the point of maximum difference, by comparing at the vertices, for each CR as shown in Tables 3 and 4. Note that both zˇ (θ)1 and zˇ (θ)3 satisfy the tolerance criteria, i.e., e 0.3%, in their respective CRs and the max{zdiff} ) 0.09. These linear approximations to the nonlinear solutions are then converted to overestimators by adding 0.09 to them. Thus, the final solution of mp-NLP is given by the following: zˇ 1(θ) ) -2.178θ1 - 4.8322θ2 - 14.17; CR1, -0.686θ1 + 0.881θ2 e 4.723, 0.63θ1 + 0.21θ2 g 1.65, θ1 e 0.75, 5.5 e θ2. zˇ (θ)2 ) ∞, CR2, -0.686θ1 + 0.881θ2 g 4.723, 0.5 e θ1 e 0.75, θ2 e 6.0. zˇ 3(θ) ) -2.8178θ1 - 5.0422θ2 - 12.51; CR3, -0.796θ1 + 0.851θ2 e 4.48, 0.63θ1 + 0.21θ2 e 1.65, 0.5 e θ1 e 0.75, 5.5 e θ2. zˇ (θ)4 ) ∞; CR4, -0.686θ1 + 0.881θ2 e 4.723, -0.796θ1 + 0.851θ2 g 4.48. This completes the first iteration of the primal subproblem.
vertex (θ1,θ2)
zˆ (θ)
zˇ (θ)3
zdiff
(%)
(0.50, 5.5) (0.50, 5.72627) (0.649, 5.866) (0.653, 5.87) (0.75, 5.575) (0.75, 5.5)
-41.74 -42.91 -44.00 -44.01 -42.75 -42.36
-41.75 -42.91 -44.01 -44.03 -42.83 -42.45
0.01 0.00 0.01 0.02 0.08 0.09
0.02396 0.00 0.023 0.045 0.189 0.212
Table 5. Deterministic Master Subproblem for Example 1 for θ ∈ CR1 zj ) miny,x,θ 3.5y1 + y2 + 1.5y3 + B2 + 1.2B3 + 1.8(A2 + A3) + 7BP - 11C1 s.t. -0.9B1 + C21/15 e 0 -A2 + exp(B2) e 1 -A3 + exp(B3/1.2) e 1 -A + A2 + A3 ) 0 -B2 - B3 - BP + B1 ) 0 BP e 1.95 y1 + y2 + y3 e 2 y1 - y2 + y3 e 1 z e -2.178θ1 - 4.8322θ2 - 14.17 C1 e 20y1, B2 e 20y2, B3 e 20y3, A2 e 20y2, A3 e 20y3 A ) θ1; C1 ) θ2 -0.686θ1 + 0.881θ2 e 4.723 0.63θ1 + 0.21θ2 g 1.65 θ1 e 0.75, 5.5 e θ2
3.1.3. Master Subproblem. (ix) For each region, CRi, i ) 1,2,3,4, obtained in the primal subproblem, the deterministic MINLP master subproblem is formulated and solved to identify the next integer solutionssee Table 5 for a sample master problem for i ) 1. The integer solution of this problem for i ) 1,3, which is given by y1 ) 1, y2 ) 1, and y3 ) 0, is recycled back to the primal subproblems for CR1 and CR3. In the rest of the regions, the master problem is infeasible and therefore the current solutions in these regions represent the final solution. The algorithm terminates when there is no feasible solution to (5) in each CR and final solution is given by current upper bounds zˆ (θ) in the corresponding CRs. (x) Alternative master subproblems corresponding to OA and GBD formulations can also be derived. The algorithm for this case terminates when the master subproblem has no feasible solution or the lower bound
3982
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999
Figure 8. Parametric solution of example 1. Table 6. Computational Requirements of Example 1 θ ) (0.50, det. formulation OA formulation GBD formulation
mp-LPsa
NLPsa
MILPsa
ITNS
total time
eqv. MINLPs
4(0.238) 8(8.274) 14(4.239)
2(0.46) 56(9.214) 22(3.240) 37(5.671)
2(0.11) 3(0.332) 17(1.302) 27(1.442)
2 3 3 5
0.57 9.784 12.816 11.352
1 17.17 22.48 19.92
5.50)b
a The figures in parentheses are the CPU times in seconds on a Sun Sparc10-51 workstation. b GAMS/DICOPT was used to solve this case.
exceeds the upper bound. The final solution is summarized in Figure 8, the computational requirements are given in Table 6, and the corresponding CRs are graphically shown in Figure 9. 3.2. Example 2. The second example is concerned with obtaining the optimal objective function as a function of uncertainty in the demand of products P1, P2, and P3. Three raw materials, A, B, and C, and one intermediate, I, are available. For each raw material there are two possible alternatives for producing intermediates which are then converted to the final productss see Figure 10. The mathematical model of this process scheme involves 31 continuous variables, 6 binary variables, and 3 uncertain parameters and is described by 6 nonlinear and 37 linear equations/inequalitiess see Table 7. This problem was solved using 3 master subproblem formulationssdeterministic, OA, and GBD. The sensitivity information was obtained by using LPs/ NLPs at fixed θ by using GAMS;69 and a maximum number of 40 major iterations was imposed for deterministic formulation while using DICOPT++70 for solv-
Figure 9. Division of critical regionssfinal solution.
ing MINLPs. While the parametric solution (for ) 0.085) is given in Table 8, the computational require-
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3983 Table 7. Mathematical Model of Example 2 6 11 11 z(θ) ) minx,y ∑i)1 FCiyi + ∑i)7 FCi + ∑i)1 OCiIS2i + 1.2A + 1.5B + 1.8C - 50P1 - 60P2 - 68P3 s.t. -ISi + KCi exp(OSi/PCi) e KCi, i ) {1, ..., 6) -PCiISi + OSi ) 0, i ) {7, ..., 11} ISi - MIiyi e 0 Ai ) IS1 + IS2 A0 ) OS1 + OS2 Bi ) IS3 + IS4 B0 ) OS3 + OS4 Ci ) IS5 + IS6 C0 ) OS5 + OS6 IS7 ) A0 + I OS7 ) D + IS9 IS8 ) D + B0 OS8 ) E + IS10 IS11 ) E + C0 Ai e 1.5; Bi e 1.6; Ci e 1.7 OS9 ) P1 ) θ1; OS10 ) P2 ) θ2; OS11 ) P3 ) θ3 0.80 e θ1 e 0.85 0.75 e θ2 e 0.80 1.05 e θ3 e 1.10 MIi ) {4,3,3.5,5.5,7.5,6,8,6,5,5,5} PCi ) {1.8,2.0,2.5,2.2,2.8,3.0,0.65,0.9,0.72,0.68,0.71} KCi ) {2.0,2.1,2.6,2.1,2.5,3.0} FCi ) {1.9,2.0,3.0,3.2,3.7,3.9,1.1,1.0,5.0,5.3,6.2} OCi ) {1.1,1.3,1.6,1.5,2.1,1.9,1.7,1.0,3.0,3.7,4.1}
ments are given in Table 9. The number of major iterations required for OA formulation is much less than that for GBD formulation; however, the solution time per major iteration for GBD is less than that for OA. 4. Concluding Remarks In this work we have presented algorithms for solving multiparametric mixed-integer nonlinear programs. The solution approaches are based upon decomposing the problem into two converging subproblemss(i) a primal and (ii) a master subproblem, which represent valid upper and lower bounds, respectively. The primal subproblem, which is formulated by fixing the integer variables, is solved by approximating the nonlinear solution with piecewise linear parametric profiles. The points where the discrepancy between linear approximation and the nonlinear solution is maximum are then
Figure 10. Flow sheet for example 2.
depicted by using the convexity property of the objective function. At these points new linear approximations are obtained, which are then used to tighten the linear approximations by rejecting the loose profiles and retaining the tight ones. The solution of a primal subproblem is given by a set of linear parametric profiles which are valid in the corresponding regions of optimality. The solution obtained in the primal subproblem is then used to formulate a master subproblem to obtain a next set of integer variables. Three master subproblem formulations were proposed in this worksthe first one is based upon introducing cuts, the second one is based upon outer-approximating the nonlinear functions at the points where linear approximations were carried out in the primal subproblem, and the third formulation relies on using the dual information obtained from the solution of the primal subproblem. The new vector of integer variables obtained in the master subproblem is then returned back to the primal subproblem, and the algorithm terminates when either no integer solution is obtained in the master subproblem or the solution of the master subproblem is greater than the solution of
Table 8. Parametric Solution of Example 2 y
z(θ)
critical region
111111
-12.01θ1 - 34.684θ2 - 41.145θ3 - 1.92
1.355θ1 + 1.561θ2 + 1.219θ3 g 3.63845
111111
-13.365θ1 - 36.245θ2 - 42.364θ3 + 1.71845
1.355θ1 + 2.88θ2 + 2.351θ3 e 5.845 5.692θ2 + 8.151θ3 g 13.24055
111111
-13.365θ1 - 36.245θ2 - 42.364θ3 + 1.71845
1.355θ1 + 1.561θ2 + 1.219θ3 e 3.63845 1.355θ1 + 8.572θ2 + 10.502θ3 g 19.08593
111110
-13.365θ1 - 30.553θ2 - 34.213θ3 - 11.5221
1.355θ1 + 2.88θ2 + 2.351θ3 e 5.845 0.312θ1 + 17.036θ2 - 1.219θ3 g 11.94844
111110
-12.01θ1 - 27.673θ2 - 31.862θ3 - 17.36748
1.355θ1 + 1.561θ2 + 1.219θ3 e 3.63845 1.355θ1 + 2.88θ2 + 2.351θ3 g 5.845 1.355θ1 + 8.572θ2 + 10.502θ3 e 19.08593 -0.158θ1 + 18.105θ2 - 10.502θ3 g 2.71712
110111
-13.067θ1 - 27.122θ2 - 41.145θ3 - 7.24786
1.355θ1 + 2.88θ2 + 2.351θ3 e 5.845 0.312θ1 + 17.026θ2 - 1.219θ3 e 11.94844 5.692θ2 + 8.151θ3 e 13.24055
110111
-12.451θ1 - 27.122θ2 - 41.203θ3 - 7.67177
1.355θ1 + 1.561θ2 + 1.219θ3 e 3.63845 1.355θ1 + 2.88θ2 + 2.351θ3 g 5.845 1.355θ1 + 8.572θ2 + 10.502θ3 e 19.08593 0.283θ1 + 17.554θ2 - 1.161θ3 e 12.418
110111
-12.168θ1 - 9.568θ2 - 42.364θ3 - 20.0846
1.355θ1 + 1.561θ2 + 1.219θ3 e 3.63845 1.355θ1 + 2.88θ2 + 2.351θ3 g 5.845 1.355θ1 + 8.572θ2 + 10.502θ3 e 19.08593 -0.158θ1 + 18.105θ2 - 10.502θ3 g 2.71712
3984
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999
Table 9. Computational Requirements of Example 2 θ ) (0.85, 0.8, det. formulation OA formulation GBD formulation
LPsa
NLPsa
MILPsa
ITNS
total time
eqv. MINLPs
215(83.162) 633(55.811)
19(4.14) 1619(407.75) 331(58.49) 725(126.928)
19(4.99) 1526(273.52) 52(48.49) 130(111.51)
19 3 17 45
9.13 681.270 190.142 294.249
1 74.62 20.83 32.23
1.1)b
a The figures in parentheses are the CPU times in seconds on a Sun Sparc10-51 workstation. b GAMS/DICOPT was used to solve this case.
the primal subproblem, i.e., when the lower bound crosses the upper bound. A graphical interpretation of main concepts was also provided. The solution algorithms were tested on two process synthesis problems and the solution steps with their graphical interpretation were also presented. The key advantage of this approach is that a complete map of the optimal objective function and operating conditions in the space of uncertain parameters can be obtained. While the computational complexity of these approaches is expected to increase with an increase in the dimensionality of the problem, there are some avenues that can be explored for reducing the computational time; for example, NLPs at the vertices can be solved independently in a parallel computing environment. While theory for solving mp-MINLP problems has been presented in this work, the current research efforts mainly focus on implementation of these algorithms, their application to larger models, and their extension for the case of nonconvex models. Acknowledgment Financial support from EPSRC (IRC Grant) is gratefully acknowledged.
Figure 11. Evaluation of max{zdiff}.
Appendix A. Vertices of a CR For the case when a single uncertain parameter is present in the model, the vertices are given by the end points of the interval of θ, whereas for the multiparametric case the vertices of a CR can be determined iteratively by systematically fixing some of the inequalities, which describe the CR, as equalities, and then solving the resulting system of linear equations. The number of linear systems of equations that has to be solved is given by tCs, where t is the number of inequalities which define the parametric region and s is the number of uncertain parameters. However, the NLPs are solved only at those vertices which satisfy all the t inequalities, which are usually less than tCs. A flowchart describing this procedure is given in Figure 11 and a graphical interpretation of feasible vertices is given in Figure 12 where, although the total number of possible vertices is 8C2 ) 28, only 8 vertices are feasible. Appendix B. Infeasibility If an infeasibility is found at any of the vertices, an infeasibility subproblem is formulated to determine a tight feasible point, θf, as follows:
min y,θ,δ
∑δ2
s.t. Eyj + g(x) e b + Fθ δ ) θ - θvinf x ∈ X; θ ∈ CRi
(9)
Figure 12. Feasible vertices of a parametric region.
where, θvinf corresponds to an infeasible vertex; δ is a vector of free variables. Another outer-approximation is performed at θf, the solution of (9), and the resulting mp-LP is solved to obtain a new lower bound, zˇ (θ)2f , as shown in Figure 13 for the single parametric case. Note that the formulation in (9) minimizes the distance between the infeasible vertex and a feasible point as shown in Figure 14, where the distance AP is minimized and a feasible point P is determined at which the model is outer-approximated and the resulting mp-LP is solved. The outerapproximation at P is then converted to an innerapproximation by subtracting from the right-hand side of the less than inequality describing the outer-approximation (or adding for the case when there is the greater than inequality) as shown in Figure 14. In this case (Figure 13), the region CR ) [θmin, θmax] is divided into subregions CR1 ) [θmin, θf], CR2 ) [θf,
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3985
Figure 15. Critical regions, CR1 and CR2. Figure 13. Outer-approximation of zˆ (θ) at a feasible point, θf.
Case 3: If the new constraint is nonredundant, then
z(θ)1 e z(θ)2, ∀ θ ∈ ∆{C B Rint, z(θ)1 - z(θ)2 e 0} z(θ)1 g z(θ)2, ∀ θ ∈ ∆{C B Rint, z(θ)1 - z(θ)2 g 0} where ∆ is an operator which removes redundant Bint represents the set of constraints constraints and CR which define CRint. Appendix D. Bounds on Error due to Linear Approximations Figure 14. Outer-approximation of a feasible region at a feasible point, θf.
θi], and CR3 ) [θi, θmax]. It may be noted that when no integer variables are present in the model, the infeasible region can be eliminated from further analysis. In the presence of integer variables, it is important to define an infeasible region, which is obtained from the solution of a primal subproblem, to explore for another integer solution. For a procedure for defining a region which has been excluded from the initial region, see Appendix E. Appendix C. Comparison of Two Parametric Solutions Acevedo and Pistikopoulos58 proposed an approach for comparing two parametric solutions, z(θ)1 and z(θ)2, which are valid in the critical regions CR1 and CR2, respectively. Their approach, which consists of two steps, is briefly described here. The first step is to define a region, CRint ) CR1 ∩ CR2, where both the parametric solutions are valid. CRint can be defined by removing all the redundant constraints from the set of inequalities which define CR1 and CR2sfor a procedure to identify redundant constraints, see Gal.41 In the second step, check if CRint ) L, then z(θ)1 and z(θ)2 are the solutions in CR1 and CR2, respectively; otherwise a new constraint, z(θ)1 e z(θ)2, is formulated and a constraint redundancy check is made for the new constraint in CRint. This constraint redundancy test results in three cases which are analyzed as follows. Case 1: e z(θ)2, Case 2: g z(θ)2,
If the new constraint is redundant, then z(θ)1 ∀ θ ∈ CRint. If the new constraint is infeasible, then z(θ)1 ∀ θ ∈ CRint.
Let the exact nonlinear solutions corresponding to two different integer solutions be denoted by zˆ (θ)1 and zˆ (θ)2. The exact hyperplane which separates two integer solutions is therefore given by zˆ (θ)1 - zˆ (θ)2 e 0. In this section we show that the error that may occur in obtaining this separating hyperplane, which we define by γ, where zˆ (θ)1 - zˆ (θ)2 e γ, due to using linear approximations is bounded by (. Let zˆ (θ)1 - γ(θ)1 and zˆ (θ)2 - γ(θ)2 denote linear underestimators for zˆ (θ)1 and zˆ (θ)2, which have been obtained for a tolerance less than or equal to . Therefore γ ) γ(θ)1 - γ(θ)2 and the maximum error, γmax, is given by
γmax ) max|γ(θ)1 - γ(θ)2| θ
s.t. 0 e γ(θ)1 e 0 e γ(θ)2 e
(10)
g0 θmin e θ e θmax . The solution of (10) is given by γmax ) ||, which implies - e γ e . The same result for overestimators can also be similarly obtained. Appendix E. Definition of Infeasible Region Given an initial region, CR1, and a feasible region, CR2, such that CR2 ⊆ CR1, a procedure is described in this section to define the infeasible region, CRinf) CR1 - CR2. For the sake of simplifying the explanation of the procedure, consider the case when only two parameters, θ1 and θ2, are present (see Figure 15), where CR1 min is defined by the inequalities: {θmin e θ1 e θmax e 1 1 , θ2 max θ2 e θ2 } and CR2 is defined by the inequalities {C1 e
3986
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999
Figure 16. Division of critical regionssstep 1. Table 10. Definition of Infeasible Region region CR3 CR4 CR5
inequalities θmin 1
C1 g 0, e θ1, θ2 e θmax 2 max C1 e 0, C2 g 0, θ1 e θmax 1 , θ2 e θ 2 min C1 e 0, C2 e 0, C3 g 0, θmin e θ1 e θmax e θ2 1 1 , θ2
Figure 17. Division of critical regionssinfeasible regions.
0, C2 e 0, C3 e 0} where C1, C2, and C3 are linear in θ. The procedure consists of considering one by one the inequalities which define CR2. Considering, for example, the inequality C1 e 0, an infeasible region is given by e θ1, θ2 e θmax CR3: C1 g 0, θmin 1 2 , which is obtained by reversing the sign of inequality C1 e 0 and removing redundant constraints in CR1 (see Figure 16). Thus, by considering the rest of the inequalities, the complete infeasible region is given by CRinf ) {CR3 ∪ CR4 ∪ CR5}, where CR3, CR4, and CR5 are given in Table 10 and are graphically depicted in Figure 17. Literature Cited (1) Grossmann, I. E. MINLP optimization strategies and algorithms for process synthesis. In Proceedings of the 3rd. International Conference on Foundations of Computer-Aided Process Design; Siirola, J. J., Grossmann, I. E., Stephanopoulos, G., Eds.; CACHE, Elsevier: New York, 1990. (2) Floudas, C. A. Nonlinear and Mixed-Integer Optimization. Oxford University Press: New York, 1995. (3) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall: Englewood Cliffs, NJ, 1997. (4) Papoulias, S.; Grossmann, I. E. A structural optimization approach in process synthesis. Part I: Utility systems. Comput. Chem. Eng. 1983, 7 (6), 695.
(5) Cerda´, J.; Westerberg, A. W. Synthesis of heat exchange networks having restricted matches using transportation problem formulations. Chem. Eng. Sci. 1983, 38, 1723. (6) Floudas, C. A.; Ciric, A. R. Strategies for overcoming uncertainties in heat exchanger network synthesis. Comput. Chem. Eng. 1989, 13 (10), 1133. (7) Yee, T. F.; Grossmann, I. E.; Kravanja, Z. Simultaneous optimization models for heat integration I: Energy and area targeting. Comput. Chem. Eng. 1990, 14 (10), 1151. (8) Yee, T. F.; Grossmann, I. E.; Kravanja, Z. Simultaneous optimization models for heat integration III: Synthesis of heat exchanger networks. Comput. Chem. Eng. 1990, 14 (10), 1185. (9) Yee, T. F.; Grossmann, I. E. Simultaneous optimization models for heat integration II: Synthesis of heat exchanger networks. Comput. Chem. Eng. 1990, 14 (10), 1165. (10) Andrecovich, M. J.; Westerberg, A. W. A simple synthesis method based on utility bounding for heat-integrated distillation sequences. AIChE J. 1985, 31 (3), 363. (11) Floudas, C. A. Separation synthesis of multicomponent feed streams into multicomponent product streams. AIChE J. 1987, 33 540. (12) Wehe, R. R.; Westerberg, A. W. An algorithmic procedure for the synthesis of distillation sequences with bypass. Comput. Chem. Eng. 1987, 11 (6), 619. (13) Floudas, C. A.; Paules, G. E., IV. A mixed-integer nonlinear programming formulation for the synthesis of heat-integrated distillation sequences. Comput. Chem. Eng. 1989, 13 (6), 531. (14) Aggarwal, A.; Floudas, C. A. Synthesis of general distillation sequences-nonsharp separations. Comput. Chem. Eng. 1990, 14 (6), 631. (15) Viswanathan, J.; Grossmann, I. E. An alternate MINLP model for finding the number of trays required for a specified separation objective. Comput. Chem. Eng. 1993, 17, 949. (16) Viswanathan, J.; Grossmann, I. E. Optimal feed locations and number of trays for distillation columns with multiple feeds. Ind. Eng. Chem. Res. 1993, 32, 2942. (17) Achenie, L. K.; Biegler, L. T. A superstructured based approach to chemical reactor networks synthesis. Comput. Chem. Eng. 1990, 14 (1), 23. (18) Kokossis, A. C.; Floudas, C. A. Optimization of complex reactor networks-I: Isothermal operation. Chem. Eng. Sci. 1990, 45 (3), 595. (19) Kokossis, A. C.; Floudas, C. A. Optimization of complex reactor networks-II: Nonisothermal operation. Chem. Eng. Sci. 1994, 49 (7), 1037. (20) Balakrishna, S.; Biegler, L. T. Targeting strategies for the synthesis and energy integration of nonisothermal reactor networks. Ind. Eng. Chem. Res. 1992, 31 (9), 2152. (21) Lakshmanan, A.; Biegler, L. T. Synthesis of optimal chemical reactor networks. Ind. Eng. Chem. Res. 1996, 35, 1344. (22) Kocis, G. R.; Grossmann, I. E. Relaxation strategy for the structural optimization for process flowsheets. Ind. Eng. Chem. Res. 1987, 26 (9). 1869. (23) Kravanja, Z.; Grossmann, I. E. PROSYN-an MINLP process synthesizer. Comp. Chem. Eng. 1990, 14, 1363. (24) Churi, N.; Achenie, L. E. K. Novel mathematical programming model for computer aided molecular design. Ind. Eng. Chem. Res. 1996, 35, 3788. (25) Churi, N.; Achenie, L. E. K. The optimal design of refrigerant mixtures for a two-evaporator refrigerant system. Comput. Chem. Eng. 1997, 21 (S), 349. (26) Duvedi, A. P.; Achenie, L. E. K. Designing environmentally safe refrigerants using mathematical programming. Chem. Eng. Sci. 1996, 51, 3727. (27) Macchietto, S.; Odele, O.; Omatsone, O. Design of optimal solvents for liquid-liquid-extraction and gas-absorption processes. Chem. Eng. Res. Des. 1990, 68, 429. (28) Odele, O.; Macchietto, S. Computer aided molecular design: A novel method for optimal solvent selection. Fluid Phase Equilibr. 1993, 82, 47. (29) Vaidyanathan, R.; El-Halwagi, M. Computer-aided synthesis of polymers and blends with target properties. Ind. Eng. Chem. Res. 1996, 35, 627. (30) Maranas, C. D. Optimal computer-aided molecular design: A polymer design case study. Ind. Eng. Chem. Res. 1996, 35, 3403.
Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3987 (31) Sahinidis, N. V.; Grossmann, I. E. MINLP model for cyclic multiproduct scheduling on continuous parallel lines. Comput. Chem. Eng. 1991, 15, 85. (32) Shah, N.; Pantelides, C. C. Optimal long-term campaign planning and design of batch operations. Ind. Eng. Chem. Res. 1991, 30, 2308. (33) Voudouris, V.; Grossmann, I. E. Mixed-integer linear programming reformulations for batch process design with discrete equipment sizes. Ind. Eng. Chem. Res. 1992, 31 (5), 1314. (34) Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A general algorithm for short term schedulling of batch operations-I. MILP formulation. Comput. Chem. Eng. 1993, 17 (2), 211. (35) Rippin, D. W. T. Batch process systems engineering: A retrospective and prospective review. Comput. Chem. Eng. 1993, 17 (S), 1. (36) Bassett, M. H.; Pekny, J. F.; Reklaitis, G. V. Decomposition techniques for the solution of large-scale scheduling problems. AIChE J. 1996, 42 (12), 3373. (37) Ierapetritou, M. G. Floudas, C. A. Effective continuoustime formulation for short-term scheduling. 1. multiple batch processes. Ind. Eng. Chem. Res. 1998, 37, 4341. (38) Pistikopoulos, E. N. Uncertainty in process design and operations. Comput. Chem. Eng. 1995, 19, (S) 553. (39) Pistikopoulos, E. N. Parametric and stochastic programming algorithms for process synthesis, design and optimization under uncertainty. Presented at Aspen World, Boston, MA, 1997. (40) Gal, T.; Nedoma, J. Multiparametric linear programming. Manage. Sci. 1972, 18, 406. (41) Gal, T. Postoptimal Analyses, Parametric Programming, and Related Topics; de Gruyter: New York, 1995. (42) Fiacco, A. V. Introduction to Sensitivity and Stability Analysis in Nonlinear Programming; Academic Press: New York, 1983. (43) Fiacco, A. V.; Kyparisis, J. Convexity and concavity properties of the optimal value function in parametric nonlinear programming. J. Optim. Theory Appl. 1986, 48, 95. (44) Ganesh, N.; Biegler, L. T. A reduced hessian strategy for sensitivity analysis of optimal flowsheets. AIChE J. 1987, 33, 282. (45) Fiacco, A. V.; Kyparisis, J. Computable bounds on parametric solutions of convex problems. Math. Program. 1988, 40, 213. (46) Fiacco, A. V. Global multiparametric optimal value bounds and solution estimates for separable parametric programs. Ann. Oper. Res. 1990, 27, 381. (47) Wolbert, D.; Joulia, X.; Koehret, B.; Biegler, L. T. A reduced hessian strategy for sensitivity analysis of optimal flowsheets. Comput. Chem. Eng. 1994, 18, 1083. (48) Seferlis, P.; Hrymak, A. N. Sensitivity analysis for chemical process optimization. Comput. Chem. Eng. 1996, 20 (10), 1117. (49) Roodman, G. M. Postoptimality analysis in zero-one programming by implicit enumeration. Naval Res. Log. Quart. 1972, 19, 435. (50) Roodman, G. M. Postoptimality analysis in zero-one programming by implicit enumeration: the mixed-integer case. Naval Res. Log. Quart. 1974, 21, 595. (51) Piper, C. J.; Zoltners, A. A. Some easy postoptimality analysis for zero-one programming. Manage. Sci. 1976, 22 (7), 759. (52) Geoffrion, A. M.; Nauss, R. Parametric and postoptimality analysis in integer linear programming. Manage. Sci. 1977, 23 (5), 453.
(53) Klein, D.; Holm, S. Integer programming postoptimal analysis with cutting planes. Manage. Sci. 1979, 25 (1), 64. (54) Holm, S.; Klein, D. Three methods for postoptimal analysis in integer linear programming. Math. Prog. St. 1984, 21, 97. (55) Marsten, R. E.; Morin, T. L. Parametric integer programming: the right-hand side case. Ann. Discr. Math. 1977, 1, 375. (56) Ohtake, Y.; Nishida, N. A branch-and-bound algorithm for 0-1 parametric mixed-integer programming. Operat. Res. Lett. 1985, 4 (1), 41. (57) Pertsinidis, A. On the parametric optimization of mathematical programs with binary variables and its application in the chemical engineering process synthesis. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA, 1992. (58) Acevedo, J.; Pistikopoulos, E. N. A multiparametric programming approach for linear process engineering problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717. (59) Pertsinidis, A.; Grossmann, I. E.; McRae, G. J. Parametric optimization of MILP programs and a framework for the parametric optimization of MINLPs. Comput. Chem. Eng. 1998, 22 S205. (60) McBride, R. D.; Yorkmark, J. S. Finding all the solutions for a class of parametric quadratic integer programming problems. Manage. Sci. 1980, 26, 784. (61) Cooper, M. W. Postoptimality analysis in nonlinear integer programming: the right-hand side case. Naval Res. Log. Quart. 1981, 28, 301. (62) Skorin-Kapov, J.; Granot, F. Nonlinear integer programming: sensitivity analysis for branch-and-bound. Operat. Res. Lett. 1987, 6 (6), 269. (63) Acevedo, J.; Pistikopoulos, E. N. A parametric MINLP algorithm for process synthesis problems under uncertainty. Ind. Eng. Chem. Res. 1996, 35 (1), 147. (64) Papalexandri, K. P.; Dimkou, T. I. A parametric mixedinteger optimization algorithm for multiobjective engineering problems involving discrete decisions. Ind. Eng. Chem. Res. 1998, 37 (5), 1866. (65) Dura´n, M. A.; Grossmann, I. E. A mixed-integer nonlinear programming algorithm for process systems synthesis. AIChE J. 1986, 32 (4), 592. (66) Dua, V.; Pistikopoulos, E. N. An algorithm for the solution of multiparametric mixed integer linear programming problems. Ann. Operat. Res. 1999, accpeted. (67) Dua, V.; Pistikopoulos, E. N. Optimization techniques for process synthesis and material design under uncertainty. Trans. IChemE 1998, 76, 408. (68) Murty, K. Computational complexity of parametric linear programming. Math. Program. 1980, 19, 213. (69) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A Users Guide, Release 2.25; Scientific Press: Redwood City, CA, 1996. (70) Viswanathan, J.; Grossmann, I. E. A combined penalty function and outer-aproximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14 (7), 769.
Received for review December 18, 1998 Revised manuscript received July 1, 1999 Accepted July 9, 1999 IE980792U
