Stochastic Programming Approach for Mixed

Based on the stochastic programming algorithms presented in Acevedo and ..... CPU times reported were obtained on a Sun SPARC Station 10-51. Top of Pa...
0 downloads 0 Views 176KB Size
Download PDF
2262

Ind. Eng. Chem. Res. 1997, 36, 2262-2270

A Hybrid Parametric/Stochastic Programming Approach for Mixed-Integer Linear 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 novel hybrid parametric/stochastic approach for the solution of stochastic mixedinteger problems is introduced. The approach, based on a two-stage stochastic programming framework recently proposed by the authors, effectively avoids the repetitive solution of the second-stage subproblems by solving instead multiparametric programs, resulting in a much more efficient evaluation of the expected value. The approach is detailed with two example problems in which computational issues are also discussed. Introduction In the area of process systems engineering, mixedinteger linear programming (MILP) models have been successfully applied 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 the last decade, several works have also been reported on the incorporation of uncertainty considerations into MILP models mainly through “multiperiod” or scenario-based formulations (Sahinidis et al., 1989; Shah and Pantelides, 1992; Subrahmanyam et al., 1994) as well as stochastic formulations (Reinhart and Rippin, 1986; Ierapetritou and Pistikopoulos, 1994; Liu and Sahinidis, 1996). Interestingly, most of these works deal explicitly with models in which the uncertainties appear in the right-hand side (e.g., production requirements, availability of raw materials, etc.) of the linear inequalities. The development of algorithms for the efficient solution of these types of problems is thus clearly of theoretical and practical significance. In this paper, a hybrid parametric/stochastic approach for the solution of stochastic mixed-integer problems is introduced. Based on the stochastic programming algorithms presented in Acevedo and Pistikopoulos (1996), the new approach relies on the solution of multiparametric problems for the evaluation of the second-stage subproblems. In this way, the computational requirements for the evaluation of the expected profit become independent of the number of integration points required by the integration method. The approach is applied to stochastic MILP with uncertainties in the right-hand side (RHS), for which an efficient solution algorithm for the corresponding multiparametric linear programming subproblems is available. The paper is organized as follows. First, the application of two-stage stochastic formulations in process synthesis and design is briefly reviewed, highlighting previous developments on which this work is based.

Then, the proposed hybrid parametric/stochastic approach is introduced, and general characteristics are discussed. Two algorithms for the incorporation of these parametric subproblems in a two-stage stochastic programming strategy are described for the case of MILP models with uncertain parameters in the RHS. The application of the hybrid approach is illustrated through some example problems, where computational issues are also discussed. Two-Stage Stochastic Programming ApproachsAn Overview Two-stage programming has been considered an effective approach to model process synthesis and design problems under uncertainty, naturally differentiating between design variables, which remain fixed once selected, and operating variables, which can be adjusted during operation to achieve feasibility (Grossmann et al., 1983). In a similar way, general process optimization problems under uncertainty can be addressed following a two-stage stochastic programming strategy, where, in the first stage, “here-and-now” decisions are made (e.g., selection of the plant structure and design parameters), whereas in the second stage the operation of the plant is optimized for specific values of the firststage variables and the uncertain parameters (wait-andsee decision). Mathematically, the two stages can be posed as follows:

Synthesis and Design (First) Stage h (y,d,zj,xj,θ)} - C(d) - cTy} R ) max {Eθ{P y,d

s.t.

h1(y,d) ) 0 g1(y,d) e 0 y ) {0, 1}m

Operating (Second) Stage P h (y,d,zj,xj,θ) ) max P(yj,d h ,z,x,θ) z

* 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)00708-7 CCC: $14.00

(1)

s.t.

