Ind. Eng. Chem. Res. 2002, 41, 67-77
67
A Hybrid Parametric/Stochastic Programming Approach for Mixed-Integer Nonlinear Problems under Uncertainty Tanya S. Hene´ , Vivek Dua, and Efstratios N. Pistikopoulos* Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College, London SW7 2BY, U.K.
This paper presents a hybrid parametric stochastic programming approach for mixed-integer convex nonlinear optimization problems under uncertainty. This approach is based upon an iterative two-stage stochastic optimization framework, where in the first stage design and integer variables are fixed and an expected profit is calculated by solving nonlinear programs (NLPs) at the integration points, in the space of uncertain parameters. Then, in the second stage, a master problem is formulated, based upon the dual information obtained by solving NLPs in the first stage, and solved to identify a new vector of design and integer variables. The solution procedure terminates when the solution of the first- and second-stage problems is within a certain tolerance. The basic idea of the hybrid approach, proposed in this work, is that in the first stage parametric programming techniques are used to obtain profit as a function of uncertain parameters. This reduces the computation of the expected profit to a function evaluation of the profit function at the integration points. Therefore, the solution of the NLPs at the integration points is avoided. The use of efficient integration techniques, which require a smaller number of integration points, in light of the hybrid approach is also discussed. P ) max f(x,y,d,θ) x,y,d
1. Introduction In many engineering problems, one usually encounters two major types of decisions. The first type of decision is usually associated with the selection or not of various components, such as processes, equipments, etc., whereas the second type of decision corresponds to variables which are continuous in nature, such as flow rates and temperatures. In an optimization framework, such problems can be formulated as follows:1,2
P ) max f(x,y,d) x,y,d s.t.
(1)
g(x,y,d) e 0 x ∈ Rn y ∈ {0, 1}m d ∈ Rl
where P is the profit, x is a vector of continuous variables, y is a vector of 0-1 binary variables, d is a vector of design variables, f is a scalar objective function which is continuously differentiable and a concave function of x, and g is a vector of constraints which are continuously differentiable and convex functions of x; f and g are linear and separable in d and y. For the case when some parameters in the model are uncertain, (1) can be reformulated as follows:3 * To whom correspondence should be addressed. E-mail:
[email protected]. Tel: +44 20 7594 6620. Fax: +44 20 7594 6606.
s.t.
(2)
g(x,y,d,θ) e 0 x ∈ Rn y ∈ {0, 1}m d ∈ Rl θ ∈ Rs
where θ is a vector of uncertain parameters which is present linearly and separately. θ is usually associated with (i) commercial parameters such as uncertainty in demands and selling prices of products and supply and cost prices of raw materials and (ii) technical parameters such as heat-transfer coefficients. The solution techniques4 for (2) can be broadly classified based upon (i) multiperiod,5-8 (ii) stochastic programming,9-13 and (iii) parametric programming principles.14,15 These three approaches for the solution of optimization problems under uncertainty are based upon different descriptions of the uncertainty and result in different types of solutions. The multiperiod approach discretizes a given time horizon into a finite number of time periods. Each time period is characterized based upon the values of the uncertain parameters forecasted in that period. The objective is to calculate the optimal net present value. In the stochastic programming approach, the uncertainty is described by a probability distribution function and the objective is to calculate the optimal “expected profit”. The aim in the parametric programming approach is to obtain the optimal objective function and the corresponding optimization variables as a function of the uncertain parameters. -
10.1021/ie0100582 CCC: $22.00 © 2002 American Chemical Society Published on Web 01/02/2002
68
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002
mulas were presented;18,19 one such, fifth degree, formula gives the expected profit as follows:
Pkexp
1
)
π
s/2
Np
BiPkmax(θi) ∑ i)1
(5)
where θi is given as follows:
θi ) µθ + x2σθvi
Figure 1. Stochastic approach for nonlinear systems.
In this work, we will focus on the stochastic programming approach and “use” some features of the parametric programming approach to improve the efficiency of the stochastic programming approach. This work16 is an extension of the work for linear systems17 to problems involving convex nonlinearities. The rest of the paper is organized as follows. In the next section, we briefly review the stochastic programming approach and highlight the key features of this approach, and in section 3, we give an outline of the parametric programming approach and motivate how this approach can be used to improve the efficiency of the stochastic programming approach. Section 4 presents the theoretical and algorithmic framework for the hybrid approach, while a comparison of the computational requirements for the stochastic and hybrid approaches is illustrated with examples in section 5 and some concluding remarks are given in section 6. 2. Stochastic Programming
where µθ and σθ are the means and standard deviation for θ and vi and Bi are given in Table 1. The number of integration points, Np, required for this formula is 2s + 2s, and this formula can be used for s g 3. Note that, for example, when four uncertain parameters are present in the model, a fifth-degree Gaussian quadrature formula will require 625 integration points and a specialized cubature formula for a similar accuracy will require only 24 points. Note that we assume that (3) is feasible for every value of θ that lies in the given ranges. This feasibility, in general, can be ensured by introducing a penalty in the objective function for the infeasibilities.20 The expected profit, Pkexp, thus calculated represents a lower bound on the solution of (2). The dual information obtained from the solution of the NLPs in the first stage is used to formulate the second-stage problem as follows:
PUB ) max µ y,d,µ s.t.
µ e Pkexp(y,d) -
(7) λhqg(y,d,xk,q,θq) ∑ q∈Q
g(y,d,xk,q,θq) e 0, ∀ q ∈ Q, ∀ k ∈ K y ∈ {0, 1}m d ∈ Rl
Two-stage programming is a widely used technique for solving stochastic optimization problems described by problem (2). The basic idea of the approach, which relies on a GBD framework, is as follows.9 In the first stage, d and y in (2) are fixed and the following
Pkmax(θq) ) max f(yk,dk,x,θq) s.t.
(6)
(3)
g(yk,dk,x,θq) e 0 x ∈ Rn
nonlinear program (NLP) is solved with ∀ q ∈ Q, where k is the number of iterations, Q is the set of integration k points, Pmax (θq) is the optimal profit at an integration point, q (see Figure 1, where asterisks denote the points, θq, where NLPs are solved for dk and yk and θL and θU are the lower and upper bounds on θ). The expected profit is then calculated as follows:
where λhq is the vector of the corrected Lagrange multipliers.9 PUB, obtained from the solution of (7), is an upper bound to the profit, the solution of (2). The design and integer variables obtained from the solution of the second-stage problem are fixed, and the first-stage problem is solved again for these fixed new vectors of d and y. The algorithm proceeds in this way, and the bounds are updated at each iteration. The algorithm terminates when the updated lower and upper bounds are within a certain tolerance. Because the solution of the first-stage problem requires the solution of a large number of NLPs at the integration points, it is a computationally expensive step in the solution procedure. In the next section, we present an alternative approach based on parametric programming. This then results in explicit expressions of profit functions in terms of the uncertain parameters without exhaustively enumerating the complete space of uncertain parameters. This then results in a more efficient calculation of the expected profit. 3. Parametric Programming
Pkexp ) Eθ{Pkmax(θq)}
(4) Pkexp
is the where Eθ is the expectancy operator and expected profit at iteration k. Gaussian quadrature integration formulas are commonly used for the placement of the integration points in (3). Recently, specialized cubature integration for-
Parametric programming is a technique for obtaining the optimal solution as a function of uncertain parameters. While the theory and algorithms of parametric programming have been developed for various classes of problems,14,21-26 here we describe an approach for convex nonlinear models.15,27 First, a multiparametric nonlinear program (mp-NLP) is formulated by fixing d
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002 69 Table 1. Cubature Formula Placement and Weights points (v)
weights
(r, 0, ..., 0) ((ξ, (ξ, ..., (ξ)
B0 B1
Parameter Values r2 ) (s + 2)/4 ξ2 ) (s + 2)/[2(s - 2)] B0 ) 4/[(s + 2)2]V] B1 ) [(s - 2)2]/[2s(s + 2)2]V V ) πs/2
and y in (2) as follows:
P(θ) ) max f(x,yk,dk,θ) x s.t.
(8)
g(x,yk,dk,θ) e 0 x ∈ Rn θ ∈ Rs
where dk and yk are the fixed design and integer variables, respectively. Note that P(θ) is a concave and continuous function of θ.21 The basic idea of the approach is then to create outer approximations of the nonlinear terms in (8) as follows:
P ˆ (θ) ) max f(x*,yk,dk,θ) + ∇xf(x*,yk,dk,θ*)(x - x*) x (9) s.t.
g(x*,yk,dk,θ) + ∇xg(x*,yk,dk,θ*)(x - x*) e 0 x ∈ Rn θ ∈ Rs
where x* is the solution of (8) for dk and yk at θ ) θ*. The formulation in (9) is a multiparametric linear program (mp-LP). The solution of such an mp-LP is given by linear objective function profiles, P ˆ (θ), constant Lagrange multiplier profiles, λˆ (θ), and the regions in the space of θ where P ˆ (θ) remains optimal and λˆ (θ) remains constant.22 These regions are known as critical regions (CR). Because f is concave and g is convex, P ˆ (θ) provides an upper bound to, P(θ), the solution of (8); see Figure ˆ (θ) in the first step. In the next 2a, where P ˆ 1(θ) denotes P step, P ˆ 1(θ) is compared to P(θ) to find the maximum ˆ 1(θ). Because P(θ) difference, ∆Pmax, between P(θ) and P is concave, ∆Pmax will lie at one of end points of the critical regions of θ, the point in the space of θ where ∆Pmax lies is denoted by θmax. See Figure 2a again where ∆Pmax lies at θL; note that the triangles in Figure 2 indicate the points where NLPs are solved. If ∆Pmax is within a certain tolerance, , the solution of (8) is given by P ˆ 1(θ); otherwise, another mp-LP of the form (9) is formulated at the point, in the space of θ, where ∆Pmax lies and solved to obtain P ˆ 2(θ); see Figure 2b. P ˆ 1(θ) and P ˆ 2(θ), obtained at different points in the space of θ, θ*, and θL in Figure 2, are then compared by using the comparison procedure described in Acevedo and Pisˆ (θ), are tikopoulos (1999),26 and tighter linear profiles, P retained. These steps are repeated until all P ˆ (θ) are within to P(θ). These parametric profiles, P ˆ (θ), provide the profit, for fixed d and y, for each value of θ that lies in the given ranges. The expected profit, Pkexp, for the first-stage problem is therefore obtained by simply
substituting the values of θ at the integration points in P ˆ (θ) as shown in Figure 3; note that the asterisks denote the integration points. One must note that while the computational requirements for calculation of Pkexp by using the traditional stochastic programming approach increase with the number of integration points, when using the parametric programming, in the first stage, the number of integration points does not affect the computational requirements. This is because the parametric programming approach does not rely on the solution of the NLPs at the integration points and instead requires simple function evaluations of the parametric profiles, P ˆ (θ), at the integration points. These ideas are summarized in the next section. 4. Hybrid Parametric/Stochastic Programming Let us first put the things we have discussed so far into perspective. Two-stage stochastic optimization is widely used to calculate the expected profit for problems of the form (2) which involve uncertain parameters.9 The basic idea of the approach is to decompose the problem into primal and master subproblems which converge to the optimal solution as the iterations proceed. This decomposition is similar to the generalized Benders decomposition (GBD).2 First, the complicating variables are identified. The complicating variables, for example, correspond to design, d, and structural, y, variables in process synthesis and design problems. In the first iteration, the complicating variables are fixed at certain values and the resulting NLPs, (3), are formulated and solved for θ fixed at the integration points. These solutions are then used to calculate the expected profit by using formulation (4). The expected profit for fixed complicating variables provides the solution of the primal subproblem. The focus of this work is on efficient solution of the primal subproblems and is discussed in the next paragraph. The solution of the primal subproblem provides a lower bound on the solution of the original problem (2). Based upon the dual information obtained by solving the primal subproblem, the master subproblem of the form (7) is formulated and solved to obtain a new value for the complicating variables which is returned back to the primal subproblem. The iterations proceed in this way until the solutions of the primal and master subproblems are within a certain tolerance. The primal subproblem consumes most of the computational time. This is due to the large number of the optimization subproblems which have to be solved. This number depends on the number of the Gaussian quadrature integration points and increases exponentially with the number of uncertain parameters.9 Three approaches have been proposed to improve the efficiency of the solution of the primal subproblem: (1) parallel computing;11 (2) cubature integration techniques18,19 [formulation (5)]; (3) hybrid parametric/stochastic programming.17 In this work the focus is on extending the third approach which is applicable for linear systems to problems which also involve convex nonlinearities. The computational requirements of the extended third approach are also compared to the second approach in section 5. The extended third approach relies on obtaining the optimal objective function value as a function of the uncertain parameters by using the mp-NLP algorithm.27 The key point is that (multi)parametric programming is a general tool and provides the optimal
70
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002
Figure 2. Parametric approach for nonlinear systems.
Figure 3. Hybrid approach for nonlinear systems.
solution as a function of the uncertain parameters without exhaustively enumerating the entire space of the uncertain parameters. The expected profit calculation, therefore, reduces to a simple function evaluation of the optimal function profile. For the case of mp-NLP, the optimal function is obtained as a set of linear parametric profiles which approximate the solution within a certain tolerance. Therefore, the computational requirements for the solution of the mp-NLPs increase with a decrease in the tolerance. This affects the overall computational burden of the complete algorithm. This is also pointed out in section 5. Based upon the developments reported in sections 2 and 3 and the description above the hybrid parametric/ stochastic algorithm for mixed-integer nonlinear problems under uncertainty can be summarized as follows: 1. Initialize the lower and upper bounds on profit, PLB ) -∞ and PUB ) +∞, and the tolerances, 1 and 2. 2. At iteration k, fix the integer and design variables yk and dk. 3. Let the current upper bound for the parametric solutions be CUB(θ) ) +∞. For fixed yk and dk, formulate and solve an mp-NLP of the form given in (8) as follows: (a) Fix θ ) θ* in (8) and solve the resulting NLP. Let the solution be given by x*. (b) Create outer approximations of the nonlinear terms in (8) and formulate and solve an mp-LP of the form given in (9). Let the solution of this problem be given by P ˆ (θ) and λˆ (θ) and the corresponding CRs. Let ˆ (θ)}. CUB(θ) ) min {CUB(θ), P (c) Solve (8) for θ fixed at the vertexes of the CRs. Let the solution vector of this problem be given by P(θv), where θv is the vector of the vertexes of the CRs.
Figure 4. Design under uncertainty algorithm for nonlinear systems using the hybrid approach.
(d) Let ∆Pmax ) max {P(θv) - CUB(θv)}. If ∆Pmax e 1, θv go to the next step; else, go to (3a), with θ* as the vertex where ∆Pmax > 1. 4. Evaluate the expected profit as given in (4) and (5) k (θq) is given by CUB(θq). where Pmax 5. Update the lower bound, PLB ) max {PLB, Pkexp}. 6. Formulate and solve a master problem of the form (7) to obtain new d and y. Let the solution of this problem be given by dk, yk, and µk. Update PUB ) µk. 7. If PUB - PLB > 2, return to (2); else, stop with PLB as the final solution. It is assumed in step 3 that (8) is feasible for any value realization of θ. For the case when this is not true, minimization of infeasibility subproblems are formulated and solved and the corresponding Lagrangian infeasibility cuts are included in the master subproblem (7). The algorithm is illustrated in Figure 4.
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002 71 Table 2. Subproblem Requirement Comparison pure stochastic MILP LP mpLP
quadrature
cubature
hybrid
7 3750 0
7 144 0
7 0 6
In this work we have assumed f to be concave, g to be convex, and both f and g to be continuously differentiable functions of x. When this assumption is not valid, a new algorithm, for example, by using some of the concepts presented in Ierapetritou and Pistikopoulos (1996)20 and Dua et al. (1999)28 could be developed. In the next section, we illustrate the concepts of the hybrid parametric/stochastic programming with examples. 5. Examples This section illustrates the steps of the proposed algorithm with examples, and a comparison of the computational requirements is also presented. The first example is a linear utility synthesis problem where electricity and steam demands are uncertain. Details of this example are given elsewhere, and here only the computational requirements are reported. The remaining three examples are nonlinear synthesis and planning examples and are used to explain the solution steps in detail. 5.1. Example 1. Consider a utility plant problem14 which has been modified to include penalties for unmet demands.16 It concerns the design of a plant for given uncertain demands of electricity and steam at medium and low pressures. The demands are all assumed to be normally distributed. The expected cost at each iteration of the two-stage stochastic approach is calculated by taking 3.09σ of the distribution (99%) into account. The model involves continuous and integer variables which are present linearly. The computational requirements for this example are as follows. In the hybrid approach, a single mpLP was required for each iteration of the GBD algorithm, each of which required approximately 0.2 CPU s per primal problem. The pure stochastic approach required approximately 0.03 CPU s per integration point, a total of approximately 0.7 CPU s per primal problem. Table 2 shows the number of subproblems required in each method. In this example, the hybrid approach is therefore computationally more efficient than the pure stochastic method, despite the small number of integration points required in the numerical integration using cubatures. 5.2. Example 2. Consider a problem where two products are produced in a single process for which a structure needs to be chosen. Two reactors of different capacities are available, but at most one will be purchased. The demand for each product is uncertain, and whenever the demand is not met, a penalty is added to the cost. This penalty is in the form of an opportunity cost, the price of the product, and is proportional to the amount of demand that goes unmet. Nonlinearities appear only in the objective function. The price of product 1 is £10/unit, and for product 2 it is £12/unit. The purchase price for each reactor is £50. Operating costs are £0.2/unit2 and £0.3/unit2. The capacities of reactor types 1 and 2 are 7 and 9 units, respectively.
Figure 5. Solution after the first linearization.
The mathematical formulation is as follows:
max Profit ) 10M1 + 12M2 - 10(θ1 - M1) - 12(θ2 M2 ) - 0.2M12 - 0.3M22 - 50(y1 + y2) s.t. M1 + M2 e 7y1 + 9y2 M1 e θ1 M2 e θ2 y1 + y2 e 1 3.98 e θ1 e 6.02 3.98 e θ2 e 6.02 M1, M2 g 0 y1, y2 ∈ {1, 0} where M1 and M2 are the amounts of products 1 and 2 produced, θ1 and θ2 are the amounts of products 1 and 2 demanded, and y1 and y2 are the binary variables indicating the existence of processes 1 and 2. The demands of products 1 and 2 are both assumed to be normally distributed with a mean value of 5 and a standard deviation of 0.33. The upper and lower bounds are therefore set to include 99% (3.09σ) of this range. The solution steps are as follows. 1. Initialization: The initial design is selected by solving the optimization problem as an MINLP29 in GAMS30 and defining the uncertain demands as variables between their upper and lower bounds. The optimal profit is found at the following point: y1 ) 0, y2 ) 1, M1 ) 3.98, M2 ) 5.02, θ1 ) 3.98, θ2 ) 5.02. A tolerance of 0.3% is initialized to be the accuracy requirement for the nonlinear function for the linear parametric solution. 2. Linearization: For the selected design, Y ) (0, 1), at the optimal solution, θ* ) (3.98, 5.02), the nonlinear terms are linearized, and the following mp-LP is formulated:
max Profit ) 18.408M1 + 20.988M2 - 10θ1 12θ2 - 39.2718 (10)
72
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002
Table 3. Profit Comparison vertex
ProfitNLP
ProfitmpLP
% diff
1 2 3 4 5 6
29.64 28.17 29.39 8.99 39.31 38.17
29.96 28.71 29.89 9.49 39.31 28.71
1.09 1.91 1.71 5.58 0.00 1.41
s.t. M 1 + M2 e 9 M1 e θ1 M2 e θ2 Figure 6. Solution after the second linearization.
3.98 e θ1 e 6.02 3.98 e θ2 e 6.02 M1, M2 g 0 3. The Parametric Solution: The parametric solution of problem 10 is obtained by using the mpLP code developed by Acevedo (1996).31 The solution is given by the following two critical regions (Figure 5). CRA1 : θ1 + θ2 e 9, θ1 g 3.98, θ2 g 3.98, M1 ) θ1, M2 ) θ2, where ProfitA1 ) 8.408θ1 + 8.988θ2 - 39.2718. CRA2 : θ1 + θ2 g 9, 3.98 e θ1 e 6.02, 3.98 e θ2 e 6.02, M1 ) 9 - θ2, M2 ) θ2, where ProfitA2 ) -10θ1 - 9.42θ2 + 126.4002. 4. Finding the Point of the Largest Deviation: The next step is to compare the linear parametric profile to the nonlinear optimal solution. This is necessary to find out how far removed from the actual solution the profile is at this iteration. Because the problem is convex (maximization of a concave objective function), the point of largest deviation will be at one of the vertexes of the CR. There are six vertexes in this solution, and at each vertex a deterministic NLP is solved in GAMS. Table 3 shows the comparison of the optimal NLP solution with the overestimators from the parametric solution. From these results it can be seen that the point of largest deviation is at vertex 4 [θ ) (6.02, 6.02); see Figure 5] where M ) (2.98, 6.02). At this point the next linearization is performed. 5. Comparison of Linearizations: The parametric solution for the linearization at M ) (2.98, 6.02) is as follows. CRB1 : θ1 + θ2 e 9, θ1 g 3.98, θ2 g 3.98, M1 ) θ1, M2 ) θ2, where ProfitB1 ) 8.808θ1 + 8.388θ2 - 37.3518. CRB2 : θ1 + θ2 g 9, 3.98 e θ1 e 6.02, 3.98 e θ2 e 6.02, M1 ) 9 - θ2, M2 ) θ2, where ProfitB2 ) -10θ1 - 10.42θ2 + 131.9202. To determine which profile gives the tightest outer approximation to the nonlinear solution, we need only to compare those regions which overlap. In this case, this means comparing critical region 1A with 1B and critical region 2A with 2B. For critical region 1, solution A is best, while in critical region 2, part of the region is best approximated by A and part by B. This result is graphically shown in Figure 6. The hyperplane (a line in two dimensions) that separates critical region 2A from 2B is defined by
ProfitA2 - ProfitB2 ) (-10θ1 - 9.42θ2 + 126.4002) (-10θ1 - 10.42θ2 + 131.9202) ) 0 θ2 ) 5.52
Figure 7. Solution after six linearizations.
6. Results: Six further linearizations and comparisons are required to provide a linear profile within 0.3% of the nonlinear solution. Figure 7 shows the critical regions that give this linear approximation. The corresponding profit functions are
CRA1 :
ProfitA1 ) 8.408θ1 + 8.988θ2 - 39.2718
CRE1 :
ProfitE1 ) 8.408θ1 + 9.612θ2 - 42.0798
CRF1 :
ProfitF1 ) 8.200θ1 + 9.300θ2 - 39.8750
CRA2 :
ProfitA2 ) -10θ1 - 9.42θ2 + 126.4002
CRB2 :
ProfitB2 ) -10θ1 - 10.42θ2 + 131.9202
CRC2 :
ProfitC2 ) -10θ1 - 8.38θ2 + 121.7202
CRD 2:
ProfitD 2 ) -10θ1 - 9.92θ2 + 129.0352
CRF2 :
ProfitF2 ) -10θ1 - 8.90θ2 + 123.9250
Using this parametric profile of the optimal solutions, the evaluation of the expected cost for this design can be calculated. Because the value of the profit in terms of the uncertain parameters has been found, a mere function evaluation at the integration points provides this result. The Gaussian quadrature formula is used in this example to calculate the expected profit. Table 4 compares this solution with the primal problem using the pure stochastic method with a five-point Gaussian quadrature formula. The number of subproblems required for the solution is comparable between the two
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002 73 Table 4. Example 2: Hybrid/Stochastic Comparison stochastic
hybrid
no. of MINLPs no. of NLPs no. of mpLPs
1 25 0
1 23 6
expected profit (£)
29.697
29.697
CF(q1,q2) )
L L θU θU 1 - θ1 2 - θ2 wq1 wq2Jq1,q2 2 2
(12)
The resulting upper bound is £29.697, and thus the convergence is achieved. The design Y ) (0, 1) is found to be optimal. For this example the computational requirements of the hybrid and the stochastic approach are comparable. The requirements for the stochastic approach would have been higher if a large number of integration points would have been used or if a high tolerance for obtaining the parametric profiles would have been specified. 5.3. Example 3. The third example is a planning problem of the following form:
min Cost ) s.t.
peni(Di - Ai) + cBB ∑ i∈P
Ai e Di,
∑ i∈P
(13)
∀i∈P
Ai2 e B Ai g 0,
∀i∈P
15 e B e 25 Figure 8. Profit as a function of the uncertain parameter θ1 for different fixed values of θ2.
methods. Of course, the amount of information from the hybrid approach is much more extensive than just the expectancy of the profit. For instance, different cuts of the solution can be taken, so that one obtains a more intuitive idea of the behavior of the profitability with respect to the uncertain parameters. For example, Figure 8 shows the profit as a function of the uncertain parameter θ1 for different fixed values of θ2. 7. Formulating the Master Problem: The Lagrange multipliers from the parametric solution correspond exactly to those from the stochastic method. The master problem is formulated similarly to the stochastic master problem as follows:
min µ s.t. µ g
10Mq1 ,q ∑ q ,q 1
1
2
(11)
+ 12Mq21,q2
2
- 10(θq11,q2 - Mq11,q2) - 12(θq21,q2 - Mq21,q2) - 0.2(Mq11,q2)2 - 0.3(Mq21,q2)2 - 50(y1 + y2) -
λq1 ,q fq1 ,q (y1,y2) ∑ q1,q2 1
2
1
2
fq11,q2(y1,y2) ) Mq11,q2 + Mq21,q2 - 7y1 - 9y2 y1 + y2 e 1 y1, y2 ∈ {1, 0} Here λq1,q2 are the corrected Lagrange multipliers, where the correction factor is given by9
where Ai is the amount of product i produced, Di the amount of product i demanded, P the set of products 1-3, B the amount of raw material purchased, peni the penalty for unmet demand, and cB the cost price of B. The planning decision to be made is to calculate the amount of raw material B to be purchased given that the demands Di are normally distributed with mean µi ) (3.5, 4.5, 4.5) and standard deviation σi ) (0.162, 0.162, 0.162). The upper and lower bounds of the parameters taking into account 3.09σ are PU ) (4, 5, 5) and PL ) (3, 4, 4). The penalties are equal to the price of the products, where peni ) (0.4, 0.3, 0.2), and the cost price of B is £0.05/unit. The solution steps are as follows. 1. Initialization: B is the complicating variable. An initial design is selected by solving the optimization problem as an NLP in GAMS and defining the uncertain demands as variables between their upper and lower bounds. The optimal cost is found at the point: B ) 22, A1 ) 3, A2 ) 3, A3 ) 2, D1 ) 3, D2 ) 4, D3 ) 4. A tolerance of 0.5% is set to be the accuracy requirement for the linear parametric profiles. 2. First Primal Problem: The parametric results are obtained by using the linear underestimators constructed at various points on the nonlinear profile by using the dual information obtained from the solution of the NLPs at these points. For the first primal problem, two such linearizations are performed to meet the accuracy requirements. These regions are as follows. CRa1: D1 e 3.3, A1 ) 3, A2 ) 3, A3 ) 2, λ1,a ) 0.050, where Costa1 ) 0.3D1 + 0.3D2 + 0.2D3 - 1.1. CRa2: D1 g 3.3, A1 ) 3.48, A2 ) 2.61, A3 ) 1.74, λ2,a ) 0.057, where Costa2 ) 0.4D1 + 0.3D2 + 0.2D3 - 1.43. Using a five-point Gaussian quadrature integration, the expected cost is £2.28. This is the same as the expected cost by using the original stochastic formulation. When using cubatures to determine the expected profit, the expected value is £2.15. 3. First Master Problem: The master problem is formulated as follows:
74
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002
min µ s.t.
µg
(14)
peni(Dqi ,q ∑ ∑ q ,q i∈P 1
1
+ q1,q2 f1,a )
stochastic
q1,q2 - Ai,a ) + cBB
2
∑∑ q ,q c∈CR 1
2
Table 5. Example 3: Hybrid/Stochastic Comparison Gaussian quadrature
cubature
hybrid
251 0
29 0
25 4
no. of NLPs no. of mpLPs
q1,q2 λqc,a1,q2f1,a
2
q ,q 2 (Ai,a ) - B, ∑ i∈P 1
2
∀ q1, q2
15 e B e 25 The resulting lower bound is £2.27, and the new value of B returned is 25. 4. Second Primal Problem: For the value of B ) 25, the optimal value is found at the following point: A1 ) 3, A2 ) 3.328, A3 ) 2.219, D1 ) 3, D2 ) 4, D3 ) 4. The parametric results for the second primal problem are given by the following two critical regions. CRb1: D1 e 3.4, A1 ) 3, A2 ) 3.328, A3 ) 2.219, λ1,b ) 0.045, where Costb1 ) 0.27D1 + 0.3D2 + 0.2D3 - 1.002. CRb2: D1 g 3.4, A1 ) 3.71, A2 ) 2.36, A3 ) 1.86, λ2,b ) 0.054, where Cos tb2 ) 0.4D1 + 0.3D2 + 0.2D3 - 1.44. The expected cost by using a five-point Gaussian quadrature integration formula is given by £2.27, which is the same as the expected cost obtained by using the original stochastic formulation, whereas when using cubatures to determine the expected profit, the expected value is £2.14. 5. Second Master Problem: The master problem is formulated as follows:
min µ s.t. µ g
∑ ∑peni(Dqi ,q 1
(15) q1,q2 - Ai,a ) + cBB
2
q1,q2 i∈P
+
q ,q λqc,a,q f1,a ∑ ∑ q ,q c∈CR 1
1
µg +
2
1
peni(Dqi ,q ∑ ∑ q ,q i∈P 1
2
1
2
q ,q λqc,b,q f1,b ∑ ∑ q ,q c∈CR 1
2
1
q ,q 2 (Ai,a ) - B, ∑ i∈P
q1,q2 f1,b
∑ i∈P
)
2
q1,q2 2 (Ai,b )
2
q1,q2 - Ai,b ) + cBB
2
q1,q2 f1,a )
1
Table 6. Synthesis/Planning Model material balances process conversions
demand limitation availability limitations production limitations objective function
B 2 + B3 + B4 ) B A2 + A3 + A4 ) A C ) 0.9(Bf + B) B2 e MP2 ln(1 + A2/k2) B3 e MP3 ln(1 + A3/k3) B4 e MP4 ln(1 + A4/k4) C e DC A e AA Bf e A B A2 g 5 A3 e 20y3 A4 e 20y4 Profit ) PrCC - (PrAA + PrBBf) - [FC1 + 15(Bf + B)] - (FC2 + 5A2) - (FC3y3 + 15A3) - (FC4y4 + 5A4) - kPrC(DC - C)
Variables binary variables representing the existence (or not) of units 3 and 4 utilization of A at processes 2-4, respectively production of B at processes 2-4, respectively total production of B amount of material B purchased production of the final product
y3, y4 A2, A3, A4 B2, B3, B4 B Bf C
2
2
1
Figure 9. Synthesis problem illustration.
DC, AA, AB MP2, MP3, MP4 k2, k3, k4 PrC, PrA, PrB FC1, FC2, FC3 k
Parameters demand of C and the availability of A and B production constants constants prices of C, A, and B fixed cost constants penalty coefficient (k ) 1)
Table 7. Parameter Values
∀ q1, q2
- B, ∀ q1, q2
process
MP
1 2 3 4
18 20 15
k
FC
20 21 26
100 80 130 150
Prices
15 e B e 25 The resulting lower bound is £2.27, and therefore the algorithm converges. The optimal value of B is 25. The computational requirement for this example by using the hybrid parametric/stochastic approach and the original stochastic approach by using Gaussian quadratures are compared in Table 5. The computational requirement of the hybrid approach is comparable to that required when using the cubature integration formula, whereas it is much smaller than that required when using the quadrature integration scheme. 5.4. Example 4. This example is a process synthesis under uncertainty problem taken from ref 19. Chemical C is manufactured from material B (in process 1), which in turn can be purchased or produced from raw material A in processes 2 and 3 or 4. Figure 9 shows a schematic
PrC PrA PrB parameter DC AA AB
550 250 300 normal PDF (µ, σ2) N(20, 4) N(15, 1) N(15, 1)
of the process. The amounts of A and B available (AA and AB) are uncertain, as is the amount of C demanded (DC). The uncertainties are normally distributed. Table 6 shows the mathematical formulation of the problem, and Table 7 gives the parameter values for the problem. The solution steps are as follows. 1. First Iteration of GBD: The initial design is obtained by solving an MINLP. The solution is given by Y ) {1, 0}. A tolerance of 2% is initialized for linear
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002 75
Figure 10. Linear approximation results for the expected profit. Table 8. Critical Regions for the Synthesis Example: First Primal Y ) (1, 0) CRa1: Profita1 ) 74.855AB + 116.827DC + 233.508 -440.02AA - 600.15AB + 666.83DC - 574.56 e 0 74.855AB - 83.173DC + 596.958 e 0 -366.79AA - 600.15AB + 666.83DC - 1664.672 e 0 CRb1: Profitb1 ) 440.02AA + 675AB - 550DC + 808.07 -440.02AA - 600.15AB + 666.83DC - 574.56 g 0 440.02AA + 675AB - 750DC + 1171.52 e 0 73.23AA - 1090.11 e 0 CRc1: Profitc1 ) 200DC - 363.45 74.855AB - 83.173DC + 596.958 g 0 440.02AA + 675AB - 750DC + 1171.52 g 0 -366.79AA - 675AB + 750DC - 2261.63 e 0 CRd1: Profitd1 ) 366.79AA + 675AB - 550DC + 1898.18 -366.79AA - 600.15AB + 666.83DC - 1664.672 g 0 73.23AA - 1090.11 g 0 -366.79AA - 675AB + 750DC - 2261.63 g 0
approximations. At each linearization, a linear parametric overestimator is constructed by using information from a GAMS optimization problem at the point of linearization. In the first primal problem, four linearizations are required to satisfy the tolerance. The left-hand side of Figure 10 shows how each additional linearization brings the expected profit estimate closer to the expected profit when calculated using the original stochastic method. The results shown here are those found when using the five-point Gaussian quadrature and the cubature integration formulas. The final solution consists of four critical regions, for which the profit as a function of the uncertain parameters is presented in Table 8. The expected profit is calculated by using the cubature formula and the quadrature formula with values shown in Table 9. The master problem is formulated by using the dual information from the parametric solution, and the upper bound to the profit, using the cubature integration formula, is £3600.58.
Table 9. Expected Profit Values: Primal 1 £
original stochastic
hybrid
% diff
cubature quadrature
3459.62 3653.27
3470.58 3668.32
0.3 0.4
6
6
% diff
Table 10. Expected Profit Values: Primal 2 £
original stochastic
hybrid
% diff
cubature quadrature
3240.63 3411.39
3268.64 3442.18
0.9 0.9
5
5
% diff
From the pure stochastic method, the upper bound is £3585.28, a discrepancy of 0.4%. The new design returned by the master problem is Y ) {0, 0}. 2. Second Iteration of GBD: In the second iteration, six linearizations are performed to meet the required tolerance. Figure 10 shows how the linearizations bring the expected profit closer to the expectancy using the stochastic method. The expected profit is shown in Table 10. The new expected profit is lower than the previous lower bound, and therefore the original design is still found to be the best so far. The master problem returns a new upper bound of £3470.58, and therefore convergence is achieved. This agrees with the result from the master problem in the original stochastic formulation. The optimal design in both the original and the hybrid approach is therefore the first solution, Y ) {1, 0}. Figure 11 shows a cut of the solution profile, where the supply of B is at its nominal value, and the supply of A varies between its upper and lower bound. This is to illustrate the additional information supplied by the hybrid approach. It can be seen that, for smaller values of the demand DC, the design Y ) (0, 0) results in a higher profit than design Y ) (1, 0), even though overall Y ) (1, 0) is the optimal solution.
76
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002
Figure 11. Parametric profiles for fixed AB in the synthesis problem. Table 11. Subproblem Requirement Comparison pure stochastic MINLP MILP NLP
Gaussian quadrature
cubature
hybrid
1 2 250
1 2 28
1 2 61
3. Computational Requirements: Table 11 compares the computational requirements for the hybrid and the original stochastic formulations. It may be noted that the parametric profiles are obtained by using the dual information obtained from the solution of the NLPs at vertexes of the critical regions, and therefore no separate mpLPs are required in the solution method. While the hybrid method required fewer subproblems than the original stochastic method with quadrature integration, the original stochastic method is more efficient when the cubature integration formulas are used. However, this is, of course, a very strong function of the tolerance specified for obtaining the linear approximations. The discrepancy between the quadrature and cubature formulas is much larger than that between the hybrid and the original stochastic methods. Considering Figure 10 again, for the first problem after the first linearization, using the cubature integration method, the result already lies between the expected profits for the different integration methods of the original solution. In the second primal, only three (instead of six) linearizations seem to be required to achieve that result. The efficiency of the hybrid approach as compared to the original stochastic approach is thus shown to be highly dependent on the tolerance requirement imposed. One may also note that the integration schemes have some inherent inaccuracies which also depend on the number of integration points, whereas for the hybrid approach, the number of integration points used is insignificant with respect to the efficiency of the method. Therefore, a
method requiring more integration points but with higher accuracy might give a better result with a smaller number of subproblems. 6. Concluding Remarks Two-stage stochastic optimization can be used to calculate the expected profit. The computational requirements of this framework are quite high, and three approachessparallel computing, cubature integration formulas, and hybrid parametric/stochastic programmingshave been proposed to reduce the computational requirements. In this paper, a hybrid parametric/ stochastic algorithm for the solution of mixed-integer convex nonlinear optimization problems under uncertainty was presented. The basic idea of the algorithm is to employ parametric programming tools for evaluating the expected profit. This is achieved by deriving a complete profile of the profit as a function of uncertain parameters and then placing the integration points in the space of uncertain parameters so as to calculate the expected profit by simple function evaluations. The proposed approach performed better than the traditional approach, which employs Gaussian quadrature integration scheme, and its performance depends on the tolerance specified for obtaining the linear approximations of the parametric profile. A comparison of the hybrid parametric/stochastic programming approach with the approach employing cubature integration formulas was also presented. The performance of the hybrid approach was comparable to that employing cubature integration scheme. The hybrid approach is independent of the number of the integration points. Therefore, the hybrid approach will perform better than the techniques which use formulas which guarantee better accuracy but require a large number of integration points. The concepts were illustrated with ex-
Ind. Eng. Chem. Res., Vol. 41, No. 1, 2002 77
amples. The future work will focus on extending these concepts to problems involving nonconvexities. Acknowledgment Financial support from EPSRC (IRC Grant) is gratefully acknowledged. 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: Amterdam, 1990; pp 105-132. (2) Floudas, C. A. Nonlinear and mixed-integer optimization; Oxford University Press: New York, 1995. (3) Pistikopoulos, E. N. Uncertainty in process design and operations. Comput. Chem. Eng. 1995, 19, S553. (4) Pistikopoulos, E. N. Parametric and stochastic programming algorithms for process synthesis, design and optimization under uncertainty. Presented at Aspen World, Boston, MA, 1997. (5) Grossmann, I. E.; Sargent, R. W. H. Optimum design of heat exchanger networks. Comput. Chem. Eng. 1978, 2, 1-7. (6) Sahinidis, N. V.; Grossmann, I. E. Reformulation of multiperiod MILP models for planning and scheduling of chemical processes. Comput. Chem. Eng. 1991, 15, 255-272. (7) Bhatia, T. K.; Biegler, L. T. Multiperiod design and planning with interior point methods. Comput. Chem. Eng. 1999, 23, 919932. (8) van den Heever, S. A.; Grossmann, I. E. An iterative aggregation/disaggregation approach for the solution of a mixedinteger nonlinear oilfield infrastructure planning model. Ind. Eng. Chem. Res. 2000, 39, 1955-1971. (9) Pistikopoulos, E. N.; Ierapetritou, M. G. Novel approach for optimal process design under uncertainty. Comput. Chem. Eng. 1995, 19 (10), 1089-1110. (10) Diwekar, U. M.; Kalagnanam, J. R. Efficient sampling technique for optimization under uncertainty. AIChE J. 1997, 43, 440-447. (11) Acevedo, J.; Pistikopoulos, E. N. Stochastic optimization based algorithms for process synthesis under uncertainty. Comput. Chem. Eng. 1998, 22, 647-671. (12) Georgiadis, M. C.; Pistikopoulos, E. N. An integrated framework for robust and flexible process systems. Ind. Eng. Chem. Res. 1999, 38, 133-143. (13) Ahmed, S.; Sahinidis, N. V.; Pistikopoulos, E. N. An improved decomposition algorithm for optimization under uncertainty. Comput. Chem. Eng. 2000, 23, 1589-1604. (14) Acevedo, J.; Pistikopoulos, E. N. A multiparametric programming approach for linear process engineering problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717-728. (15) Dua, V. Parametric programming techniques for process engineering problems under uncertainty. Ph.D. Thesis, Imperial College, London, U.K., 2000. (16) Hene´, T. S. A hybrid parametric stochastic programming approach for mixed-integer nonlinear optimization problems under uncertainty. Master’s Thesis, Imperial College, London, U.K., 1998.
(17) Acevedo, J.; Pistikopoulos, E. N. A hybrid parametric/ stochastic programming approach for mixed-integer linear problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 22622270. (18) Epperly, T. G. W.; Ierapetritou, M. G.; Pistikopoulos, E. N. On the global and efficient solution of stochastic batch plant design problems. Comput. Chem. Eng. 1997, 21, 1411-1431. (19) Bernardo, F. P.; Pistikopoulos, E. N.; Saraiva, P. M. Integration and computational issues in stochastic design and planning optimization problems. Ind. Eng. Chem. Res. 1999, 38, 3056-3068. (20) Ierapetritou, M. G.; Pistikopoulos, E. N. Global optimization for stochastic planning, scheduling and design problems. In Global Optimization in Engineering Design; Grossmann, I. E., Ed.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; pp 231-287. (21) Fiacco, A. V. Introduction to Sensitivity and Stability Analysis in Nonlinear Programming; Academic Press: New York, 1983. (22) Gal, T. Postoptimal analyses, parametric programming, and related topics; de Gruyter: New York, 1995. (23) Acevedo, J.; Pistikopoulos, E. N. A parametric MINLP algorithm for process synthesis problems under uncertainty. Ind. Eng. Chem. Res. 1996, 35, 147-158. (24) 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. (25) Papalexandri, K. P.; Dimkou, T. I. A parametric mixed integer optimization algorithm for multi-objective engineering problems involving discrete decisions. Ind. Eng. Chem. Res. 1998, 37 (5), 1866-1882. (26) Acevedo, J.; Pistikopoulos, E. N. An algorithm for multiparametric mixed integer linear programming problems. Oper. Res. Lett. 1999, 24, 139-148. (27) Dua, V.; Pistikopoulos, E. N. Algorithms for the solution of multiparametric mixed-integer nonlinear optimization problems under uncertainty. Ind. Eng. Chem. Res. 1999, 38, 3976-3987. (28) Dua, V.; Papalexandri, K. P.; Pistikopoulos, E. N. A parametric mixed-integer global optimization framework for the solution of process engineering problems under uncertainty. Comput. Chem. Eng. 1999, 23, S19-22. (29) Viswanathan, J.; Grossmann, I. E. A combined penalty function and outer-approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14 (7), 769-782. (30) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: a users guide; Scientific Press: Washington, D.C., 1996. (31) Acevedo, J. Parametric and stochastic programming algorithms for process synthesis under uncertainty. Ph.D. Thesis, Imperial College, London, U.K., 1996.
Received for review January 17, 2001 Revised manuscript received October 3, 2001 Accepted October 16, 2001 IE0100582
