A Parametric Mixed-Integer Optimization Algorithm ... - ACS Publications

1866

Ind. Eng. Chem. Res. 1998, 37, 1866-1882

A Parametric Mixed-Integer Optimization Algorithm for Multiobjective Engineering Problems Involving Discrete Decisions Katerina P. Papalexandri* and Theodora I. Dimkou FORTHsChemical Process Engineering Research Institute, 57001 Thessaloniki, Greece

An iterative, decomposition-based algorithm is proposed in this paper, for the solution of convex parametric MINLP problems and, in extension, the identification of the noninferior solution set in multiobjective problems involving continuous and discrete decisions. The parametric optimal solution is constructed via an upper- and lower-bound procedure. Parametric upper bounds are identified with the solution of parametric NLP problems, while the parametric lower bound is updated in each iteration via the solution of a parametric MILP Master problem, which involves only the binary variables of the initial problem. Convergence properties and computational requirements are discussed in example problems from process synthesis under uncertainty and simultaneous product/process design. 1. Introduction In the contemporary process industry conventional process economics, i.e., capital and operating cost, have become only part of the directives in process decision making. Environmental impact, robust and flexible production planning, process controllability, product performance, etc., formulate distinct objectives, acting usually in a conflicting manner. Within a multiobjective decision framework, the best-compromise solution strongly depends on the relative importance of the different objectives to the decision maker (DM). It will, however, belong to the set of noninferior solutions, where one objective can only improve at the expense of other objectives. Research work on multiobjective optimization focuses on either (i) the identification of the best-compromise solution interactively with the decision maker (e.g., Loganathan and Sherali, 1987; Jahn and Merkel, 1992) or (ii) the identification of the noninferior solution curve and its properties (Khahn, 1992; Maeda, 1994; Li, 1996). Reviews of multicriteria engineering optimization problems and solution approaches can be found in Clark and Westerberg (1983) and Sawaragi et al. (1985). Most recent systems engineering applications of multiobjective optimization methods appear in process design and control (Meadowcroft et al., 1992; Luyben and Floudas, 1994), waste treatment (Ciric and Huchette, 1993; Ciric and Jia, 1994), and concurrent process and product design (Yoshimura and Kondo, 1996). When the multiobjective problem involves continuous and discrete decisions, as is usually the case in process decision making, it can be formulated as follows:

min [f1(x,y), f2(x,y), ..., fn(x,y)] x,y

s.t.

h(x,y) ) 0 g(x,y) e 0

(P1)

x ∈ X, y ∈ Y ) {0, 1}m * Author to whom correspondence should be addressed. E-mail: [email protected]. Telephone: +30-31-980 149. Fax: +30-31-980 180.

Then, the identification of the noninferior solutions can be formulated as a multiparametric mixed-integer nonlinear programming problem (pMINLP):

min f1(x,y) x,y

s.t.

h(x,y) ) 0 g(x,y) e 0 f2(x,y) e p2 l fn(x,y) e pn

(P)

x ∈ X, y ∈ Y ) {0, 1}m up lo up p2 ∈ [plo 2 , p2 ], ..., pn ∈ [pn , pn ]

where the parameter ranges [plo, pup] result from the optimization of each scalar objective separately, as described in section 2. The parametric programming theory has mostly focused on linear problems (Gal, 1979). Extensions to the nonlinear and mixed-integer case are limited. Fiacco (1983) and Fiacco and Ishizuka (1990) review algorithms for nonlinear parametric optimization. Special cases of mixed-integer parametric optimization have been addressed by McBride and Yorkman (1980), by Skorin-Kapov and Granot (1987), and recently by Acevedo and Pistikopoulos (1997). The proposed algorithms for parametric MINLP problems have been limited to the scalar parametric case. They are based on decomposing the continuous from the discrete optimization part, employing principles from the outer approximation method (OA) of Kocis and Grossmann (1987). Pertsinidis (1992) in his Ph.D. thesis presented two algorithms for the solution of RHS-parametric MINLP problems: (i) an extension

S0888-5885(97)00720-3 CCC: $15.00 © 1998 American Chemical Society Published on Web 04/16/1998

Ind. Eng. Chem. Res., Vol. 37, No. 5, 1998 1867

of the Jenkins (1982) algorithm for linear problems and (ii) an alternative algorithm based on the outer approximation method. Focusing on waste treatment issues within process synthesis problems, Ciric and Huchette (1993) proposed an iterative algorithm to determine the sensitivity of maximum net profits of a process to changes in the waste treatment cost through the solution of a twoobjective optimization problem. In an outer optimization step, candidate discrete regions of the noninferior solution curve are identified, solving an outer approximation of the net-profit minimization problem (MILP Master) at one point of the waste-treatment-cost range. For each discrete vector, the continuous noninferior curve is found by applying the sequential approximation method of Ciric and Jia (1994). The overall noninferior solution curve is, thus, updated in an iterative manner, until the Master problem becomes infeasible. Acevedo and Pistikopoulos (1996) extended the work of Pertsinidis and proposed the solution of a parametric MILP Master problem to construct a lower bound to the parametric solution of a pMINLP and reduce, thus, the size of the discrete search space. Their algorithm is based on the construction of upper- and lower-bound curves to the parametric solution, alternating between a parametric NLP problem (upper bound) and a parametric MILP problem which is an outer approximation of the original problem (lower bound). As linearizations of the full model at several parameter points are accumulated in the MILP Master problem, this may become of considerable size. In this work an alternative Master problem is proposed for a decomposition-based identification of the noninferior solution set through a parametric optimization problem. Lagrangian information from a parametric continuous optimization step is utilized, similarly to the generalized Benders decomposition method for the scalar MINLP case (Geoffrion, 1972; Paules and Floudas, 1989), to derive candidate discrete vectors and discrete regions of the set. As only Lagrangian-cut constraints are considered in the mixed-integer optimization step (Master problem), smaller-in-size problems are, in general, required to be solved. In the sequel, the theoretical developments are first presented and an iterative algorithm is proposed for the solution of RHSscalar parametric MINLP problems and, in extension, the identification of the noninferior solution curve of two-objective MINLP problems (which are the most common in process synthesis and design). Extensions to the multiparametric/multiobjective case are discussed. Numerical examples from process synthesis and simultaneous product/process design are solved to illustrate the applicability and the computational requirements of the method.

2. Mathematical Developments The noninferior solution set of the multiobjective MINLP problem (P1) can be identified from the para-

metric solution of problem (P):

z(p) ) min f1(x,y) x,y

h(x,y) ) 0

s.t.

g(x,y) e 0 f2(x,y) e p2 l fn(x,y) e pn

(P)

x ∈ X, y ∈ Y ) {0, 1}m up lo up p2 ∈ [plo 2 , p2 ], ..., pn ∈ [pn , pn ]

where

plo i ) min fi(x,y) x,y

s.t.

h(x,y) ) 0 g(x,y) e 0 x ∈ X, y ∈ Y ) {0, 1}m

and

pup i ) max (pii′) i′*i

pii′ ) fi(x*,y*) with s.t.

fi′(x*,y*) ) min fi′(x,y) x,y

h(x,y) ) 0 g(x,y) e 0 x ∈ X, y ∈ Y ) {0, 1}m

for i ) 2, ..., n, the objectives of the initial problem (P1). Let us define the following:

p ) [p2, ..., pn] parameter vector F ) [f2, ..., fn]

function vector of the parametric constraints

up lo up P ) [plo 2 , p2 ] × ... × [pn , pn ] parameter space

In the following analysis we assume that (P1) and (P) are convex and y appears linearly in the objective functions and constraints; i.e., the following conditions hold: (i) f1, h, g, and F are linear with respect to y. In most process system problems, binary variables are used linearly in process models (e.g., to denote the existence of a unit). (ii) h is affine with respect to x. (iii) f1, g, and F are continuously differentiable and convex with respect to x.