h2(yj,d h ,z,x,θ) ) 0 h ,z,x,θ) e 0 g2(yj,d

(2)

∀θ ∈ J(θ) ∀θ ∈ J(θ)

where y is the vector of binary 0-1 variables defining © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997 2263

the flowsheet structure (e.g., existence of units), d is the vector of design variables (physical parameters as volumes, sizes, etc.), z and x are the vectors of control and state variables (operating conditions); θ represents the vector of uncertain parameters, described through a general probability function J(θ), and Eθ is the expectancy operator. The set of equalities h1 and h2 denote process equations (heat and mass balances, equilibria relations), the set of inequalities g1 and g2 correspond to design specifications and logical constraints for the first and second stages, respectively; P represents a scalar objective function, typically an economic performance index which must be maximized or minimized, such as a general profit function; C(d) represents the capital cost function of the plant, and c denotes a fixed charge cost vector. The overbar implies that the corresponding variables are held constant. Pai and Hughes (1987) proposed modeling and computational strategies for the solution of problems similar to (1) for the case of design under uncertainty (i.e., for fixed values of y). Reinhart and Rippin (1986) considered the design of flexible batch plants with uncertain demands through a similar two-stage programming formulation. Several integration methods for the approximation of the expectancy were proposed based on the linearity of the model and for the case of normal distribution functions for the description of uncertainty. Wellons and Reklaitis (1989) also presented a two-stage programming formulation for the design of multiproduct batch plants. In their work, the constraints were characterized as “hard” and “soft” (following the earlier work of Johns et al., 1978) to denote whether a constraint should be satisfied for all the possible values of the uncertain parameters or partial feasibility could be allowed. Straub and Grossmann (1990, 1993) presented a strategy for the evaluation of the stochastic flexibility of a design, i.e., the maximum ability of a design to accommodate uncertainty given by some probability distribution function, and its optimization. Their procedure, based on a generalized Benders decomposition (Geoffrion, 1972) type of approach, decomposes the problem by selecting the design variables as “complicating variables” and the second-stage subproblem in a series of smaller subproblems with which the feasible region of the design can be evaluated. To approximate the multiple integral representing the stochastic flexibility, the authors applied a numerical integration scheme in which the integration points are placed inside the depicted feasible region. On the basis of that work, Pistikopoulos and Ierapetritou (1995) proposed a two-stage formulation for process design problems under uncertainty where an optimal design is selected based on the maximization of an expected profit function. Their two-stage design strategy allows for the selection of an optimal design in the first stage, whereas in the second stage, the expected profit of the plant is approximated through a Gaussian quadrature method by optimizing its operation for different realizations (integration points) of θ. In this way, the trade-offs between process flexibility and economic optimality can be properly exploited, with which an optimal degree of stochastic flexibility can be obtained. Their proposed iterative solution procedure involves three major steps: 1. The feasible region of a fixed vector of structural and design variables (yj, d h ) is depicted through a series of “feasibility subproblems”, where the range of feasibility of the plant is evaluated for each uncertain parameter as follows:

max {θui - θli}

(3)

z,θ

s.t.

u h ,zu,xu,θ h j)1,...,i-1,θj)i,...,N )e0 fu(yj,d θ l fl(yj,d h ,zl,xl,θ h j)1,...,i-1,θj)i,...,N )e0 θ

θmin e θl e θ e θu e θmax where f represents all the constraints (equalities and inequalities) of the original model. 2. By placing the integration points (θqi) within the feasible region, the optimal profit of the plant is evaluated (“profit subproblems”) at each θqi, by solving

max Pq(yj,d h ,z,x,θq)

(4)

z

s.t.

f(yj,d h ,z,x,θq) e 0

q ) q1, ..., Q

A Gaussian quadrature formula is then used to approximate the expected profit of the given design, from which a valid lower bound on the final solution is defined. 3. Based on the information obtained through the solution of the subproblems in (3) and (4), the following master problem is generated and solved to obtain a new structure/design:

RUB ) max µ

(5)

y,d,µ

µ e Rk(y,d) -

s.t.

∑

ηkf(‚) -

Q){q1,...,qNθ} Nθ

{ ∑ i)2

∑

Q){q1,...,qi-1}

l,k l u,k u u λu,k i f (‚) + λi f (‚)} - λ1 f (‚) l λl,k 1 f (‚)

k ) 1, ..., K

h1(y,d) ) 0 g1(y,d) e 0 where k ) 1, ..., K represents the number of iterations. The solution of this master problem, RUB, yields an upper bound to the final solution, which is compared to the lower bound defined by the solution of the subproblems to evaluate whether convergence has been achieved or not. This work was later extended by Ierapetritou et al. (1996) to allow for the incorporation of different types of uncertainties and design objectives (in terms of flexibility requirements) in a unified formulation. Recently, Acevedo and Pistikopoulos (1996) studied the computational performance of alternative numerical and statistical integration schemes proposed for the evaluation of the expectancy term in (1) within a solution procedure as outlined above. The integration methods considered in their work were as follows: Scheme I: numerical integration evaluated within the entire space of the uncertain parameters. In this case, the integration points are placed “a priori” according to general Gaussian quadrature formulas, thus avoiding the first step of the original algorithm, i.e., the identification of the feasible region as in (3). The expected profit is then approximated by solving the profit subproblems in (4) at each of these points; if a subproblem is found infeasible, the point is discarded.

2264 Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997

Scheme II: numerical integration with a Gaussian quadrature formula within the simultaneously depicted feasible region, as in the original algorithm of Pistikopoulos and Ierapetritou (1995) outlined above. Scheme III: a (pseudo) Monte Carlo integration scheme, where the integration points are sampled from the probabilistic distribution functions of the uncertain parameters. Advantages and disadvantages of each integration scheme for the solution of different types of problems are discussed, as well as their implementation on different computational platforms. Despite these advances, the main difficulty for the application of such two-stage stochastic programming algorithms to large-scale problems remains the rapid explosion of the number of subproblems that must be solved in order to derive an accurate approximation of the expected value of the objective function considered. In an attempt to address this problem, a parametric programming approach, integrated within the two-stage stochastic programming algorithms, is proposed in the next section. This hybrid solution strategy avoids the repetitive optimization of subproblems at different values of the uncertain parameters, resulting in a procedure for which the solution times are practically independent of the number of integration points (samples) required. Hybrid Parametric/Stochastic ApproachsGeneral Concepts The two-stage programming algorithm described in the previous section for the solution of stochastic mixedinteger problems requires the solution of a potentially large number of optimization subproblems to evaluate the expected profit. These subproblems are formed by fixing the value of the uncertain parameters at specific points, provided by the integration scheme selected. The actual number of subproblems to be solved varies depending on the integration scheme and, most critically, the number of uncertain parameters in the model. Although several strategies have been suggested to reduce the computational requirements for the optimization of two-stage stochastic models (see for example Reinhart and Rippin, 1986; Straub and Grossmann, 1993; Acevedo and Pistikopoulos, 1996), the efficient evaluation of the expectancy still presents a major obstacle for the practical application of these formulations. On the other hand, the main impetus for solving an optimization problem parametrically is to avoid the unnecessary repetition of solutions of the problem for different values of the uncertain parameters. Instead, the parametric solution provides explicit functions representing the optimal solution and optimal objective value for all the values of the uncertain parameters for which an optimal solution exists in a given range. Since the structural and design variables are fixed in the “design stage” at each iteration of the stochastic algorithm in the previous section, the parametric solution of the “operating stage” subproblem, if available, represents the optimal operating policy for the selected plant for any value of the uncertain parameters in the range of values under study. Furthermore, the optimal value functions of this parametric solution could be integrated directly to obtain the expected profit of the plant or can be used to evaluate the values of the optimal profit and Lagrange multipliers necessary to formulate a new master problem in any of the integration schemes described in the previous section.

In the next section, the direct incorporation of the parametric solution into the stochastic programming algorithm is described for the specific case of MILP formulations with uncertainties in the right-hand side (RHS). Two example problems are then presented to illustrate the steps of the proposed approach. Hybrid Parametric/Stochastic Algorithms for MILP Problems with RHS Uncertainties Consider the case of mixed-integer linear models with uncertainties in the RHS. Then, it becomes clear from (2) that the second-stage problem results in a multiparametric linear program (mpLP), for which the algorithm of Gal and Nedoma (1972) (see also Gal, 1995) can be applied for its solution. The algorithm for mpLPs (described in Appendix I) is based on the evaluation of all possible optimal bases of the problem by incorporating the coefficients of the uncertain parameters into the Simplex tableau. The solution is given by a set of optimal value functions in terms of the vector of uncertainties, the corresponding optimal solution functions, and a set of constraints defining the space of optimality of each function. The optimal solution of the problem for specific values of the uncertain parameters can then be obtained by evaluating the optimal value function corresponding to the optimality (critical) space to which the set of values belongs. This then implies that the expected profit in the stochastic algorithm could be approximated simply by evaluating the parametric solution at the integration points/samples required, without actually re solving the optimization subproblems. In the following, we incorporate such a solution procedure for mpLPs into the two-stage programming algorithm presented in the first section, for the different integration methods proposed to approximate the expected profit, namely, Schemes I/III, which do not require the determination of the feasible region for a given design, and Scheme II, in which the feasible region is explicitly identified. Integration without Explicit Evaluation of the Feasible Region In this case, the second-stage involves profit optimization subproblems which can be reformulated as the following multiparametric LP:

max P(yj,d h ,z,x) z

s.t.

(6)

f(yj,d h ,z,x) e b + Fθ θmin e θ e θmax

where the functions P and f are linear, with f representing all the constraints (equalities and inequalities) of the original model. θ is the uncertain vector parameter and F is a constant conforming matrix. Problem (6) can then be solved parametrically through the algorithm given in Appendix I. If an integration scheme is selected where the integration points are placed within the entire uncertain parameter space, either through integration Scheme I involving Gaussian quadratures or Scheme III involving a Monte Carlo type of algorithm, the solution of the subproblems can then be obtained by evaluating the solution of (6) at the values of θ given by the corresponding numerical quadratures or the sampling pro-

Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997 2265