1868 Ind. Eng. Chem. Res., Vol. 37, No. 5, 1998

(iv) The set

Let us define

Zy,p ) {z ) [z1z2] ∈ Rq: h(x,y) ) 0, g(x,y) e z1, F(x,y) - p e z2, for some x ∈ X} is closed for each fixed y ∈ Y and p ∈ P. (v) For each fixed y ∈ Y and p ∈ P, the resulting problem (P) (a) has a finite solution and an optimal multiplier vector for the equalities and the inequalities or (b) is unbounded (i.e., its objective function value goes to -∞). The main idea of the proposed algorithm is to decompose the initial pMINLP problem (P) into two subproblems: (a) A parametric NLP problem (Primal problem, P2) which results from fixing the binary vector y and providing, thus, a parametric upper bound to the optimal solution of (P). (b) A parametrix MILP (Master problem), which involves only the binary variables y and is a relaxation of (P), so that a valid lower bound to z(p) is obtained. a. Parametric Upper Bound. Fixing the binary vector y, the solution of the resulting parametric NLP problem provides clearly an upper bound to z(p):

x

x

g(x,y) e 0 F(x,y) e p x∈X and V ) {(y, p): h(x,y) ) 0, g(x,y) e 0, F(x,y) - p e 0 for some x ∈ X}. v(y,p) is the projection of (P) on the (p,y) space. (P3) can be written as

z(p) ) min v(y,p) y

(y, p) ∈ (Y × P) ∩ V

s.t.

inf f1(x,y) x

(P2)

F(x,yj) e p

s.t. h(x,y) ) 0 ) v(y,p) ) g(x,y) e 0 F(x,y) e p x∈X

x∈X

[ sup

inf L(x,y,p,µ,λ,τ)], (y, p) ∈ (Y × P) ∩ V

µ,λg0,τg0 x∈X

p∈P

where

Based on the convexity assumptions and for a given tolerance , zup(p) can be found in a series of upper- and lower-bound approximations, as a piecewise linear function of the parameters (Fiacco, 1983). The solution procedure for the RHS-parametric case is discussed in Appendix A. While the theoretical background for the solution of (P2) is provided for the general case, this solution can be efficiently obtained (in terms of computational requirements) only for a small number of parameters. b. Parametric Lower Bound. The derivation of the Master problem is based on nonlinear duality theory, similar to the Master problem of the generalized Benders decomposition (Geoffrion, 1972; Paules and Floudas, 1989). Problem (P) can be written as follows:

z(p) ) min inf f1(x,y) y

L(x,y,p,µ,λ,τ) ) f1(x,y) + µTh(x,y) + λTg(x,y) + τT(F(x,y) - p) Due to the strong duality theorem (Bazaraa et al., 1993), problem (P) is equivalent to

z(p) ) min y

sup

inf L(x,y,p,µ,λ,τ)

µ,λg0,τg0 x∈X

z(p) ) min µB y,µB

µB g inf L(x,y,p,µ,λ,τ), ∀ µ, ∀ λ g 0, x∈X

∀ τ g 0 (M)

h(x,y) ) 0 (y, p) ∈ (Y × P) ∩ V

g(x,y) e 0

x ∈ X, y ∈ Y ) {0, 1}m p∈P

(P′)

Using the definition of supremum as the lowest upper bound and introducing a scalar µB, we have

s.t.

x

F(x,y) e p

(P4)

[ ]

h(x,yj) ) 0 g(x,yj) e 0

s.t.

h(x,y) ) 0

s.t.

We consider the dual representation of v(y,p) (Bazaraa et al., 1993):

zup(p) ) min f1(x,yj) s.t.

v(p,y) ) inf f1(x,y)

(P3)

Problem (M) (Master problem) is equivalent to (P) but involves an infinite number of constraints. Considering only the constraints that correspond to specific values of µ, λ g 0 and τ g 0 (i.e., dropping out constraints), we obtain a relaxation of (M) that when


solved parametrically for p will provide a parametric lower bound to the optimal solution of (P):

The corresponding relaxed Lagrangian-cut constraints are considered in the relaxed Master problems (Paules and Floudas, 1989), which can be written as

zlo(p) ) min µB

zlo(p) ) min µB

y,µB

y,µB

s.t.

kT

k

kT

k

k

µB g f1(x ,y) + µ h(x ,y) + λ g(x ,y) + T

τk (F(xk,y) - p) (M′) (y, p) ∈ (Y × P) ∩ V where xk, µk, λk, and τk correspond to solutions of problem (P2) for specific yk and pk. c. Infeasible Primal Subproblems. Problem (P2) may be infeasible for a fixed yj, in the whole parameter range P or part of P. In the general multiparametric case, the parameter region R ⊆ P, where (P2) is feasible, can be found from the solution of a series of parametric problems, as discussed in Appendix A. For the scalar parametric case, the p-range where (P2) is feasible, [p1, p2], is given by

p1 ) min θ, x,θ

s.t.

p2 ) max θ x,θ

h(x,yj) ) 0 g(x,yj) e 0 F(x,yj) e θ x∈X θ ∈ [plo, pup]

Note that, from the definition of pup, (P2) is feasible at pup, when the scalar optimization problem

min f1(x,yj) x

s.t.

h(x,yj) ) 0 g(x,yj) e 0

is feasible for yj and then only the minimization problem is required. For the part of P where (P2) is infeasible, a relaxed problem can be solved:

min (s1 + s2 + R)

x,s1,s2,R

h(x,yj) + s1 - s2 ) 0 g(x,yj) e R F(x,yj) - p e 0 x∈X s1, s2, R g 0 p ∈ [plo, pup]

T

T

µB g f1(xk,y) + µk h(xk,y) + λk g(xk,y) + T

τk (F(xk,y) - p)

(k feasible p-point, p ∈ R)

µTk h(xk,y) + λTk g(xk,y) + τTk (F(xk,y) - p) e 0 (k infeasible p-point, p ∈ P) (M1) y∈Y Clearly, with the accumulation of Lagrangian constraints ((P2) solutions for fixed (y, p)) in (M1), zlo(p) becomes tighter (increases). Note that problem (RP2) is not solved parametrically, as only p-points are of interest, to generate relaxed Lagrangian-cut constraints. In the general case, (M1) is a RHS-multiparametric MILP, involving only binary variables and µB (see for solution methods Appendix B). Furthermore, for each point solution of (P2) (or (RP2)) one Lagrangian constraint is considered. This results in a much smaller problem compared to the Master problem proposed by Acevedo and Pistikopoulos (1996), which involves linearizations of all the problem constraints and both the continuous and binary variables. Due to this difference, as in the case of deterministic MINLPs, the OA-based Master problem may predict tighter lower bounds. However, specific properties of the Lagrangian cuts (see Remarks in section 3) may lead to better convergence. In general, the trade-off between problem size and quality of bounds depends on the particular problem. The theoretical developments for the parametric upper and lower bounds are valid for the general multiparametric case. However, as the number of objectives increases, the computational requirements for the solution of the subproblems (P2) and (M1) increase significantly. In the following sections, the scalar parametric case (two-objective MINLPs) will be examined in detail. 3. Identification of Noninferior Curve in Two-Objective MINLPs

x∈X

s.t.

s.t.

(RP2)

The bounding properties of problems (P2) and (M1) can be exploited to identify the parametric solution of an MINLP problem (and, in extension, the noninferior solution curve of the initial multiobjective problem (P1)). For the two-objective case, the iterative algorithm of Figure 1 is proposed, including the following steps: Step 0. Initialization. (a) Solve each scalar-optimization problem (MINLP) separately to identify the extreme points of the noninferior solution curve and the range of the second objective function (initial parameter interval). The initial best upper bound is zup(p) ) fup 1 . (b) Define a tolerance and an initial binary vector y0 (e.g., the optimal binary vector for the first objective). Step 1. Primal Problem. In each iteration k, for each one of the current sets of parameter intervals, i, and the corresponding integer vector, yk,i, do the following:


(a) Determine the feasible parameter range of problem (P2) (see Appendix A). (b) Solve the parametric NLP (within tolerance ), to determine (i) z* k,i as a piecewise linear function of p and (ii) the optimal solution and multiplier vectors for parameter points in the interval. (c) Compare z*k,i(p) to the current upper bound and update zup k (p). (d) If (P2) is infeasible in some part of the interval, solve the relaxed problem (RP2) at infeasible parameter points, to obtain the multiplier vector for the relaxed Lagrangian constraints. Step 2. Master Problem. For each one of the current sets of parameter intervals, do the following: (a) Formulate the Master problem, adding the Lagrangian constraints from the last Primal step. Note that, for convex problems, the Lagrangian constraints exclude integer vectors that have already been examined. However, integer cuts may be included in the Master problem for those vectors, to avoid repetition due to numerical discrepancies in the exact calculation of (especially large) Lagrange multipliers. (b) Solve the updated parametric MILP Master problem. If (M1) is infeasible in some part of the interval, the optimal solution has been found in this part and is given from the current parametric upper bound and the corresponding vector. Step 3. Convergence Check. (a) For each interval resulting from the solution of (M1), compare the solution of (M1) to the current upper up lo - zk,i e in some part of the interval, bound. If zk,i then the optimal solution has been found in this part up and the corresponding vector. and is given from zk,i (b) Check if convergence has been achieved in the whole parameter range. If yes, stop. Otherwise, go to (c). (c) Update the set of intervals where the optimal solution has not been found and the candidate binary vectors. Go to Step 1. In the initialization step, the parameter bounds can be set arbitrarily (without solving the scalar MINLPs) and updated in each major iteration, solving separately the scalar NLPs for the current integer vector. In step 2b, as Lagrangian constraints are added in each iteration, the obtained solution of the Master problem is a monotonic nondecreasing sequence of lower bounds. The convergence criterion (comparison of bounds) results from the bounding properties of steps 1 and 2. The algorithm will converge in a finite number of iterations, since the number of discrete alternatives is finite. Remarks. (1) From the Kuhn-Tucker optimality conditions of the Primal problem (P2) and for a fixed vector yj, we have

∇xf1(x*,yj) + µ*T∇xh(x*,yj) + λ*T∇xg(x*,yj) + τ*∇xf2(x*,yj) ) 0 λ*Tg(x*,yj) ) 0 τ*(f2(x*,yj) - p) ) 0 It results that the optimal multiplier vector [µ*, λ*, τ*] will be a continuous function of p, in each p-interval where the set of active constraints is constant. These intervals can be identified by solving for the maximum

Figure 1. Iterative algorithm for the solution of pMINLP.

p-range, where a set of active constraints defines the optimal solution (Kuhn-Tucker conditions). In those intervals where the parametric constraint is not active, the corresponding Lagrangian constraint will not be a function of p, and one solution point (Lagrangian constraint) suffices in the Master problem to exploit sensitivity information. In general, including in the Master problem Lagrangian constraints from p-points where the active constraints change leads usually to good lower bounds and candidate vectors. (2) Similar to the observation of Sahinidis and Grossmann (1991), from the strong duality theorem at the optimal solution of (P), at a p-point we have

L(x*,y*,µ*,λ*,τ*,p) ) f1(x*,y*,p) and also (saddle-point theorem):

L(x*,y*,µ*,λ*,τ*,p) e L(x*,y,µ*,λ*,τ*,p), ∀ y ∈ [0, 1] If at the k major iteration of the algorithm (see Figure 1) the pNLP is solved for y ) y* in a p-interval and Lagrangian cuts for p-points of the interval are included in the Master problem

zlo(p) g L(x*,y*,µ*,λ*,τ*,p) ) f1(x*,y*,p) Thus, the algorithm will converge in k-iteration for the corresponding p-interval or part of the interval where the Lagrangian cuts have been considered.


Exploiting this, we can start the algorithm (a) for an interval around plo with the optimal y-vector of the f1-optimization problem and (b) for an interval around pup with the optimal y-vector of the scalar optimization problem that corresponds to pup. Thus, we may “cut-off” (i.e., find the optimal binary vector) parts of the parameter range in the first iteration. (3) Extending the algorithm of Figure 1 to the multiobjective case requires the solution of multiparametric problems at steps 2 and 3. Due to convexity, in each parameter region (hyperspace) the solution of the Primal problem (upper bound) can be approximated by a set of hyperplanes, which will also give a valid upper bound to the overall solution of the mpMINLP. The solution of the mpMILP Master problem will also be given by a set of hyperplanes. The efficient and exact solution of multiparametric NLP and MILP problems, as well as efficient comparison of parametric bounds, remains an open issue. In Appendices A and B, some special cases are examined and valid approximations are suggested that require, however, intensive computations.

Solving the MINLP problems for f1 and f2 separately, we find

flo 1 ) -57

f2 ) 329 ) fup 2

flo 2 ) -0.593 75

for y* ) (0, 0, 0)

f1 ) -1.279 95 ) fup 1 for y* ) (0, 0, 1)

The noninferior solution curve will be the solution with respect to p of the following problem:

z(p) ) min x12 - x2 + x3 + 3y1 + 2y2 + y3 x,y

s.t. f2 ) 2x12 + x32 - 3x1 + x2 - 2y1 + y2 - 2y3 e p g1 ) 3x1 - x2 + x3 + 2y1 e 0 g2 ) 4x12 + 2x1 + x2 + x3 - 40 + y1 + 7y2 e 0 g3 ) -x1 - 2x2 + 3x3 + 7y3 e 0 g4 ) -x1 - 10 + 12y1 e 0

4. Examples

g5 ) x1 - 5 - 2y1 e 0 Example 1. Consider the following two-objective MINLP problem:

g6 ) -x2 - 20 + y2 e 0

min [f1, f2]

g7 ) x2 - 40 - y2 e 0

x,y

g8 ) -x3 - 17 + y3 e 0

where

f1 ) x12 - x2 + x3 + 3y1 + 2y2 + y3 f2 ) 2x12 + x32 - 3x1 + x2 - 2y1 + y2 - 2y3

g9 ) x3 - 25 - y3 e 0 y1, y2, y3 ∈ {0, 1} p ∈ [-0.593 75, 329]

s.t. g1 ) 3x1 - x2 + x3 + 2y1 e 0 g2 ) 4x12 + 2x1 + x2 + x3 - 40 + y1 + 7y2 e 0 g3 ) -x1 - 2x2 + 3x3 + 7y3 e 0 g4 ) -x1 - 10 + 12y1 e 0 g5 ) x1 - 5 - 2y1 e 0 g6 ) -x2 - 20 + y2 e 0 g7 ) x2 - 40 - y2 e 0 g8 ) -x3 - 17 + y3 e 0 g9 ) x3 - 25 - y3 e 0 y1, y2, y3 ∈ {0, 1}

Applying the -constraint method (Clark and Westerberg, 1983), f* 1(p) is found (see Figure 2a) by solving 40 MINLPs. Applying the algorithm of Figure 1 we have the following: Iteration 1. Let us start with a y1 ) (0, 1, 1). The corresponding pNLP is infeasible for p ∈ [-0.593 75, 0.406 25]. For p ∈ [0.406 25, 329] a first approximation zj1 to zup 1 (see Figure 2a) is found from

p1 ) 0.406 25 pup ) 329

z* 1 ) 1.0781 z* 2 ) -54

zj1 ) -0.6176p + 1.1462 LB1 ) -68207p + 27710

LB2 ) -54

An intermediate point pint ) 0.407 057 is defined from LB1 ) LB2. Continuing the partitioning procedure, zup 1 (p) is found as illustrated in Figure 2a (see Table 1 for detailed expressions).