cedure. The resulting hybrid parametric/stochastic algorithm consists of the following steps: 1. Select an initial structure and design (yj, d h ). 2. Solve the parametric LP in (6) to obtain the optimal profit as a function of the vector of uncertain parameters, θ. 3. Define the fixed values of θ for the approximation of the expected profit either by specifying quadratures or by sampling. 4. Approximate the value of the expectancy by evaluating the profit function P h from the parametric solution obtained in step 2. 5. Form and solve a master problem as in (5) to obtain a new structure and design. 6. If the convergence criterion is not met, return to step 2. Note that, in this case, the Lagrange multipliers (dual values) of the active constraints, necessary to formulate the master problem, are given by the optimal tableau(s) of the parametric solution; these values do not depend on θ and therefore are constant at each parametric region of optimality (critical region). Furthermore, the parametric solution of the subproblem in (6) implicitly provides the feasible region of operation of the plant, given by the union of all critical regions of (6); this information, however, cannot be directly passed to the master problem in the present form of the algorithm. Finally, given that the execution time is mostly dependent on the solution of the parametric program, one could, in principle, select a very large number of points to evaluate the expectancy; this would allow the utilization of numerical integration techniques even if the number of uncertain parameters is large. Of course, this may not be necessary if a solution of the same accuracy can be obtained with a smaller number of points using efficient sampling techniques within a Monte Carlo type of integration. Integration with Explicit Evaluation of the Feasible Region If Scheme II is employed, i.e., a procedure involving the explicit evaluation of the feasible region, parametric problems can be constructed for the solution of the feasibility and profit subproblems as follows. The feasibility subproblems in (3) can be reformulated as the following multiparametric LP:

max {θui - θli} z,θi,θk

s.t.

(7)

u h ,zu,xu,θui ,θk)i+1,...,N ) e b + Fθj)1,...,i-1 fu(yj,d θ l fl(yj,d h ,zl,xl,θli,θk)i+1,...,N ) e b + Fθj)1,...,i-1 θ

0 e g + Gθj)1,...,i-1 θmin e θl e θ e θu e θmax where θj, j ) {1, ..., i - 1}, is the parametrization subvector, θi is the variable whose feasible range is being evaluated, and θk, k ) {i + 1, ..., Nθ}, are variables with bounds defined by the maximum and minimum possible values of the corresponding uncertain parameters. The new constraint on θj represents the bounds on the parameters from the previous i - 1 feasibility subproblems.

Note that, in contrast with the algorithm in the first Nθ i-1 section where a total of ∑i)1 Q deterministic feasibility subproblems must be solved, the parametric approach requires only the solution of problem (7) for i ) {1, ..., Nθ}; that is, Nθ parametric LPs are solved, each involving i - 1 uncertain parameters. Once the feasibility subproblems have been solved, problem (6) is solved to obtain the optimal profit of the plant as a function of θ. The quadratures for the evaluation of the expected profit can then be placed inside the feasible region since the information from the feasibility subproblems will be passed onto the new master problem. The resulting hybrid parametric/ stochastic algorithm with explicit evaluation of the feasible region then involves the following steps: 1. Select an initial structure and design (yj, d h ). 2. For i ) 1-Nθ: (a) Solve the feasibility mpLP subproblem in (7). Include θi in the subvector of parameters θj. (b) Evaluate the solution in step 2a at the fixed integration points of the parameters in θj. (c) Include the expressions of θui and θli as functions of θj as new bounding constraints (g + Gθj) for the next feasibility (multiparametric) subproblem. 3. Solve the multiparametric LP in (6) to obtain the optimal profit as a function of θ. 4. Approximate the value of the expectancy by evaluating P h from the parametric solution obtained in step 3 at the fixed values of θ (quadratures or samples). 5. Form and solve a master problem as in (5) to obtain a new structure and design for the process. 6. If the convergence criterion is not met, return to step 2. Note that in step 2 of the algorithm, i.e., the solution of the feasibility subproblems for the first uncertain parameters, it may be computationally more efficient to solve the deterministic problems for fixed values (quadratures) of the i - 1 parameters of (7). This, however, depends on the number of integration points selected and the ratio between the solution times of the parametric problem and the deterministic one. An Extension for Problems Involving Continuous Design Variables In recent computational studies for solving process synthesis problems under uncertainty via two-stage programming strategies (Acevedo and Pistikopoulos, 1996), it was observed that the (structural) 0-1 variables converge faster to their optimal values than the continuous design variables; i.e., a large percentage of the major iterations was consumed to adjust the design variables of only a relatively few structures. In such cases, if a parametric approach is employed for the solution of the subproblems, it seems advantageous to include the design variables as parameters varying between upper and lower bounds. For the case of the profit maximization subproblems, for instance, such a reformulation would lead to the following multiparametric problem:

max P(z,x) z

s.t.

f(yj,z,x) e Fθ - F′θd θmin e θ e θmax dmin e θd e dmax

(8)

2266 Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997

where θd represents the extra parameters introduced corresponding to the design variables and F′ the matrix of their coefficients in the original constraints. The new formulation has the advantage that its solution represents the optimal profit for a given structure for any value of the design variables (within the specified limits) as a function of the uncertain parameters. As a result, the total number of parametric problems that are required (in the case of integration without explicit evaluation of the feasible region) is equal to the number of different structures (integer solutions) that appear throughout the iterations of the stochastic algorithm, which is typically much less than the number of iterations. The penalty, on the other hand, is the increase in the computational complexity of the multiparametric subproblems due to the increase in the number of parameters. Benefits from this extended formulation will only be realized then if the ratio between the solution times of the two hybrid formulations is less than the number of iterations with repeated structures in the original stochastic programming algorithm. Note that this reformulation does not alter the way in which the two-stage stochastic algorithm progresses: the expected profit is still evaluated according to the integration scheme selected for the fixed values of the structural and design variables given by the master problem. The master problem, in turn, corresponds exactly to the formulation in (5). In the next section, the performance of the proposed hybrid parametric/stochastic algorithms is analyzed through the solution of two example problems. Numerical Examples Two examples from the open literature are revisited here to compare the stochastic algorithm to the hybrid formulations presented in this work. In both cases, a Monte Carlo integration scheme is considered for the approximation of the expected value. For the case of the purely stochastic algorithm, the implementation is based on GAMS/OSL (Brooke et al., 1988) for the solution of the subproblems and master problems. For the hybrid formulation, an experimental implementation of the mpLP algorithm in Appendix I was used for the subproblems and GAMS/OSL for the master problem. CPU times reported were obtained on a Sun SPARC Station 10-51. Example 1sSynthesis of a Utility Plant The first example, taken from Biegler et al. (1997), 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 1, 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; letdown 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. 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,b), is given in Table 1. The multiparametric MILP consists of 30 constraints, 8 integer variables, 24 continuous variables, and 4

Table 1. Mathematical Model for the Utility Plant Example 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′) mass balances HP ) HPMP1 + HPLP1 + HPMP2 + HPLP2 + HPLD MPProd ) MP + HPMP1 + HPMP2 + HPLD MPLP1 - MPLP2 - MPLD LPProd ) HPLP1 + HPLP2 + MPLP1 + MPLP2 + MPLD energy balances EHPT ) EHP1 + EHP2 EMPT ) EMP1 + EMP2 EHP1 ) 42.6HPMP1 + 121HPLP1 EHP2 ) 42.6HPMP2 + 121HPLP2 EMP1 ) 78.4MPLP1 EMP2 ) 78.4MPLP2 demands MPProd g 20 + 10θ1 LPProd g 80 + 10θ4 EHP1 + EMP1 g 7000 + 1000θ2 EHP2 + EMP2 g 4000 + 1000θ3 logical constraints 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 y ) (YT1, YT2, YT3, YT4, YT1′, YT2′)

cost functions

Table 2. Statistical Parameters for the Utility Plant Example demand

mean value

variance

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

25 tons/h 7500 kW 4500 kW 85 tons/h

1.67 167 167 1.67

uncertain parameters, defining the range of values that the electricity and steam demands may take. Normal distribution functions are assumed to describe these uncertainties with statistical parameters (mean and variance) as given in Table 2. The stochastic MILP was first solved following the two-stage approach, formulating the problem as the minimization of the expected cost. The optimal structure, found after 15 iterations, is given by the highpressure boiler and turbines T2 and T3, without expansions, with an expected cost of $1487302.88/yr, evaluated with 100 samples. Because of the restrictions in the maximum flow through the different equipment, many structures found throughout the iterations were infeasible for any value of the uncertain parameters. To avoid the optimization at each of the samples in these cases and therefore reduce the final solution time, a “feasibility” subproblem was implemented to check first whether a feasible region existed or not. (This does not represent a problem in the hybrid parametric/stochastic formulation since the infeasible structures are detected in the first phase of the multiparametric LP algorithm.) Following the proposed hybrid parametric/stochastic approach, the solution of this problem, for the first iteration, consists of the following steps: 1. An initial structure is selected including the highpressure boiler and turbines T1, with expansion T1′, and T4 (y ) {1, 0, 1, 0, 0, 1, 1, 0}). 2. The solution of the multiparametric LP, problem (6), resulting from fixing the integer variables with the

Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997 2267

Figure 1. Superstructure for the utility plant synthesis.

selected initial structure, is given by:

Cost1 ) 14.1 + 0.622θ1 + 1.043θ2 + 0.938θ3 CR1 )

{

3.52θ1 - 8.26θ2 - 8.26θ3 + 10θ4 e 3.87 θe1 ∀θ

Cost2 ) 13.73 + 0.96θ1 + 0.25θ2 + 0.145θ3 + 0.96θ4 CR2 ) -5.43θ1 + 12.75θ2 + 12.75θ3 - 15.43θ4 e -5.97 θe1 ∀θ