Figure 3. Iterations 2 and 3 in example 1.

For the Master problem we examine the following cases: (a) We solve the Master problem considering only the Lagrangian constraints at plo and pup. Then, Figure 2. Iteration 1 in example 1.

The Lagrangian constraints for the Master problem are lo

at p ) -0.593 75 0.667(1.0729 - 2y1 + y2 2y3 - p) + 0.333(-6.333 + 7y3) e 0 (relaxed) at p1 ) 0.406 25 -1.922 + 3y1 + 2y2 + y3 + 2302080(1.40624 - 2y1 + y2 - 2y3 - p) + 1151040(-6.999 + 7y3) e µ at pup ) 329

-57 + 3y1 + 2y2 + y3 + (1 - y2) + (-1 + y3) e µ

zlo 1 ) -57 and the candidate vector is y2 ) (0, 0, 0). (b) We find the p-points where the set of active constraints of the pNLP problem changes (see Table 1) and consider all the corresponding Lagrangian constraints. The obtained lower bound is illustrated in Figure 2b and is clearly a tighter one (see Table 1 for detailed expressions). The candidate vector is again y2 ) (0, 0, 0). (c) We consider the Lagrangian constraints from an equal number of p-points that lie within one interval, where the active constraints remain the same. The obtained lower bound (see Figure 2b) is close to that of

{


Table 1. Parametric Bounds in Example 1

fup 1 (p) )

{

-58.2130p + 24.7271 -10.1330p + 5.1555 -3.6927p + 2.4321 -1.0275p + 1.0636 -0.94966p + 0.8836 -0.29627p - 22.3999 -0.1190p - 32.1446 -0.0393p - 42.3480 -54.0

Iteration 1. y* ) (0, 1, 1) p ∈ [0.40625, 0.4070] p ∈ [0.4070, 0.4228] p ∈ [0.4228, 0.4813] p ∈ [0.4813, 2.3120] p ∈ [2.3120, 35.6335] flo 1 (p) ) p ∈ [35.6335, 55] p ∈ [55, 128.01] p ∈ [128.01, 296] p ∈ [296, 329]

{

-p - 1 -0.0610p - 41.6603 -0.03125p - 46.7521 -57

p ∈ [-0.59375, 43.3041] p ∈ [43.3041, 170.8564] p ∈ [170.8564, 327.9338] p ∈ [327.9338, 329]

-1.0129p - 1.5233 -0.7658p - 10.6372 -0.42754p - 23.55995 -0.06105p - 40.78247 -0.03124p - 45.8147 -0.03125p - 45.68956

p ∈ [-0.59375, 36.88252] p ∈ [36.88262, 38.20101] p ∈ [38.20111, 46.99312] p ∈ [46.99322, 168.85634] p ∈ [168.85644, 172.85634] p ∈ [172.85644, 329]

active sets at 0.75, 29.75, 101.2082, 106.8423, 295.8715, 296

fup 2 (p) )

fup 3 (p) )

{

-0.7865p - 0.2059 -113.3737p - 28.3505 -51.9769p - 13.0139 -29.0339p - 7.2936 -15.8187p - 4.0186 -8.5695p - 2.2563 -4.5396p - 1.3399 -2.3004p - 0.9440 -1.0021p - 0.9273 -0.4257p - 20.1866 -0.0475p - 41.3668

Iteration 2. y* ) (0, 0, 0) p ∈ [-0.25, -0.24998] p ∈ [-0.24998, -0.24979] p ∈ [-0.24979, -0.24933] p ∈ [-0.24933, -0.24786] p ∈ [-0.24786, -0.24304] flo 2 (p) ) p ∈ [-0.24304, -0.22741] p ∈ [-0.22741, -0.17680] p ∈ [-0.17680, -0.01285] p ∈ [-0.01285, 33.41078] p ∈ [33.41078, 56] p ∈ [56, 329]

{

{

Iteration 3. y* ) (0, 0, 1) for p ∈ [-0.59375, 39.604] and y* ) (0, 0, 0) for p ∈ [39.604, 329] -71.9041p - 43.6150 p ∈ [-0.59375, -0.59321] -p - 1 p ∈ [-0.59375, 33.8825] -11.4365p - 7.7450 p ∈ [-0.59321, -0.58013] -0.7658p - 8.9340 p ∈ [33.8826, 35.8825] -3.9719p - 3.4145 p ∈ [-0.58013, -0.52576] -0.7658p - 8.4020 p ∈ [35.8826, 38.4647] -1.1514p - 1.9316 p ∈ [-0.52576, 0.45098] -0.4273p - 21.4150 p ∈ [38.4648, 44.3625] flo 3 (p) ) -0.9698p - 2.0136 p ∈ [0.45098, 39.604] -0.1205p - 34.3252 p ∈ [44.3626, 104.3531] -0.4257p - 20.1866 p ∈ [39.604, 56] -0.0610p - 40.5382 p ∈ [104.3532, 172.8563] -0.0475p - 41.3668 p ∈ [56, 329] -0.0312p - 45.6895 p ∈ [172.85644, 329]

case b around the considered points and looser in the rest of the p-range. The candidate vector is y2 ) (0, 0, 0). (d) We consider the Lagrangian constraints corresponding to more than one p-point in each interval with the same active constraints. The obtained lower bound is illustrated in Figure 2b and is not much tighter than the one in case b. Naturally, the more Lagrangian constraints we include in the Master problem, the tighter the obtained bound. However, we see that the p-points where the set of active constraints changes (and, thus, the function of the Lagrangian constraint with respect to p) suffice to give good lower bounds. Since up zlo 1 (p) < z1 (p) ∀ p ∈ [-0.593 75, 329]

convergence is not achieved in the first iteration. Iteration 2. The parametric NLP problem for y2 ) (0, 0, 0) yields a better upper bound for the whole p-range; thus, the “best” current structure is updated (see Figure 3a and Table 1). The second Master problem (including the Lagrangian constraints for the intermediate points where the set of active constraints changes) shows that (see Figure 3a) (i) in [-0.593 75, 61.388] convergence of bounds is not achieved and y3 ) (0, 0, 1) is the new candidate binary vector for the interval and (ii) in [61.388, 329] convergence is achieved and y2 ) (0, 0, 0) is the optimal binary vector (see remark 2 in section 3). Iteration 3. The pNLP problem is solved for y3 ) (0, 0, 1) in [-0.593 75, 61.388]. The solution is proven to be a better upper bound in the interval [-0.593 75, 39.604], while for [39.604, 61.388] y2 ) (0, 0, 0) remains the “best” binary vector (see Figure 3b). The parametric upper bound is updated accordingly (see Table 1).

up In the third Master problem, we find zlo 3 (p) g z3 (p), ∀ p ∈ [-0.593 75, 329], and overall convergence is achieved. The optimal parametric solution (i.e., the noninferior solution curve of the initial two-objective MINLP problem) is the current upper bound and is illustrated in Figure 3b. Note that a discontinuity appears at p* ) 39.604. The optimal solution has been found in three major iterations, solving 36 NLPs, 6 MILPs, and 3 pLPs. Applying the OA-based algorithm of Acevedo and Pistikopoulos (1996), for the same initial binary vector, the optimal parametric solution is found in four major iterations (38 NLPs, 8 MILPs, and 4 pLPs). Furthermore, if we start the algorithm with y0 ) (0, 0, 1) or y0 ) (0, 0, 0), i.e., the optimal binary vector at one of the p-bounds, the noninferior solution curve is found in two major iterations (see remark 2 in section 3). Example 2. Process Synthesis/Planning under Uncertainty. The optimal process configuration for a desired product and the corresponding production specifications clearly depend on the market requirements for the product and the extent to which the decision maker wants to cover the market demand. Changes in market requirements affect not only the optimal production rate (raw material consumption, etc.) but also discrete decisions about the process, e.g., shutting down a particular process. Thus, the trade-offs between process cost/profit and market demands are part of a mixed-integer optimization problem. The process superstructure of Figure 4 is considered (Kocis and Grossmann, 1987; Acevedo and Pistikopoulos, 1996), which comprises four processes used to produce C. Raw material A is transformed into B through processes P1, P2, and/or P3. A restricted amount of fresh B can be added to complement production requirements, and the mixture is further processed