{

3. The uncertain parameters are sampled (100 samples) and the expected value is approximated to obtain an expected cost of $1490238.96/yr. 4. The solution of the master problem yields a new structure given by the vector y ) {1, 0, 1, 0, 0, 1, 0, 0} with which a new major iteration is started. The same optimal structure is found after 15 iterations. The parametric solution for this structure is given by

Cost1 ) 14.1 + 0.96θ1 + 1.369θ2 + 0.25θ3

{

-12.75θ2 + 10θ4 e 9.286 θe1 ∀θ

CR1 )

Cost2 ) 11.61 + 0.96θ1 + 0.145θ2 + 0.25θ3 + 0.96θ4

{

-10θ1 + 23.47θ3 - 10θ4 e -6.1 CR2 ) 12.75θ2 - 10θ4 e -9.286 θe1 ∀θ Cost3 ) 11.03 + 0.145θ2 + 2.5θ3 CR3 )

{

10θ1 - 23.47θ3 + 10θ4 e 6.1 θe1 ∀θ

Using the hybrid parametric/stochastic approach, the solution of the multiparametric profit subproblems required a total of approximately 1.6 CPU s. In this implementation, the solution of the parametric problems

took 6-7 times longer than their corresponding base cases (first phase of the multiparametric LP algorithm) solved as deterministic LPs. In comparison, the solution of the profit subproblems (including the feasibility subproblems) using the purely stochastic programming algorithm took approximately 6.15 CPU s. Using 20 samples to approximate the expected cost, 14.69 CPU s and using 50 samples, and 28.61 CPU s using 100 samples. The solution of the master problems in all cases took approximately 4.4 CPU s in total. In addition to savings in the solution time, another advantage of the approach is that a full sensitivity analysis is obtained for the different structures/designs examined, from which robust solutions can be easily identified, i.e., structures which are less sensitive to variations. Similarly, the most important uncertainties (and their corresponding ranges of values) can also be depicted. Example 2sA Capacity Planning Problem The capacity planning example presented in Acevedo and Pistikopoulos (1996) is revisited here to examine the performance of the two hybrid formulations proposed in the previous sections. The superstructure of the problem, shown in Figure 2, comprises 11 processes for the manufacturing of 5 products from 5 possible raw materials. The mathematical model, given in Table 3 with constants as given in Table 4, involves 11 integer variables, 57 continuous variables, 68 constraints, and 10 uncertain parameters given by the demands of the products and the availability of the raw materials. Linear yield relations have been assumed, resulting in a stochastic MILP model. The starting point selected was given by processes 2, 4, 5, 7, 8, and 10, i.e., y ) (0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0), with capacities Q ) {0, 2, 0, 2, 2, 0, 3, 1, 0, 2, 0}. Five thousand samples and 21 iterations of the stochastic programming algorithm were needed to converge to the optimal solution given by y ) (1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1) with capacities (Q ) {1.5, 0, 0, 1.9, 1.6, 0, 3.3, 1.5, 1.6, 0, 1.9}), for a total of approximately 6000 CPU s.

2268 Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997

Figure 2. Superstructure for the capacity planning example. Table 3. Mathematical Model of the Capacity Planning Example 5 (PricePiPi) - ∑5j)1(CostRMjRMj) - ∑11 P ) Eθ[∑i)1 k)1(OCkISk)] (DCkQk) - (FCkyk) OSk ) PckISk ∀k Pi e Di(θ) RMj e MaxRMj(θ) ISk - MIkQk e 0 ∀k Qk - MaxQkyk e 0 ∀k price of product i cost of raw material j uncertain demand of product i uncertain maximum availability of raw material j capacity of process k input stream to process k output stream to process k yield constant for process k flow to volume relation for process k maximum capacity of process k

objective function yield relations desired production availability of logical constraints Pricei Costj Di MaxRMj Qk ISk OSk Pck MIk MaxQk Table 4. Constants of the ModelsExample 2 Pc Mi Pc Mi

Proc. I

Proc. II

Proc. III

Proc. IV

Proc. V

Proc. VI

1.5 18

1.25 20

1 15

0.8 20

2 20

1.25 21

Proc. VII

Proc. VIII

Proc. IX

Proc. X

Proc. XI

0.75 15

1.5 15

0.9 25

0.75 15

0.85 20

Following the hybrid approach, the problem was first solved following the algorithm proposed in the previous section, taking 5000 samples for the approximation of the expected profit, evaluated through the solution of problem (6) for fixed values of the structural and design variables, at each major iteration. The same optimal solution, consisting of processes 1, 4, 5, 7, 8, 9, and 11, was obtained after 21 iterations. The base case for these problems took approximately 0.05 CPU s; three or four optimal bases were found for the different designs in the parametric solution, consuming from 1.3 to 1.7 CPU s for its completion. The total solution time for the parametric profit maximization subproblems was approximately 31 CPU s, representing a 2 orders of

magnitude decrease in the computational time required compared to the original algorithm! The problem was then solved using the extended formulation proposed in the previous section, taking again 5000 samples for the approximation of the expected profit, evaluated this time through the solution of problem (8) for fixed values of the structural variables. In order to avoid numerical problems, the new parametric subproblems were “preprocessed” every time a new structure was encountered, defining only as parameters the capacities of existing processes. On average then, six or seven new parameters were introduced; the range of these parameters was taken from 0.5 to 4, which covered every design selected by the algorithm. The optimal solution was again obtained after 21 iterations with the same optimal structure and capacities. The base case for this formulation took 0.07 CPU s, but this time five to seven optimal bases were found for the different structures, taking from 1.8 to 2.5 CPU s for the solution of each parametric subproblem, an increase of approximately 70% on average with respect

Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997 2269

to the first formulation. The advantage, however, is that the solution of only eight of these problems was required for the evaluation of the expected profit of the different structures/designs from every iteration, for a total solution time of approximately 19 CPU s (a reduction of a further 60%). As in the first example, another advantage of the hybrid approach is the sensitivity analysis information simultaneously obtained. In the case of the extended formulation, this analysis includes a complete profile of the effect of changes in the design on the optimal operation and profit. General Remarks From the computational results presented, a number of issues are worth highlighting. First, note that the implementation of the multiparametric LP used in these studies is based on a rather “rough” implementation of the Simplex algorithm; more importantly, however, it does not include an efficient way to identify the neighbors of an optimal basis (see Gal, 1995). Instead it relies on the solution of successive optimization (auxiliary) problems. This not only increases the solution time but also results in the examination of some bases which are not part of the final solution and are later discarded, e.g., spaces given by a single hyperplane, generated by the exchange of slack variables (the time spent on these calculations was not included in the reported CPU time). Despite this, the advantages of the new hybrid approach, compared to previous formulations, are apparent, and some clear trends have been observed from the solution of these examples. First, although further computational experimentation is clearly needed, the solution times in the examples presented have been significantly reduced (in some cases by 1 or 2 orders of magnitude). A more advanced implementation of the algorithms (currently under development) is required for the evaluation of their performance to large-scale problems. Another important trend is the modest increase of the solution time as a function of the number of uncertain parameters, which can be observed in the second example where the two formulations were compared. In this case, while the incorporation of the capacities as uncertain parameters represented an increase of more than 60% in the number of uncertain parameters (from the initial 10 uncertainties to 16 or 17), the increase in the average computational effort required for the solution of each parametric subproblem was less than 70%. As a result, a reduction of the overall solution time of more than 60% was obtained for this example when using the extended formulation. Finally, it is important to note that the proposed hybrid approach can be easily extended to allow uncertainties in the cost coefficients by including a solution algorithm for RIM (cost coefficients and right-hand-side) parametrizations of linear problems, such as the one given by Gal (1995), and ultimately to mixed-integer nonlinear formulations. The latter, however, would require the solution of parametric nonlinear programs for the evaluation of the second-stage subproblems, which still represents a very difficult problem to solve (Fiacco, 1983). Concluding Remarks A new hybrid parametric/stochastic approach for the solution of stochastic mixed-integer optimization problems was presented. The approach follows the steps of

the two-stage stochastic programming algorithm presented in earlier works, coupled with a multiparametric solution for the second-stage subproblems with which the evaluation of the expected profit can be performed much more efficiently. Alternative algorithms have been presented to incorporate the multiparametric problems into a stochastic optimization framework with different integration strategies, for the case of MILP models with RHS uncertainties. An extended formulation of the subproblems was also introduced for problems involving continuous design variables. Two example problems were presented to compare the hybrid approach with the original stochastic algorithm and highlight the advantages of the proposed formulations. 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 Variables y ) 0-1 integer (structural) variables d ) continuous design variables z ) continuous control variables x ) continuous state variables yj ) fixed values of the 0-1 integer variables d h ) fixed values of the design variables P(y,d,z,x,θ) ) profit function P(yj,d h ,z,x,θ) ) profit function for fixed structure and design Nθ ) number of uncertain parameters Q ) number of integration points per parameter Greek Symbols θ ) vector of uncertain parameters θd ) vector of design variables defined as uncertain parameters

Appendix I. Algorithm for mpLP Problems (Gal, 1995) The parametric/stochastic algorithm proposed in this work relies on the solution of multiparametric linear subproblems (mpLP) at each major iteration of the stochastic programming algorithm. These subproblems are of the following general form:

z(θ) ) min cTx x

s.t.

(9)

Ax e b + Fθ Gθ e g

x e 0;

x ∈ Rn

θ ∈ Rs where A is an m × n matrix, F is an m × s matrix, and G is an r × s matrix; c, x, b, and g are conforming vectors. The elements of A, F, G, c, b, and g are constant, while vector θ is a vector of uncertain parameters varying in Ξ (Ξ ) {θ|Gθ e g; θ ∈ Rs}). All vectors are column vectors. An algorithm for solving parametric LP problems is given by Gal and Nedoma (1972). The basic steps of the solution procedure are given here; further details

2270 Ind. Eng. Chem. Res., Vol. 36, No. 6, 1997

on the theory and the algorithm itself can also be found in the excellent book by Gal (1995). 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∈Ξ

where K is the set of all possible vector parameters for which problem (9) has a finite optimal solution. 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 (9). If no feasible point can be found, then stop; the 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 bases 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 basis 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 (9). Literature Cited Acevedo, J.; Pistikopoulos, E. N. Computational Studies of Stochastic Optimization Algorithms for Process Synthesis under Uncertainty. Comput. Chem. Eng. 1996, S20, 1-778. Andrecovich, M. J.; Westerberg, A. W. A Simple Synthesis Method Based on Utility Bounding for Heat-Integrated Distillation Sequences. AIChE J. 1985, 31, 363-375. Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods for Chemical Process Design; Prentice Hall: New York, 1997. Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User’s Guide; Scientific Press: Redwood City, CA, 1988. Fiacco, A. V. Introduction to Sensitivity and Stability Analysis in Nonlinear Programming; Academic Press: New York, 1983. 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-422. Geoffrion, A. M. Generalized Benders Decomposition. J. Optim. Theory Appl. 1972, 10, 237-260.

Grossmann, I. E.; Haleman, K. P.; Swaney, R. E. Optimization Strategies for Flexible Chemical Process. Comput. Chem. Eng. 1983, 7, 439-462. Ierapetritou, M. G.; Pistikopoulos, E. N. Novel Optimization Approach of Stochastic Planning Models. Ind. Eng. Chem. Res. 1994, 33, 1930-1942. Ierapetritou, M. G.; Acevedo, J.; Pistikopoulos, E. N. An Optimization Approach for Process Engineering Problems Under Uncertainty. Comput. Chem. Eng. 1996, 20, 703-709. Johns, W. R.; Marketos, G.; Rippin, D. W. T. The Optimal Design of Chemical Plant to Meet Time-Varying Demands in the Presence of Technical and Commercial Uncertainty. Trans. Inst. Chem. Eng. 1978, 56, 249-257. 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-227. Liu, M. L.; Sahinidis, N. V. Optimization in Process Planning under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 4154-4165. Pai, C. C. D.; Hughes, R. R. Strategies for formulating and solving Two-Stage problems for process design under uncertainty. Comput. Chem. Eng. 1987, 11, 695-706. Papoulias, S.; Grossmann, I. E. A Structural Optimization Approach in Process Synthesis. Part I: Utility Systems. Comput. Chem. Eng. 1983a, 7, 695-734. Papoulias, S.; Grossmann, I. E. A Structural Optimization Approach in Process Synthesis. Part II: Heat Recovery Networks. Comput. Chem. Eng. 1983b, 7, 707-721. Pistikopoulos, E. N.; Ierapetritou, M. G. A Novel Approach for Optimal Process Design under Uncertainty. Comput. Chem. Eng. 1995, 19, 1089-1110. Raman, R.; Grossmann, I. E. Symbolic Integration of Logic in Mixed-Integer Linear Programming. Comput. Chem. Eng. 1993, 17, 909-927. Reinhart, H. J.; Rippin, D. W. T. The Design of Flexible Batch Chemical Plants. Presented at the AIChE Annual Meeting, Miami, FL, 1986; Paper 50e. 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. 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, 967985. Straub, D. A.; Grossmann, I. E. Design Optimization of Stochastic Flexibility. Comput. Chem. Eng. 1993, 17, 339-354. Subrahmanyam, S.; Pekny, J. F.; Reklaitis, G. V. Design of Batch Chemical Plants under Market Uncertainty. Ind. Eng. Chem. Res. 1994, 33, 2688-2701. 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-1326. Wellons, H. S.; Reklaitis, G. V. The Design of Multiproduct Batch Plants under Uncertainty with Staged Expansion. Comput. Chem. Eng. 1989, 13, 115-126.

Received for review November 5, 1996 Revised manuscript received February 21, 1997 Accepted February 24, 1997X IE960708F X Abstract published in Advance ACS Abstracts, April 15, 1997.