1874 Ind. Eng. Chem. Res., Vol. 37, No. 5, 1998 Table 2. Detailed Model in Example 2 min C ) 250A + 300Bf + (5A1 + 15A2 + 5A3 + 15B) 550C + 80y1 + 130y2 + 150y3 + 100 s.t.

Figure 4. Superstructure for example 2.

in P4 to form C. The aim is to analyze the effect of changes in market demand (D) for product C on the process flowsheet and production specifications (i.e., a combination of processes, C production rate, and raw material consumption) that minimize a process cost function (or equivalently maximize a net profit function). The synthesis/planning problem is modeled as a scalar parametric MINLP problem. Binary variables

Figure 5. Parametric upper and lower bounds in example 2.

C ) 0.9B B1 e 18 ln(1 + A1/20) B2 e 20 ln(1 + A2/21) B3 e 15 ln(1 + A3/26) A ) A1 + A2 + A3 B ) B1 + B2 + B3 + Bf Bf e 15 CeD A1 - 25y1 e 0 A2 - 20y2 e 0 A3 - 20y3 e 0 A1, A2, A3 g 0

are employed to denote the existence of P1, P2, and P3 in the final flowsheet. The detailed model is given in Table 2 and involves material balances, production restrictions for each process, and design constraints. The

Ind. Eng. Chem. Res., Vol. 37, No. 5, 1998 1875 Table 3. Parametric Bounds in Example 2a Iteration 1 best upper bound for (1, 1, 1) -495.915p + 8243.261 zup ) -209.458p + 2513.837

pNLP solution for y1 ) (1, 1, 1) -495.915p + 8243.261 z* ) -209.458p + 2513.837

p ∈ [17, 20] p ∈ [20, 25]

pMILP solution-lower bound -511.688p + 8361.067 zlo ) -200p + 2127.380

p ∈ [17, 20] p ∈ [20, 25]

pNLP solution for y2 ) (1, 1, 0) -515.527p + 8502.289 z* ) -210.283p + 2419.361

Iteration 2 best upper bound for (1, 1, 0) p ∈ [17, 19.928] -515.527p + 8502.289 zup ) p ∈ [19.928, 25] -210.283p + 2419.361


p ∈ [17, 20] p ∈ [20, 25]

{ { { {

{

proposed structure (1, 1, 0)

{

pNLP solution for y3 ) (1, 0, 1)

p ∈ [17, 19.928] p ∈ [19.928, 25]

proposed structure (1, 0, 1) Iteration 3 best upper bound for (1, 1, 0) -515.527p + 8502.289 zup ) -210.283p + 2419.361

z* ) -200p + 3147.046

p ∈ [17, 25]


p ∈ [17, 20.182] p ∈ [20.182, 25]

{

p ∈ [17, 20] p ∈ [20, 25]

{

p ∈ [17, 19.928] p ∈ [19.928, 25]

objective function is a net cost function, involving the fixed and operating cost for each process, cost of raw materials, and the revenue from product C. A fixed flow of A is assumed (A ) 33) and market demand in C varies

17 e D e 25 The algorithm of Figure 1 is applied to solve the parametric MINLP, starting with

plo ) 17 pup ) 25 y1 ) (1, 1, 1) In iteration 1, the obtained upper and lower bounds are shown in Figure 5a. A new structure is proposed, y2 ) (1, 1, 0) (see Table 3 for detailed expressions). The parametric upper bound is updated in the Primal step of iteration 2 (see Figure 5b and Table 3). The lower bound obtained at the second Master pMILP problem is very close to the upper bound (see Figure 5b). With a convergence tolerance ) 1%, we find that, in [20, 25], y2 ) (1, 1, 0) is proven the optimal structure, while for [17, 20], y3 ) (1, 0, 1) is proposed. The upper bound is not improved in iteration 3 and y2 ) (1, 1, 0) remains the best structure. In the Master step, we find zlo(p) > zup(p) and y2 (i.e., the selection processes P1 and P2) is proven to be the optimal structure for the whole p-range. The parametric solutions are shown in Figure 5c (see Table 3 for detailed expressions). In the strict case, when < 1%, the optimal parametric curve requires the solution of 8 NLPs, 9 MILPs, and 3 pLPs and 2.8 CPU s on a DEC A500 workstation. The required number of NLPs also depends on the given tolerance for the approximation of the nonlinear functions with piecewise linear functions. Note that the lower bound obtained from the pMILP Master problem, after the optimal structure has been examined (iteration 2), is very close to the best upper bound. With a sufficient convergence tolerance to alleviate numerical problems in the calculation of Lagrange multipliers, we can “cut-off” a significant part


of the parametric range in iteration 2 (see also remark 2 in section 3). Example 3. Process Synthesis/PlanningsII. Similar to example 2, a synthesis/planning problem is examined (Kocis and Grossmann, 1987; Acevedo and Pistikopoulos, 1996), with two desired products PR1 and PR2. The process superstructure, involving eight processes, is illustrated in Figure 6. The detailed model, with mass balances, yield and availability constraints, bounds, and production specifications, is given in Table 4. We assume that the two product demands vary simultaneously:

D1 ) 14 - 7p D2 ) 18.5 - 9p p ∈ [0, 1] and analyze the effect of p on the optimal total process cost (where revenue from the products is also accounted for) and process configuration. The resulting pMINLP is solved by employing the algorithm of Figure 1, starting with an initial configuration

y1 ) (1, 1, 1, 1, 1, 1, 1, 1) that involves all processes.

1876 Ind. Eng. Chem. Res., Vol. 37, No. 5, 1998 Table 4. Detailed Model in Example 3a MINLP Model objective function mass balances

yield relations desired production availability flow bounds

cost ) -1775PR1 - 1500PR2 + (200A + 550Bf + 300C) + (5IS12 + 18IS22 + 7IS32 + 15IS42 + 8IS52 + 18IS62 + 17IS72 + 10IS82 + (100y1 + 800y2 + 6800y3 + 1000y4 + 500y6 + 1100y7 + 1000y8) A ) IS1 + IS2 + IS3 OS1 + OS2 + OS3 + Bf ) IS4 C ) IS5 + IS6 OS5 + OS6 ) IS8 + (IS7 - OS4) OSi ) PciISi i ) {4, 7, 8} OSi e Pci ln(1 + ISi/ki) i ) {1, 2, 3, 5, 6} PR1 g D1(p), PR2 g D2(p) Bf e 30 ISi - MIiyi e 0 ∀ i process D1 ∈ [7, 14], D2 ∈ [9.5, 18.5] Model Constants

P1 P3 P5 P7 a

Pc

k

MI

18 25 22 0.65

20 26 21

40 35 55 60

P2 P4 P6 P8

Pc

k

MI

20 0.9 28 0.9

21

30 80 75 600

25

Where Pci and ki ) yield constants and MIi ) maximum allowable input to process.

Table 5. Optimal Parametric Solution in Example 3 interval

cost vs p

(y1

[0, 0.0853] [0.0853, 0.1710] [0.1710, 0.3670] [0.3670, 0.6544] [0.6544, 0.7770] [0.7770, 1]

-9368.1918p - 339.6257 -8052.7031p - 451.8497 -5945.2753p - 812.3487 -3983.3350p - 1532.4863 -5679.2705p - 618.5626 -3696.4540p - 2159.2346

optimal structure y3 y4 y5 y6 y7

y2

y 8)

0

1

0

1

1

1

1

1

0

1

0

1

0

1

1

1


Figure 7. Optimal solution in example 3.

The optimal parametric solution is found in 14 major iterations (with a convergence tolerance = 0), requiring the solution of 137 NLPs, 56 MILPs, and 15 pLPs and 44.8 CPU s on a DEC A500 workstation, and is illustrated in Figure 7. The configuration involving processes P2 and P4-P8 is the optimal one for p e 0.6544 (i.e., D1 g 9.42 and D2 g 12.11); for 0.6544 e p e 1 (i.e., when market demands drop below 35% of the expected decrease), process P5 becomes redundant and the configuration involving P2, P4, and P6-P8 is more economical. The corresponding cost functions are given in Table 5. Example 4. Simultaneous Product and Process Design. A custom-made composite product is manufactured with the impregnation of chemicals with

desired properties into a porous wood matrix. Consider the manufacturing alternatives in Figure 8 (Ozyurt et al., 1996). Alternative monomers may be used as impregnants within either aqueous mixtures (LS) or nonaqueous solvents (SS), e.g., CO2. The solventimpregnant mixture penetrates the porous material in a conventional liquid impregnation step (process P4) or at supercritical conditions (process P5) that allow for higher mass-transfer rates, selective distribution of chemicals, etc. The solvent is recycled, and the impregnants, which include monomers, are further treated in a polymerization step, which may be aided by radiation (process P6) or catalysts (process P7). Two alternative monomers (either B or C) can be used with the nonaqueous solvent, which must be preprocessed (processes P2 and P3, respectively). The nonaqueous solvent is produced from A (e.g., propane in the case of CO2) in P1 or bought. A simplified process model is given in Table 6. The maximum mechanical strength of the final composite product is of interest. Clearly, this depends on

Ind. Eng. Chem. Res., Vol. 37, No. 5, 1998 1877 Table 6. Process Model in Example 4 min cost )

100y1 + 0.3is1 + 1500y2 + 150is2 + 1200y3 + 160is3 + 800y4 + 50ls + 200(p4 - 1) + 10t4 + 1200y5 + 5ss + 200(p5 - 1) + 15t5 + 1000y6 + 700is6 + 800y7 + 200is7 x 6 + x7

max X ) s.t. mass balances

mechanical strength

production spec’s process bounds

os1 - 0.7is1 ) 0 os2 - 0.85is2 ) 0 os3 - 0.65is3 ) 0 os4 - 0.85is4 ) 0 rs2 - 0.4os4 ) 0 rs2 + ls - is4 ) 0 os5 - 0.99is5 ) 0 rs1 - 0.2os5 ) 0 os1 + rs1 + os2 + os3 + ss - is5 ) 0 is6 + is7 - os4 - os5 ) 0 os6 - 0.97is6 ) 0 os7 - 0.9is7 ) 0 x4 - 2.75 ln(p4) + 3/(450 + t4) - 1500(1 - y4) e 0 x4 - 20y4 e 0 x5 - 2.98 ln(p5) + 3/(450 + t5) - 1500(1 - y5) e 0 x5 - 150y5 e 0 -x6 - x7 + 2.8 e 0 x6 - 25y7 e 0 x7 - 7y7 e 0 x6 - 0.95(x4 + x5) e 0 x7 - 0.8(x4 + x5) e 0 -os6 - os7 + 30 e 0 is1 - 1000y1 e 0 is2 - 1000y2 e 0 is3 - 1000y3 e 0 is4 - 1000y4 e 0 is5 - 1000y5 e 0 is6 - 1000y6 e 0 is7 - 1000y7 e 0 20y4 e t4 e 100y4 1 + y4 e p4 e 5y4 + 1 20y5 e t5 e 200y5 1 + 30y5 e p5 e 150y4 + 1

a Where is ) inlet flow to process i, os ) outlet flow from i, rs ) recycle flow, x ) mechanical strength of composite out of i, p , t ) i i i i i pressure and temperature, and ss, ls ) solvent flows.

Table 7. Optimal Parametric Solution in Example 4 interval

cost vs p

[2.8, 3.0489] [3.0489, 3.937] [3.937, 5.402] [5.402, 8.610] [8.610, 8.618] [8.618, 11.953] [11.953, 16.717] [16.717, 16.81] [16.81, 17.255] [17.255, 19.472] [19.472, 21.118] [21.118, 21.432] [21.432, 21.862]

344.4757p + 9312.0261 488.7835p + 8994.0013 82.99p + 10728.396 181.957p + 10193.75 73333.8065p - 619644.5 17301.6835 0.7236p + 17293.0542 2227.7831 - 19937.0038 2555.737p - 25450.068 10.549p + 18472.1456 307.086p + 12697.85 671.124p + 5009.873 2400.394p - 32051.96

the manufacturing conditions, selection of chemicals, entrainers, etc., that affect also the total process cost. In order to exploit process-product trade-offs, the process configuration and the final product (product specifications) are optimized simultaneously. The problem is formulated as a two-objective MINLP. Assuming convex relations between mechanical strength and temperature and pressure, a simplified convex model is considered (see Table 6). The noninferior solution curve is calculated as the parametric solution of the cost optimization problem with respect to final mechanical strength (x). The parameter bounds are calculated through the solution of the single MINLP problems:

Xlo ) 2.8 and Xup ) 25.873

optimal structure y2 y3 y4 y5 y6

(y1

y7)

0

0

0

1

0

0

1

0

0

0

1

0

1

1

0

0

1

0

1

1

1

0

0

1

1

1

1

1

Starting with an initial configuration involving processes P3-P7, the optimal parametric solution (noninferior curve) is found in seven major iterations of the proposed algorithm (see Figure 9 and Table 7) requiring 14 CPU s on a DEC A500 workstation. As the specifications for final mechanical strength improve, the total process cost increases and processes with better characteristics are selected. Note that for a X* ) 8.618 there is a drastic increase in the total cost, as process P5 becomes necessary. If the marketable value of mechanical strength is close to that point, the decision maker can trade-off product performance with significant savings in manufacturing cost (nearly 28.75% cost reduction by not selecting P5). Example 5. To illustrate the computational requirements of the proposed algorithm in MINLP problems


Figure 9. Optimal solution in example 4.

with more than two objectives, consider the following problem:

min [f1, f2, f3] x,y

where

f1 ) x12 + x22 - 10x1 - x2 - y1 - 2y2 f2 ) (4x12 + 3x22 - x1 - 5x2 - y1 + 10y2 - 10)/3 f3 ) (2x12 + 7x1 - 14x2 + 2y1 + 2y2 - 6)/2 s.t.

-x1 + 3x2 - y1 + 0.5 e 0 y1, y2 ∈ {0, 1}

Solving each scalar MINLP problem, we find

flo 1 ) -28.25

Figure 10. Iteration 1 in example 5.

The following Master problem is formulated considering Lagrangian-cut constraints at p-points of P:

min µ

fup 1 ) 5.2021

flo 2 ) -4.021 36

fup 2 ) 30.75

flo 3 ) -3.506 95

fup 3 ) 37.5

The parametric solution of the corresponding pMINLP is sought in the parameter space:

(f2, f3) ∈ P ) [-4.021 36, 30.75] × [-3.506 95, 37.5] Starting with y1 ) (1, 1), the feasible region of the first Primal NLP is the solution of a scalar parametric problem (see Appendix A). A good approximation is given by (see Figure 10b)

R1 ) [-0.544 38, 30.75] × [-2.718 35, 37.5] The parameter space R1 is discretized into triangles (see Appendix A and Figure 13) and the parametric solution is found as a piecewise linear function of f2 and f3 (see Figure 10a):

z1,i ) ai + bif1 + cif3 For a sufficient approximation of zup 1 (f2,f3), the solution of 180 NLPs is required.

µ,y1,y2

s.t.

µ g -36.9158 - 7.4378y1 + 30.6483y2 (f2, f3) ∈ R1 9.7945f2

µ g -17.39564 - 14.0453y1 + 2.8048y2 9.60965f3 (f2, f3) ∈ R1 µ g -24.75 - y1 - 2y2 µ g -18.2270 - 0.9118y1 - 1.4523y2 - 0.1438f2 (f2, f3) ∈ R1 0.1361f3 µ g -18.52 - 0.7830y1 - 1.8915y2 0.2169f3

(f2, f3) ∈ R1

µ g -19.039 - 1.1y1 - 0.9912y2 0.3026f2

(f2, f3) ∈ R1

-4.02136 - 0.33y1 + 3.33y2 - f2 e 0

(f2, f3) ∈ P

-0.6520 - 0.4y1 + 0.15y2 - 0.3f3 e 0

(f2, f3) ∈ P

y1 + y2 e 1 y1, y2 ∈ {0, 1}


Note that the application of the -constraint method (Clark and Westerberg, 1983) results in the same solution requiring 400 MINLPs. 5. Concluding Remarks

Figure 11. Iteration 2 in example 5.

and is a two-parametric MILP problem with overlapping critical regions of the pLP subproblems; thus, the algorithm of Acevedo and Pistikopoulos (1997) cannot be employed. zlo is approximated by employing a discretization scheme (see Appendix B), involving the solution of 10 MILPs. The obtained lower bound is illustrated in Figure 10a. The first iteration requires overall 15.4 CPU s on a DEC A500 workstation. Comparison of the current upper and lower bounds (see also Figure 10a) shows that y1 ) (1, 1) is the optimal binary vector in the space where the maximum inscribed rectangle is defined by

7.279 e f2 e 30.75 7.333 e f3 e 37.5 For the rest of P, divided into three new rectangles, we consider y2 ) (1, 0), as the candidate vector (resulting as the candidate vector in the larger part of P). The second pNLP is solved similarly, and z*2 is approximated with the solution of 99 NLPs, as another set of planes (see Figure 11a). Comparison of zup 2 shows that (1, 0) is a better 1 and z* binary vector in the shaded region shown in Figure 11b (defined by a set of parametric constraints). up The second pMILP Master problem gives zlo 2 g z2 in the whole parameter space. The second iteration requires 11.9 CPU s.

The solution of parametric MINLP problems arises in multiobjective optimization problems and trade-off analysis, where continuous and discrete decisions are simultaneously considered to exploit better process potential. An iterative algorithm is proposed in this paper for the solution of a class of such problems (convex case) and, in extension, the identification of discontinuous, convex noninferior solution sets in multiobjective problems. An upper- and lower-bound procedure is developed, without requiring the repetitive solution of MINLP problems. Decomposing the continuous from the discrete optimization part, the lower-bound curve is identified with the solution of a pMILP problem that involves only binary variables and a small number of constraints (one Lagrangian cut for each point solution of the Primal problem). The mathematical background is developed for the general multiparametric case, where the parameters appear in the model in a convex way. Computationally efficient approximation methods are given for the scalar parametric case (two-objective problems). For the multiparametric problem presented in this work, the proposed method, although computationally intensive, still required the solution of less NLP and MILP problems than earlier methods. However, the efficient generation and comparison of upper and lower bounds are challenging questions and are the subject of our current work. The convexity assumption limits the application of the proposed algorithm. Still, a large class of problems can be addressed, as in high-level decision making simple process models are usually employed that can be convexified with proper mathematical manipulations. However, the extension of global optimization techniques to the nonconvex parametric case remains a major research challenge. Acknowledgment Financial support from the European Communities, under Contract JOE3-CT95-0017 (IDEES), is gratefully acknowledged. Appendix A. Solution of pNLP Problems a. Scalar Parametric Case. Consider the problem (P2), for a fixed vector yj and one parameter p:

zup(p,yj) ) min f1(x,yj) x

s.t.

h(x,yj) ) 0 g(x,yj) e 0 f2(x,yj) e p x∈X p ∈ [plo, pup]

(S1)


Figure 13. Partitioning of the two-parameter region. Figure 12. Scalar parametric NLP bounding procedure.

(a) Solving (S1) at pi and pi+1 ∈ [plo, pup], to find z*i and z*i + 1, respectively, the connecting line segment

b. Multiparametric Case. Consider problem (P2) for the general multiparametric case:

zup(p) ) min f1(x,yj) x

zji(R,yj) ) Rz* i + (1 - R)z* i+1, R ∈ (0, 1) is a valid parametric upper bound to zup(p,yj), since (S1) is convex. (b) A lower bound to zup(p,yj) is given by

s.t.

h(x,yj) ) 0 g(x,yj) e 0 F(x,yj) e p

zi(p,yj) ) max (LBi, LBi+1)

(P2)

x∈X

where

p∈P LBi ) z* j) + ∇pz* j)(p - pi) i (pi,y i (pi,y

LBi+1 ) z* j) + ∇pz*i (pi+1,yj)(p - pi+1) i+1(pi+1,y ∇pz* is evaluated through the Lagrange multipliers of the parametric constraints. The above bounds can be further sharpened by defining an interior point pint ∈ [plo, pup] from

LBi ) LBi+1 and repeating steps a and b in the subintervals [plo, pint] and [pint, pup] (see Figure 12). Thus, for a given tolerance , a sufficient approximation to the solution of (S1) can be found by zji’s, dividing the initial range [plo, pup] into subintervals, where

zji - zi e The convergence properties of such a procedure are given by Fiacco and co-workers (Fiacco, 1983; Fiacco and Ishizuka, 1990). The approximation is a piecewise linear function of p. Due to the convexity assumptions, it will be an upper bound to zup(p,yj) and, thus, a valid upper bound to the solution of the initial parametric MINLP problem (P1).

where P is a hyperspace that can be sufficiently approximated by a set of hyperrectangles. Solving (P2) at K + 1 vertices of a hyperrectangle pi (with K parameters), from the (z*i,k, pi,k) solutions, we can obtain a hyperplane

zji ) ai +

∑k bkipk

that, due to the convexity assumptions, is a valid upper bound to zup(p). In each one of the considered vertices, the tangent hyperplanes

zi,k ) zji,k * +

∑k ∇p zji,k* (pk - pi,k) k

define valid lower bounds for zup(p). These upper bounds, zji, can be refined by discretizing further the ith hyperrectangle. For the two-parameter case a partitioning scheme is illustrated in Figure 13, where sufficient approximation of zup(p1,p2) in each triangle is based on comparison of zji and the actual solution of (P2) at the center of each triangle. c. Feasible Region of Primal pNLP Problems. For a fixed binary vector yj, let R be the p-space, where the corresponding Primal pNLP (P2) is feasible. We distinguish two cases:


the RHS of inequality constraints (λk, τk g 0), whereas the scalar variable µ is the only continuous variable:

(i) The scalar optimization problem

min f1(x,yj)

zlo(p) ) min µ

x

µ,y

h(x,yj) ) 0

s.t.

s.t.

g(x,yj) e 0

y ∈ {0, 1}m

(P0)

p∈P

x∈X is infeasible. Then problem (P2) is also infeasible and R ) 0. (ii) Problem (P0) is feasible for some x ∈ X. up Then, from the definition of pup 2 , ..., pn , problem (P2) up will be feasible for some x ∈ X at pk ) [pup 2 , ..., pn ]. The region R ⊆ P, where (P2) is feasible, is the convex region defined from the parametric solutions of the following NLP problems:

Fi(p′) ) max δ x,δ

h(x,yj) ) 0

s.t.

Aµ + By e D + Rp

where the entries of matrix A are ai ∈ {0, -1}. a. Scalar Parametric Case. For scalar parametric problems, the algorithm of Pertsinidis (1992) can be employed. Starting from the upper or lower parameter bound, the optimal binary vector is updated by scanning the parameter range along a single direction: (a) Solve the deterministic MILP at plo to find a binary vector y′. (b) Fixing the binary vector y′, solve the resulting scalar pLP (Gal, 1979). The parametric solution z′(p) represents an upper bound to zlo(p). (c) Find a minimum p-value for which a new binary vector renders a better objective value than the one predicted by the pLP, solving the following MILP problem:

g(x,yj) e 0

pbp ) min θ x,y,θ

fi(x,yj) e pup i - δ fi′(x,yj) e pi′,

(P2)

i′ ) 2, ..., n, i′ * i

lo 0 e δ e pup i - pi

solved for i ) 2, ..., n. An underestimation of R can be found from the following problem, as the maximum hyperrectangle within P, where (P2) is feasible:

x,δi

s.t.

∑i δi

h(x,yj) ) 0 g(x,yj) e 0 F(x,yj) e pup i - δi

i ) 2, ..., n

x∈X lo 0 e δi e pup i - pi

Ax + By e Rθ c1Tx + c2Ty e z′(θ) y ∈ {0, 1}m

x∈X

max

s.t.

i ) 2, ..., n

Note that the above problem is similar to the flexibility index evaluation problem of Swaney and Grossmann (1985). Appendix B. Solution of pMILP Problems In the proposed pMILP Master problem (problem (M1)) the optimization parameters appear linearly in

θ∈P (d) If the above problem is feasible, update y′ and repeat steps a-c in the remaining interval [pbp, pup] until the whole parameter range is scanned. b. Multiparametric Case. For mpMILPs Acevedo and Pistikopoulos (1997) have proposed an extension of the branch & bound (B&B) method to the multiparametric case, utilizing the algorithm for the solution of multiparametric LP problems in each node (Gal and Nedoma, 1972): (a) Solve the fully relaxed problem as a mpLP to provide a lower bound to the optimal solution. (b) Select the branching variables and develop the B & B enumeration tree. At each node: (i) Solve a mpLP. (ii) Compare the parametric solution to the current upper bound for the parameter space under consideration. (iii) Discard a parameter region if (1) mpLP is infeasible (2) an integer solution is found (3) the parametric solution is greater than the upper bound for this region. (c) Update the parameter space, where the optimal solution is not found and the mpLP is not proven infeasible, and the corresponding parametric upper and lower bounds. Continuing with the B & B tree (by fixing binary variables), repeat steps b and c until in each parameter region the optimal solution is found or mpLP is infeasible.


A major requirement of the algorithm is that critical regions of the resulting mpLPs are not overlapping. When there are overlapping critical regions, the algorithm may not identify all the solutions. Alternatively, in the general case, in order to identify all the solutions, the parameter space can be discretized into hyperrectangles (see Figure 13), with the same binary vector proposed in each vertex, zlo(p) is found from the solution of mpLP in each hyperrectangle. Obviously, in such a case, the number of required MILPs may increase significantly. Literature Cited Acevedo, J.; Pistikopoulos, E. N. A Parametric MINLP Algorithm for Process Synthesis Problems under Uncertainty. Ind. Eng. Chem. Res. 1996, 35, 147-158. Acevedo, J.; Pistikopoulos, E. N. A multiparametric programming approach for linear process engineering problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717-728. Bazaraa, M. S.; Sherali, H. D.; Shetty, C. M. Nonlinear Programming: Theory and Algorithms; John Wiley: New York, 1993. Ciric, A. R.; Huchette, S. G. Multiobjective Optimization Approach to Sensitivity Analysis: Waste Treatment Costs in Discrete Process Synthesis and Optimization Problems. Ind. Eng. Chem. Res. 1993, 33, 2636-2646. Ciric, A. R.; Jia, T. Economic Sensitivity Analysis of Waste Treatment Costs in Source Reduction Projects: Continuous Optimization Problems. Comput. Chem. Eng. 1994, 18, 481495. Clark, P. A.; Westerberg, A. W. Optimization for Design Problems having more than one Objective. Comput. Chem. Eng. 1983, 7, 259-283. Fiacco, A. V. Introduction to Sensitivity and Stability Analysis in Nonlinear Programming; Academic Press: New York, 1983. Fiacco, A. V.; Ishizuka, Y. Sensitivity and Stability Analysis for Nonlinear Programming. Ann. Oper. Res. 1990, 27, 215-237. Gal, T. Post-Optimal Analyses: Parametric Programming and Related Topics; McGraw-Hill: London, U.K., 1979. Gal, T.; Nedoma, J. Multi-Parametric Linear Programming. Mang. Sci. 1972, 18, 406. Geoffrion, A. M. Generalized Benders Decomposition. JOTA 1972, 10, 237-260. Jahn, J.; Merkel, A. Reference Point Approximation Method for the Solution of Bicriterial Nonlinear Optimization Problems. JOTA 1992, 74, 87-103. Jenkins, L. Parametric Mixed Integer Programming: An Application to Solid Waste Management. Mang. Sci. 1982, 28, 12701284.

Khahn, P. Q. Proper Solutions of Vector Optimization Problems. JOTA 1992, 74, 105-130. Kocis, G. R.; Grossmann, I. E. Relaxation Strategy for the Structural Optimization of Process Flowsheets. Ind. Eng. Chem. Res. 1987, 26, 1869-1880. Li, D. Convexification of a Noninferior Frontier. JOTA 1996, 88, 177-196. Loganathan, G. V.; Sherali, H. D. A Convergent Interactive Cutting-Plane Algorithm for Multiobjective Optimization. Oper. Res. 1987, 35, 365-377. Luyben, M. L.; Floudas, C. A. Analyzing the Interactions of Design and Control. 1. A Multiobjective Framework and Application to Binary Distillation Synthesis. Comput. Chem. Eng. 1994, 18, 933-969. Maeda, T. Constraint Qualifications in Multiobjective Optimization Problems: Differentiable Case. JOTA 1994, 80, 483-500. McBride, R. D.; Yorkman, J. S. Finding all the Solutions for a Class of Parametric Quadratic Integer Programming Problems. Mang. Sci. 1980, 26, 784-795. Meadowcroft, T. A.; Stephanopoulos, G.; Brosilow, C. The Modular Multivariable Controller. 1. Steady-State Properties. AIChE J. 1992, 38, 1254-1278. Ozyurt, B.; Mogili, P.; Mierau, B.; Sunol, S. G.; Sunol, A. K. A Hierarchical Approach to Simultaneous Design of Products and Processes. Comput. Chem. Eng. 1996, 20, S73-S78. Paules, G. E.; Floudas, C. A. APROS: Algorithmic Development Methodology for Discrete Continuous Optimization Problems. Oper. Res. 1989, 37, 902-915. Pertsinidis, A. On the parametric optimization of mathematical programs with binary variables and its applications in the chemical engineering process synthesis. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA, 1992. Sahinidis, N. V.; Grossmann, I. E. Convergence Properties of Generalized Benders Decomposition. Comput. Chem. Eng. 1991, 15, 481-491. Sawaragi, Y.; Nakayama, H.; Tanino, T. Theory of Multiobjective Optimization; Academic Press: London, 1985. Skorin-Kapov, J.; Granot, F. Nonlinear Integer Programming: Sensitivity Analysis for Branch and Bound. Oper. Res. Lett. 1987, 6, 269-274. Swaney, R. E.; Grossman, I. E. An Index for Operational Flexibility in Chemical Process Design. AIChE J. 1985, 31, 621-630. Yoshimura, M.; Kondo, H. Product Design based on Concurrent Processing of Design and Manufacturing Information by Utility Analysis. Concurrent Eng. Res. Appl. 1996, 4, 379-388.

Received for review October 14, 1997 Revised manuscript received February 4, 1998 Accepted February 5, 1998 IE970720N

A Parametric Mixed-Integer Optimization Algorithm ... - ACS Publications

Recommend Documents