Nonconvex Generalized Benders Decomposition with Piecewise

Article pubs.acs.org/IECR

Nonconvex Generalized Benders Decomposition with Piecewise Convex Relaxations for Global Optimization of Integrated Process Design and Operation Problems Xiang Li, Yang Chen, and Paul I. Barton* Process Systems Engineering Laboratory, Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, United States ABSTRACT: This paper considers the global optimization of challenging stochastic or multiperiod mixed-integer nonconvex programs that arise from integrated process design and operation. The difficulties of the problems are large scale and the nonconvexity involved. Recently, a novel decomposition method, called nonconvex generalized Benders decomposition (NGBD), has been developed to solve this problem to global optimality finitely, and this method shows dramatic computational advantages over traditional branch-and-bound based global optimization methods because it can exploit well the decomposable structure of such problems. Since the convergence rate of NGBD is largely dependent on the tightness of the convex relaxations of the nonconvex functions, the efficiency of NGBD can be improved by generating tighter convex relaxations. Building on the success of piecewise linearization for bilinear programs in the process systems engineering literature, this paper develops a piecewise convex relaxation framework, which can yield tighter convex relaxations for factorable nonconvex programs, and integrates this framework into NGBD to expedite the solution. Case studies of a classical literature problem and an industry-level problem show that, while NGBD can solve problems that are intractable for a state-of-the-art global optimization solver, integrating the proposed piecewise convex relaxation into NGBD helps to reduce the solution time by up to an order of magnitude.

■

INTRODUCTION Many optimal process design and operation problems can be posed mathematically as mixed-integer nonlinear programming (MINLP) problems. When optimizing the design decisions and the operational decisions simultaneously, the MINLP model may exhibit the following structure:

associated with xh for scenario h. The global optimization of Problem P is usually very challenging because of its large scale (when s is large) and the nonconvex functions involved. As a nonconvex MINLP, Problem P can be solved rigorously by general-purpose deterministic global optimization methods, such as branch-and-reduce,5 SMIN-αBB, and GMIN-αBB,6 and nonconvex outer approximation7 (which is an extension of traditional outer approximation methods8,9 to programs with nonconvex functions participating). However, when the number of scenarios addressed in the problem increases, the solution times with these methods may increase dramatically,10,11 because the sizes and/or the numbers of the subproblems to be solved in these methods increase significantly. The solution of Problem P can be significantly expedited by exploiting its decomposable structure. Problems with such structures, such as stochastic linear/convex programs, can be decomposed via Benders decomposition (BD)12 (also called the L-shaped method in the stochastic programming literature13) and generalized Benders decomposition (GBD),14 into a sequence of subproblems whose sizes are independent of the number of scenarios addressed. The solution of these subproblems leads to a sequence of nonincreasing upper bounds and a sequence of nondecreasing lower bounds which converge to an optimum of the problem. BD and GBD are based on transforming the original problem into a dual space and this transformation is rigorous only when strong duality holds. Therefore, these methods cannot guarantee a global solution

s

∑ wh(chTy + fh (xh)) x ,..., x , y min s

1

s.t.

h=1

gh(xh) + Bh y ≤ 0,

xh ∈ Xh ,

∀ h ∈ {1, ..., s},

∀ h ∈ {1, ..., s},

y∈Y

(P)

where Xh = {xh ∈ Πh ⊂  : ph(xh) ≤ 0}, Y = {y ∈ {0,1} : Ay ≤ d}, Πh is convex, and functions f h: Πh → , gh: Πh → m and ph: Πh → mp are continuous. It is assumed at least one function in Problem P is nonconvex. The binary variables y represent the design decisions, such as whether or not to build a unit or pipeline in the system. The continuous variables xh represent the operational decisions, such as long-term operating conditions of the system, for scenario h. A scenario may represent a possible realization of uncertain parameters whose outcome is known only after the design decisions are made; then wh denotes the probability of the outcome of each scenario and Problem P is usually referred to as a stochastic MINLP.1,2 Or a scenario may be associated with a set of parameter values in an operating period which are known to vary over different operating periods; then wh denotes the frequency of the occurrence of each operating period and P is often referred to as a multiperiod MINLP.3,4 The objective of Problem P is a cost function, which is usually associated with an economic metric, over the scenarios, where cTh y represents the cost associated with y for scenario h and f h(xh) represents the cost nx

ny

© 2012 American Chemical Society

Received: Revised: Accepted: Published: 7287

June 13, 2011 February 28, 2012 April 11, 2012 April 11, 2012 dx.doi.org/10.1021/ie201262f | Ind. Eng. Chem. Res. 2012, 51, 7287−7299

Industrial & Engineering Chemistry Research

Article

where Dh = {(xh,eh) ∈ Πh × Θh: up,h(xh,eh) ≤ 0, uq,h(xh,eh) ≤ 0}, Θh is convex, and functions uf,h: Πh × Θh → , ug,h: Πh × Θ h → m, up,h: Πh × Θ h → mp, and ue,h: Πh × Θh → me are convex on Π h × Θ h. So Problem LBP is a convex MINLP or a mixed-integer linear program (MILP). This problem involves auxiliary variables eh and constraints ue,h(xh,eh) ≤ 0 that may be required to construct smooth relaxations. Several convex relaxation techniques are available to generate Problem LBP, e.g., McCormick relaxation and outer linearization24 for factorable nonconvex functions, which usually introduce additional variables and constraints for differentiable relaxations, and αBB for twicedifferentiable nonconvex functions, 25 which does not require additional variables and constraints. Readers can refer to Gatzke et al.26 for more discussions on the convex relaxation techniques. If Problem LBP is assumed to satisfy a constraint qualification (which implies strong duality) for any y ∈ Y for which Problem LBP is feasible, then it can be equivalently transformed into a master problem via the principles of projection and dualization. 15 The master problem contains an infinite number of constraints and is usually difficult to solve directly. Instead, it is solved via solving a sequence of Primal Bounding Problems (PBP), Feasibility Problems (FP), and Relaxed Master Problems (RMP), which are much easier to solve. The primal bounding problem is constructed at each iteration k by restricting the binary variables to specific values, say y = y(k), in the lower bounding problem, whose solution yields a valid upper bound for the lower bounding problem (and hence the master problem). Furthermore, it can be naturally decomposed into s subproblems of the following form:

for Problem P because of the usual loss of strong duality when nonconvex functions are participating. Recently a novel decomposition method, called nonconvex generalized Benders decomposition (NGBD),15 has been developed to solve Problem P to guaranteed global optimality, and case study results10,15 have demonstrated its dramatic computational advantages over the state-of-the-art branch-andreduce optimizer, BARON.16 In NGBD, a lower bounding problem, for which the dual transformation is rigorous, is generated by convex relaxation of the nonconvex functions in Problem P. The solution of this lower bounding problem not only yields a sequence of nondecreasing lower bounds, but also gives information to identify subproblems that generate a sequence of nonincreasing upper bounds. A global optimum of the problem is verified once the upper and lower bounds converge. The major advantage of NGBD is that the sizes of the subproblems to be solved are all independent of the number of scenarios addressed. In NGBD, the lower bounding problem serves as a surrogate for Problem P for the purpose of valid decomposition, and the tightness of the convex relaxation (i.e., the closeness of the lower bounding problem to Problem P) determines the quality of the information generated through the decomposition. Therefore, the tighter the convex relaxation is, the faster the NGBD may converge. However, NGBD does not reduce the convex relaxation gap during the solution procedure because it does not branch on the variables in the full search space, so it may have to visit most or even all of the binary variable realizations when the relaxation gap is large (although finite termination is still guaranteed). Recently, it has been recognized in the process systems engineering literature that piecewise linear relaxation enables much tighter relaxation of bilinear programs and can expedite global optimization significantly in the branch-and-bound framework17−20 (while the notion of piecewise linear relaxation dates back to the 1980s21). This paper generalizes the notion of piecewise relaxation to factorable functions in the context of McCormick relaxation,22,23 and integrates the so-obtained piecewise McCormick relaxation technique into the NGBD algorithm to reduce the gap between Problem P and the lower bounding problem. The remaining part of the paper is organized as follows: a brief introduction to NGBD is given in Section 2 and a piecewise McCormick relaxation approach for factorable programs is developed in Section 3. Then, Section 4 describes the additional subproblems to be solved when integrating piecewise convex relaxation into NGBD and the updated NGBD algorithm is given. The computational advantage of the NGBD with piecewise convex relaxation is demonstrated through a pump network configuration problem and an energy polygeneration synthesis problem in Section 5, and the paper ends with conclusions in Section 6.

objPBP (y(k)) = min wh(chTy(k) + uf , h(xh , eh)) h

s.t.

(PBP-hk) (k)

where objPBPh(y ) denotes the optimal objective value of Problem PBP-hk, h = 1, ..., s. Obviously, ∑sh=1objPBPh(y(k)) = objPBP(y(k)), where objPBP(y(k)) is the optimal objective value of the primal bounding problem for y = y(k). When the primal bounding problem is infeasible for y = y(k), a corresponding feasibility problem is solved, which also can be decomposed into s subproblems of the following form: objFP (y(k)) = h

OVERVIEW OF NONCONVEX GENERALIZED BENDERS DECOMPOSITION In NGBD, Problem P is first relaxed by convex relaxation of the nonconvex functions into the following lower bounding problem:

where objFPh(y(k)) denotes the optimal objective value of Problem FP-hk, h = 1, ..., s. Obviously, ∑sh=1objFPh(y(k)) = objFP(y(k)), where objFP(y(k)) is the optimal objective value of the feasibility problem. The relaxed master problem is constructed at each iteration k by relaxing the master problem with a finite number of constraints that are derived according to the solution information of all the

∀ h ∈ {1, ..., s},

∀ h ∈ {1, ... , s},

y∈Y

z h ∈ Zh (FP-hk)

∑ wh(chTy + uf ,h(xh , eh))

(xh , eh) ∈ Dh ,

ug , h(xh , eh) + Bh y(k) ≤ zh ,

(xh , eh) ∈ Dh ,

x1,..., xs , e1,..., es , y h = 1

ug , h(xh , eh) + Bh y ≤ 0,

min wh||zh||

xh , eh , zh

s.t.

s

s.t.

ug , h(xh , eh) + Bh y(k) ≤ 0,

(xh , eh) ∈ Dh

■

min

xh , eh

(LBP) 7288

dx.doi.org/10.1021/ie201262f | Ind. Eng. Chem. Res. 2012, 51, 7287−7299


Article

or an indication of the infeasibility of Problem P. The algorithm flowchart of NGBD is illustrated in comparison with the NGBD with piecewise convex relaxation, by Figure 1 in a subsequent section. More details of NGBD and its theoretical foundation can be found in Li et al.15

previously solved primal bounding and feasibility problems, as follows: min η η,y

s

■

η ≥ objPBP(y(j)) + (∑ (whchT + (λh(j))T Bh ))(y − y(j)),

s.t.

h=1 k

PIECEWISE MCCORMICK RELAXATION OF FACTORABLE PROGRAMS This section develops a piecewise McCormick relaxation framework for generating tighter lower bounding problems for factorable programs in NGBD. Note that this presentation follows the approach of explicitly introducing extra variables and constraints into Problem LBP for each factor, in contrast to the approach described in McCormick22 and Scott et al.23 Factorable programs are programs with factorable functions. A function f: S ⊂ nx →  is called factorable if there exist factors q1, ..., qL such that the function can be represented by a finite sequence of addition, multiplication, and univariate functions in the following form:

s

0 ≥ objFP(y(i)) + (∑ (μh(i))T Bh )(y − y(i)),

∀j∈T ,

h=1

∀ i ∈ Sk ,

∑ r ∈ {r : yr(t ) = 1, r = 1,..., ny}

yr −

× yr ≤ |{r: yr(t ) = 1}| − 1,

∑ r ∈ {r : yr(t ) = 0, r = 1,..., ny}

∀ t ∈ T k ∪ Sk ,

y ∈ Y, η ∈ 

(RMP-k)

where the index sets T k = {j ∈ {1, ..., k }: Problem LBP is feasible for y = y(j)} S k = {i ∈ {1, ..., k }: Problem LBP is infeasible for y = y(i)} (k) λ(j) = h are the Lagrange multipliers for Problem PBP-hk (for y (j) k (j) y ), which form an optimality cut for iteration j (∀j ∈ T ). μh are the Lagrange multipliers for Problem FP-hk (for y(k) = y(i)), which form a feasibility cut for iteration i (∀i ∈ Sk). The last group of constraints represent a set of canonical integer cuts that prevent the previously examined integer realizations from becoming a solution.27 The solution of the relaxed master problem yields a valid lower bound for the master problem (and therefore Problem P) augmented with the integer cuts. Notice that the relaxed master problem is a MILP whose size is determined by the current iteration number but is independent of the number of scenarios. In case Tk = ⌀, the relaxed master problem is unbounded, then the following feasibility version of it is solved to make the algorithm proceed: ny

(1)

ql = xl ,

l ∈ {1 , ..., nx}

(2)

ql = qi + qj ,

(3)

ql = qiqj ,

(4)

ql = Ul(qi),

(5)

f (x) = qL(x)

∀ (l , i , j) ∈ ΩA

∀ (l , i , j) ∈ ΩM

∀ (l , i) ∈ ΩU

min ∑ yi y

i=1 s

s.t.

0 ≥ objFP (y(i)) + (∑ (μh(i))T Bh )(y − y(i)), h=1

∀ i ∈ Sk ,

∑ r ∈ {r : yr(t ) = 1, r = 1,..., ny}

≤ |{r: yr(t ) = 1}| − 1,

yr −

∑ r ∈ {r : yr(t ) = 0, r = 1,..., ny}

∀ t ∈ Sk ,

y ∈ Y,

Here, ΩA, ΩM, and ΩU are index sets that contain indices of the factors associated with additions, multiplications, and univariate functions, respectively. If upper and lower bounds on each factor are given or inferred from the known bounds on other factors (say, via interval analysis28), a valid convex underestimator and concav overestimator of f can be generated by performing McCormick relaxation for each multiplication and nonaffine univariate function with associated variable bounds. Relaxation of Multiplications. For multiplication, say z = xy, with known upper and lower bounds on x and y, say, xlo, xup, ylo, yup, the McCormick relaxation can be written as

yr

η∈

(FRMP-k)

On the other hand, a restriction of Problem P, which is called the Primal Problem (PP), is constructed at iteration l by restricting y = y(l) in Problem P, whose optimal objective value is a valid upper bound for Problem P. The primal problem can be further decomposed into s subproblems of the following form: objPP (y(l)) = min wh(chTy(l) + fh (xh)) h

xh

s.t.

gh(xh) + Bh y(l) ≤ 0,

z r ≥ x upy + xy up − x upy up , xh ∈ Xh

z r ≥ x loy + xy lo − x upy lo ,

(PP-hl)

z r ≤ x upy + xy lo − x upy lo ,

where objPPh(y(l)) denotes the optimal objective value of Problem PP-hl, h = 1, ..., s. Obviously, ∑sh=1objPPh(y(l)) = objPP(y(l)), where objPP(y(l)) is the optimal objective value of the primal problem. The nonconvex nonlinear programming (NLP) problem (PP-hl) can be solved to ε-optimality in finite time by state-of-the-art branch-andbound global optimization solvers, such as BARON, and the solution can be significantly accelerated by adding an additional cut derived from the solution of the previously solved subproblems.15 The aforementioned subproblems are solved iteratively in NGBD until a termination with global optimum of Problem P

z r ≤ x loy + xy up − x loy up , x lo ≤ x ≤ x up , y lo ≤ y ≤ y up

(1)

r

where z denotes the value of the relaxed multiplication. Define the relaxation gap of a function over [xlo,xup] × [ylo,yup] to be the maximum difference of the function and its relaxation over this domain, then the relaxation gap of multiplication diminishes to zero 7289



Article

as either |xup − xlo| or |yup − ylo| approaches zero, because according to eq 1 (and z = xy):

divided into K pieces. Notice that only the domain of one variable in each multiplication needs to be partitioned, according to eq 3. Relaxation of Univariate Functions. For a (nonaffine) univariate function, say z = U(x), defined on [xlo,xup], its relaxation zr can be expressed as:

z r − z ≥ x upy + xy up − x upy up − xy = (x up − x)(y − y up ) ≥ (x up − x lo)(y lo − y up ), z r − z ≥ x loy + xy lo − x loy lo − xy = (x lo − x)(y − y lo ) ≥ (x lo − x up)(y up − y lo ),

z r ≤ U conc(x , x lo , x up), z r ≥ U conv(x , x lo , x up),

z r − z ≤ x upy + xy lo − x upy lo − xy = (x up − x)(y − y lo ) ≤ (x up − x up)(y up − y lo ), z r − z ≤ x loy + xy up − x loy up − xy = (x lo − x)(y − y up ) ≤ (x lo − x up)(y lo − y up )

so

x lo ≤ x ≤ x up

where Uconc(x,xlo,xup) denotes a concave overestimator of U(x) over [xlo,xup] and U(x,xlo,xup) denotes a convex underestimator of U(x) over [xlo,xup], and it is assumed that a constant C exists such that (2)

|U conc(x , x lo , x up) − U conv(x , x lo , x up)| ≤ C|x up − x lo| ,

|z r − z| ≤ |x up − x lo||y up − y lo |

sup x lo ≤ x ≤ x up , y lo ≤ y ≤ y up

∀ x ∈ [x lo , x up]

(3)

Therefore, the relaxation gap can be decreased by partitioning the x domain, [xlo,xup] into K subdomains, i.e., picking K + 1 points x1 < x2 < ··· < xK+1 such that x1 = xlo, xK+1 = xup and performing McCormick relaxation on each individual subdomain. Furthermore, each subdomain is assigned a binary variable δk to determine if the McCormick relaxation on this subdomain is used or not, and ∑Kk=1δk = 1 is enforced. Thus the McCormick relaxation of the multiplication can be upgraded into: K

zr =

K

sup |z − z r| ≤ lo

x ≤x≤x

≥x

yk + xky

up

k + 1 up

y δk ,

z kr ≥ x kyk + xky lo − x ky lo δk , zkr

≤x

k+1

lo

yk + xky − x

k+1

lo

up

δ kx ≤ x k ≤ δ kx

x ≤ x ≤ x up

K

zr =

k + 1 lo

y δk ,

(7)

K

∑ zkr , x = ∑ xk , k=1

zkr ≤ x kyk + xky up − x ky up δk , k

sup |U conc(x , x lo , x up) − U conv(x , x lo , x up)| lo

Therefore, the relaxation gap for the univariate function diminishes to zero as |xup − xlo| approaches zero. So, the gap can be decreased by partitioning the x domain as before and performing McCormick relaxation for each individual subdomain, as follows:

k=1

−x

up

≤ C|x up − x lo|

K

k=1

k+1

(6)

From eqs 5 and 6,

∑ zkr , x = ∑ xk , y = ∑ yk , k=1

z kr

(5)

zkr

,

≤U

k=1

conc, k

(xk , δk , x k , x k + 1),

zkr ≥ U conv, k(xk , δk , x k , x k + 1)

δky ≤ yk ≤ δky ,

K

δ kx k ≤ x k ≤ δ kx k + 1 ,

K

∑ δk = 1,

∑ δk = 1, k=1

k = 1, ..., K (8)

k=1

k = 1, ..., K

where the binary variables and partitioning points are as defined before, Uconc,k(x,1,xk,xk+1) and Uconv,k(x,1,xk,xk+1) denote a concave relaxation and a convex relaxation of U(x) on [xk,xk+1 ], respectively, and the following condition,

(4)

Proposition 1. If zr is feasible for the constraints in eq 4, then it is feasible for the constraints in eq 1 as well. Proof. See Appendix for the proof. The piecewise McCormick relaxation eq 4 characterizes a disjunctive polyhedral set that contains the value of the multiplication over its domain, and its continuous relaxation leads to the convex hull of the disjunctive set.29 So this formulation is favorable for mixed-integer programming. Let nM be the total number of multiplications and nPM denotes the total number of variables whose domains need to be partitioned for the multiplications, then upgrading relaxation eq 1 into eq 4 incurs KnPM additional binary variables, 3(K − 1)nM continuous variables, and (8K − 1)nM linear constraints provided all the partitioned domains are

U conc, k(0, 0, x k , x k + 1) = 0, U conv, k(0, 0, x k , x k + 1) = 0 k = 1, ..., K

(9)

is enforced such that zrk = 0 if δk = 0. Notice that eq 9 may not hold for an arbitrary convex or concave relaxation. In this case, Uconc,k or Uconv,k can be further relaxed into 7290



Article

where δh are nδ binary variables introduced for the piecewise McCormick relaxation, ẽh involves eh in Problem LBP as well as additional continuous variables for the piecewise relaxation, and Dpw h is defined by constraints for the piecewise relaxation. Let nP be the total number of variables whose domains need to be partitioned, then Problem LBP-PCR has sKnP additional binary variables, s[(3K − 1)nM + KnU] continuous variables, and s[(8K − 1)nM + 2(K − 1)nU] linear constraints compared to Problem LBP, provided all the partitioned domains are divided into K pieces. Problem LBP-PCR is a convex MINLP or a MILP for which commercial solvers are available, such as DICOPT31 (for convex MINLP) and CPLEX32 (for MILP). If Problem LBP-PCR is a convex MINLP that is expensive to solve, it can be further relaxed into a MILP (which is usually much easier to solve) via outer linearization.24

affine functions through outer linearization, 24 and then eq 8 becomes: K

zr =

K

∑ zkr , x = ∑ xk , k=1

i = 1, ..., K conc, k

zkr ≥ αjconv, kxk + βjconv, k δk ,

j = 1, ..., K conv, k

≤

αiconc, kxk

k=1

+ βiconc, k δk ,

zkr

K

δ kx k ≤ x k ≤ δ kx k + 1 ,

∑ δk = 1,

k = 1, ..., K

k=1

(10)

■

where Kconc,k and Kconv,k denote the numbers of linearization points for outer linearization of Uconc,k and Uconv,k to construct the affine over- and under-estimators, respectively. Readers can refer to Tawarmalani and Sahinidis24 for the way to select linearization points. Now the new concave and convex relaxations are α conc,k xk + βconc,k δk and α conv,k xk + i i j conv,k βj δk, respectively, and they obviously satisfy eq 9.

NONCONVEX GENERALIZED BENDERS DECOMPOSITION WITH PIECEWISE CONVEX RELAXATION The new lower bounding problem (Problem LBP-PCR) obtained by the proposed piecewise McCormick relaxation technique can be exploited to improve the NGBD method. First, new subproblems that offer improved information for NGBD are developed by manipulating Problem LBP-PCR. Then the updated NGBD algorithm integrating the new subproblems is described, which still retains the convergence property. New Subproblems. Fixing y = y(k) in Problem LBP-PCR leads to a convex MINLP or MILP, which can be further decomposed into s subproblems in the following form:

Proposition 2. Suppose for any xk ∈ [xk,xk+1], U conc , k(xk , δk , x k , x k + 1) ≤ U conc(xk , x lo , x up), U conv , k(xk , δk , x k , x k + 1) ≥ U conv(xk , x lo , x up) k = 1, ..., K

(11)

r

Then if z is feasible for the constraints in eq 8, it is also feasible for the constraints in eq 5. Proof. See Appendix for the proof. The piecewise McCormick relaxation eq 8 characterizes a disjunctive convex set that contains the univariate function value over its domain, and this formulation covers the generalized disjunctive programming formulation30 whose continuous relaxation leads to the convex hull of the disjunctive set. So again this piecewise relaxation scheme benefits the solution of the resulting mixed-integer program. Let nU be the total number of univariate functions and nPU be the number of variables whose domains need to be partitioned for the univariate functions, then upgrading relaxation eq 5 into eq 8 incurs KnPU additional binary variables, KnU continuous variables and 2(K − 1)nU linear constraints provided all the partitioned domains are divided into K pieces. Lower Bounding Problem with Piecewise McCormick Relaxation. With the proposed piecewise McCormick relaxation formulations, a factorable program can be relaxed by relaxing the multiplications and univariate functions in its factorable representation. Then a new lower bounding problem, which is a tighter (or equal) relaxation compared to Problem LBP, according to Proposition 1 and Proposition 2, can be developed as follows:

s

objPBP − PCR (y(k)) = h

where objPBP−PCRh(y ) denotes the optimal objective value of Problem PBP-PCR-hk, h = 1, ..., s, and ∑sh=1objPBP−PCRh(y(k)) = objPBP−PCR(y(k)), where objPBP−PCR(y(k)) denotes the optimal objective value of Problem LBP-PCR for y = y(k). With the piecewise convex relaxation, Problem PBP-PCR-hk is a tighter relaxation of Problem FP-hk compared to Problem PBP-hk, so objPP (y(k)) ≥ objPBP − PCR (y(k)) ≥ objPBP (y(k)) h

h

h

(12)

Summing eq 12 over h = 1, ..., s yields objPP(y(k)) ≥ objPBP − PCR (y(k)) ≥ objPBP(y(k))

(13)

Therefore, Problem PBP-PCR-hk leads to a better estimate of the optimal objective value of Problem P for y = y(k) compared to Problem PBP-hk. Problem PBP-PCR-hk can also be “dualized” into the following problem:

s

(xh , eh̃ , δh) ∈ Dhpw × {0, 1}nδ , y∈Y

ugpw, h(xh , eh̃ , δh) + Bh y(k) ≤ 0, (PBP-PCR-hk)

∑ wh(chTy + upw f , h(xh , eh̃ , δh))

ugpw, h(xh , eh̃ , δh) + Bh y ≤ 0,

h=1

(k)

δ1,..., δs , y

s.t.

∑ wh(chTy(k) + upw f , h(xh , eh̃ , δh))

(xh , eh̃ , δh) ∈ Dhpw

x1,..., xs , h=1 e∼ ,..., e∼ , 1

xh , eh , δh

s.t.

s

min

min ∼

(k) objDPBP − PCR (y(k) , λh̃ ) =

∀ h ∈ {1, ..., s},

h

min wh(chTy(k) + upw f , h(xh , eh̃ , δh))

xh , e∼h , δh

(k)

+ (λh̃ )T (ugpw, h(xh , eh̃ , δh) + Bh y(k))

∀ h ∈ {1, ..., s},

s.t.

(LBP-PCR)

(xh , eh̃ , δh) ∈ Dhpw

(DPBP-PCR-hk) 7291



Article

where λ̃(k) h denotes the KKT multipliers obtained at the solution of Problem PBP-PCR-hk, and h = 1, ..., s. Let ob(̃ k) ̃ jDPBP−PCR(y(k),λ(̃ k)) = ∑sh=1objDPBP−PCRh(y(k),λ(k) = h ), where λ (k) (k) ̃ ̃ (λl , ..., λs ), then construct the following enhanced relaxed master problem

NGBD-PCR Algorithm. Initialize: 1 Iteration counters k = 0, l = 1, and the index sets T 0 = 0, T 0pw = 0, T̃ 0 = 0, S 0= 0, U 0= 0. 2 Upper bound on the problem UBD = +∞, upper bound on the lower bounding problem UBDPB = +∞, lower bound on the lower bounding problem LBD = −∞. 3 Set tolerances εh and ε such that ∑sh=1εh ≤ ε. 4 Integer realization y(1) is given. repeat if k = 0 or (Problem ERMP-k/FRMP-k is feasible and LBD < UBDPB and LBD < UBD − ε) then repeat Set k = k + 1 1) Solve the decomposed primal bounding subproblem PBP-hk for each scenario h = 1, ..., s sequentially. If Problem PBP-hk is feasible with Lagrange multipliers λ(k) h for all the scenarios, add an optimality cut to Problem k k−1 ERMP-k according to λl(k), ..., λ(k) ∪ {k}. s , set T = T 2) If Problem PBP-hk is feasible for all the scenarios, solve subproblem PBP-PCR-hk for each scenario h = 1, ..., s sequentially. If Problem PBP-PCR-hk is feasible with ̃ for all the scenarios, set Tkpw = KKT multipliers λ(k) h Tk−1 ∪ {k}. In this case, if objPBP−PCR(y(k)) < UBDPB, pw update UBDPB = objPBP−PCR(y(k)), y* = y(k), k* = k, solve ̃ subproblem DPBP-PCR-hk with λ(k) h for each scenario h = 1, ..., s sequentially, add an optimality cut to Problem ERMP-k, and set T̃ k = T̃ k−1∪{k}. 3) If Problem PBP-hk is infeasible for one scenario, say ĥ, stop solving it for the remaining scenarios and set Sk = ̂ Sk−1 ∪ {k}. Then, set μ(k) h = 0 for h = 1, ..., h − 1, and solve the decomposed feasibility subproblem FP-hk for h = ĥ, ..., s and obtain the corresponding Lagrange multipliers μ(k) h . Add a feasibility cut to Problem ERMP(k) k/FRMP-k according to μ(k) 1 , ..., μs . k 4) If T = 0, solve Problem FRMP-k; otherwise, solve Problem ERMP-k. In the latter case, set LBD to the optimal objective value of Problem ERMP-k if Problem ERMP-k is feasible. In both cases, set y(k+1) to the y value at the solution of either problem. until LBD ≥ UBDPB or (Problem ERMP-k/FRMP-k is infeasible). end if if UBDPB < UBD − ε then 1 Solve the decomposed primal subproblem PP-h* (i.e., for y = y*) to εh-optimality for each scenario h = 1, ..., s sequentially. Set Ul = Ul−1{k*}. If Problem PP-h* is feasible with optimum x*h for all the scenarios and objPP(y*) < UBD, update UBD = objPP(y*) and set y*p = y*, x*p,h = x*h for h = 1, ..., s. 2 If Tkpw\Ul = ⌀, set UBDPB = +∞.

min η η,y

s

s.t.

(j) (j) η ≥ objDPBP − PCR (y(j) , λ ̃ ) + (∑ (whchT + (λh̃ )T Bh))(y − y(j)), h=1 k

∀ j ∈ T̃ , s

η ≥ objPBP(y(j)) + (∑ (whchT + (λh(j))T Bh))(y − y(j)), h=1

∀ j ∈ Tk, s

0 ≥ objFP(y(i)) + (∑ (μh(i))T Bh)(y − y(i)), h=1 k

∀i∈S ,

∑ r ∈ {r : yr(t ) = 1, r = 1,..., ny}

yr −

∑ r ∈ {r : yr(t ) = 0, r = 1,..., ny}

yr ≤ |{r: yr(t ) = 1}| − 1,

∀ t ∈ T k ∪ Sk , y ∈ Y, η ∈ 

(ERMP-k)

where T̃ k ⊂ Tkpw = {j: Problem LBP-PCR is feasible for y = y(j)} since j ∈ Tkpw is necessary for constructing Problem DPBP-PCR-hk. Also notice that T̃ kpw ⊂ Tk due to the tighter relaxation, so T̃ k ⊂ Tk as well. There are different ways to determine the set T̃ k, i.e., to determine the integer realizations for which subproblems DPBP-PCR-hk are solved to derive additional optimality cuts. Since Problem DPBP-PCR-hk is usually expensive to solve (compared to other subproblems), it is better to keep the size of set T̃ k small. In this paper, a heuristic is adopted which solves subproblems DPBP-PCR-hk only when the corresponding subproblems (Problem PBP-PCR-hk) updates the current upper bound for Problem LBP-PCR. (See the detailed algorithm in the next subsection.) The enhanced relaxed master Problem ERMP-k differs from Problem RMP-k in that it has a group of additional optimality cuts from the solution of the new subproblems PBP-PCR-hk and DPBPPCR-hk. Proposition 3. Compared to Problem RMP-k, Problem ERMPk is a tighter (or equal) relaxation of Problem P when Problem P is augmented with the relevant integer cuts excluding the previously examined integer realizations. Proof. See Appendix for the proof. If Tk = 0, Problem ERMP-k is unbounded, so the feasibility relaxed master Problem FRMP-k introduced in Section 2 is solved instead. In principle, the feasibility subproblem FP-hk can also be upgraded with piecewise convex relaxations to generate additional valid feasibility cuts for the relaxed master problem, but our computational experience indicates that these cuts do not significantly accelerate the convergence, which may be because the integer cuts in the problem (which prevent visiting an integer realization twice) are already strong enough. So these feasibility cuts are not addressed in this paper. For ease of discussion, the NGBD incorporating the piecewise convex relaxations as described in this subsection is called NGBD-PCR in this paper.

3 If Tkpw\Ul ≠ ⌀, pick i ∈ arg minj∈Tkpw \Ul{objPBP−PCR(y(j))}. Update UBDPB = objPBP‑PCR(y(i)), y* = y(i), k* = i. Set l = l + 1. end if until UBDPB ≥ UBD − ε and (Problem ERMP-k/FRMP-k is infeasible or LBD ≥ UBD − ε). An ε-optimal solution of Problem P is given by (y*p,x*p,1,...,x*p,s) or Problem P is infeasible. 7292



Article

Figure 1. Flowchart of the NGBD-PCR algorithm.

Discussion. Figure 1 illustrates the NGBD-PCR algorithm and the two major differences from the NGBD algorithm are highlighted in gray, i.e., solving the new subproblems PBPPCR-hk and DPBP-PCR-hk and solving the enhanced relaxed master Problem ERMP-k instead of Problem RMP-k. Another important change (which is not shown in the figure) is that the selection of the integer realization y(k) for constructing primal subproblems (PP-hl) is based on objPBP−PCR(y(k)) instead of objPBP(y(k)). According to eq 13, the new selection criterion is more likely to locate a global optimum earlier. Theorem 1. If all the subproblems can be solved to ε-optimality in a finite number of steps, then the NGBD-PCR algorithm terminates in a finite number of steps with an ε-optimal solution of Problem P or an indication that Problem P is infeasible. Proof. The solution procedure in NGBD-PCR is the same as that in NGBD except that NGBD-PCR solves a finite number of additional subproblems (since set T̃ k ⊂ Y is finite) and it solves a different relaxed master problem. Since all these subproblems can be solved finitely, and the new relaxed master problem still prevents visiting the same integer realization twice, the convergence property holds for NGBD-PCR as well according to the proof of NGBD convergence in Li et al.15

■

GAMS, their CPU times are simply the times reported by GAMS. Two CPU times were recorded for NGBD and NGBD-PCR in all the case studies. One is called “total solver time”, which means the sum of the times for the NLP and LP/MILP solvers to solve the subproblems in NGBD/NGBD-PCR. The other is called “GAMS time”, which means the overall GAMS running time for NGBD/NGBD-PCR to locate an ε-optimal solution, including the “total solver time” and the time for auxiliary operations in GAMS (e.g., exchanging data between subproblems and compiling new subproblems). Notice that the “total solver time” only slightly underestimates the total solution time required by NGBD/NGBD-PCR, because it primarily ignores the overhead for exchanging the data among the subproblems, which is a small fraction of the “total solver time” if the algorithm is implemented in a general purpose programming language (such as C and Fortran). The “GAMS time” significantly overestimates the total solution time, because GAMS has to recompile the subproblems frequently when realizing NGBD/ NGBD-PCR and the total compilation time may dominate the overall GAMS running time. Therefore, the “total solver time” may be deemed a more precise estimate of the NGBD/NGBDPCR solution time and the “GAMS time” may be deemed a (loose) upper bound on the solution time. Since this paper does not focus on scenario generation for stochastic programs, all the uncertain parameters in the case studies are assumed to obey normal distributions and are independent of each other. A naive sampling rule was used to generate scenarios for normally distributed uncertain parameters, which is detailed in the Appendix. Case 1: Pump Network Configuration. This is a stochastic MINLP problem for the optimal configuration of a centrifugal pump network. The objective of the optimization is to minimize the expected annualized cost while achieving a prespecified pressure rise with a given total flow rate. The superstructure of the pump network is shown in Figure 2a. The deterministic version of the problem was initially presented in Westerlund et al.34 and then updated in Adjiman et al.6 with a set of additional linear constraints for tighter relaxation in global optimization. The problem was further reformulated to reduce the number of

CASE STUDIES

The purpose of the case studies is to compare the performance of NGBD-PCR with NGBD and commercial global and local MINLP solvers for problems in the form of Problem P. BARON was used as the global MINLP solver and DICOPT the local solver for the comparison. The case study problems were solved on a computer allocated a single 2.66 GHz CPU, 512 MB memory and running Linux Kernel. GAMS 23.5.233 was used to formulate the problems, program the NGBD and NGBD-PCR algorithms, and interface the solvers for the subproblems. NGBD and NGBD-PCR employed CPLEX 12.2.0 to solve LP/MILP subproblems and BARON 9.0.6 to solve the primal subproblem PP-hl to global optimality. BARON 9.0.6 employed CPLEX 12.2.0 as the linear program (LP) solver and CONOPT 3.14 as the local NLP solver. As DICOPT and BARON are commercial solvers available in 7293



Article

nonlinear functions while exhibiting the structure of Problem P in Li et al.15 (where only the two levels of the network were

addressed), and this formulation is adopted in this paper. Therefore the target stochastic MINLP model is as follows:

⎡ 2 ⎛ ⎞⎤ 2 2 2 2 j+k−2 j+k−2 y fix ⎢ ⎜ ∑ wh⎢∑ ⎜Ci ∑ ∑ 2 yj ,k ,i + Ci ∑ ∑ 2 Pj , k , i , h⎟⎟⎥⎥ ⎠⎦ h=1 j=1 k=1 j=1 k=1 ⎣ i=1 ⎝ s

min

yjnp , y ns , z , y , P , Δpi , h , vi̇ , h , xi , h , ωi , hPjy, k , i , h , Δpky, i , h , vjẏ , i , h ,i k ,i i j,k ,i i ,h

s.t.

⎛ ωi , h ⎞3 ⎛ ωi , h ⎞2 ⎛ ωi , h ⎞ 2 Pi , h − αi⎜ ⎟ − βi ⎜ ⎟ vi̇ , h − γi⎜ ⎟vi̇ , h = 0, ⎝ ωmax ⎠ ⎝ ωmax ⎠ ⎝ ωmax ⎠

⎛ ωi , h ⎞2 ⎛ ωi , h ⎞ ⎛ ωi , h ⎞ 2 Δpi , h − ai , h⎜ ⎟ − bi⎜ ⎟vi̇ , h − ci⎜ ⎟vi̇ , h = 0, ⎝ ωmax ⎠ ⎝ ωmax ⎠ ⎝ ωmax ⎠ 3

∑ xi , h = 1, i=1 2

∑ 2j − 1vj̇y, i , h − xi , hVtot = 0, j=1 2

ΔPtotzi −

∑ 2k− 1Δpky, i , h

= 0,

k=1

0 ≤ xi , h ≤ 1, 0 ≤ Pi , h ≤

0 ≤ vi̇ , h ≤ Vtot ,

Pimax ,

Δpi , h ≤ ziΔPtot ,

0 ≤ ωi , h ≤ ωmax , Pi , h ≤ ziPimax ,

0 ≤ Δpi , h ≤ ΔPtot , vi̇ , h ≤ ziVtot ,

xi , h ≤ zi ,

2

ωi , h ≤ ziωmax , 2

2

∑ 2j − 1yjnp, i ≤ 3zi ,

∑ 2k− 1yknp, i ≤ 3zi ,

j=1

k=1

2

2

∑ zi ≥ 1,

∑ 2j − 1yjnp, i ≥ zi ,

∑ 2k− 1ykns, i ≥ zi ,

i=1

j=1

k=1

yj , k , i ≤

yjnp , ,i

yj , k , i ≤

ykns, i ,

0 ≤ Pjy, k , i , h ≤ Pimaxyj , k , i ,

∀ i ∈ {1, 2, 3},

)

)

vi̇ , h − Vtot 1 − yjnp ≤ vjẏ , i , h ≤ vi̇ , h , ,i

0 ≤ Δpky, i , h ≤ ΔPtotykns, i ,

{ }

(

Pi , h − Pimax 1 − yj , k , i ≤ Pjy, k , i , h ≤ Pi , h ,

(

0 ≤ vjẏ , i , h ≤ Vtotyjnp , ,i

yjnp ∈ 0, 1 , ,i

yj , k , i ≥ yjnp + ykns, i − 1, ,i

(

)

Δpi , h − ΔPtot 1 − ykns, i ≤ Δpky, i , h ≤ Δpi , h ,

{ }

ykns, i ∈ 0, 1 , ∀ j ∈ {1, 2},

{ }

zi ∈ 0, 1 ,

{ }

yj , k , i ∈ 0, 1 ,

∀ k ∈ {1, 2},

∀ h ∈ {1, ..., s}

pressure rise with rotation speed and rate of flow going through the pump. This model was obtained from regression,34 so the nominal values of the regression parameters αi, βi, γi, ai, bi, ci may not reflect the performance of the real system. Therefore, in this case study, parameters a1, a2, and a3 are assumed to be uncertain and independently obey normal distributions with means of 629.0, 215.0, 361.0 and standard deviations of 62.90, 21.50, 36.10, respectively. The case study addresses 5 realizations for each uncertain parameter so the problem includes 125 scenarios in total. The probability of scenario h, wh, and the associated values of uncertain parameters are calculated according to the sampling rule explained in the Appendix.

Here, index i indicates different levels of the pump network, and index h indicates different scenarios. Variable zj decides whether pump level i is developed or not, ynp j,i decides the number of parallel ns lines on level i (which is ∑2j=12j−1ynp j,i ), and yk,i decides the number of pumps on each pump line on level i (which is ∑2k=12k−1yns k,i). The continuous variables include the fraction of total flow going to level i, xi, the flow rate on each line at level i, vi, the rotation speed of pumps on level i, Ωi, the power requirements at level i, Pi, and the pressure rise at level i, Δpi. More details of the problem can be found in Li et al.15 The first two groups of equations in the problem denote the pump performance model, relating power consumption and 7294



Article

Figure 2. The superstructure and design results of the pump network configuration problem.

Details for other parameters can be found in Westerlund et al.34 and Adjiman et al.6 The resulting problem contains 27 binary variables (12 of them are dummy variables for reformulation so there are only 215 elements in the set Y) and 7125 continuous variables. The convex and concave envelopes of the quadratic and cubic functions in the model include nonlinear functions, so the resulting Problem PBP-PCR-hk in NGBD-PCR is a convex MINLP, which is much more expensive to solve compared to other subproblems. So each of these nonlinear functions was linearized with 5 linearization points to yield MILP subproblems and the details of the outer linearization are given in the Appendix. The relative and absolute termination criteria for optimization were set to be 10−3 and the initial value of all the integers was 0. First, the case study compares the designs with the deterministic model (in which s = 1 and the uncertain parameters are represented with their mean values) and the stochastic model. The designed structures for the pump work with the two models are shown in Figure 2, parts b and c, respectively. The nominal optimal objective values of the designs with the deterministic and the stochastic models are 128.89 FIM and 131.12 FIM, respectively. However, the design with the deterministic model is infeasible for some realizations of the uncertain parameters, while the design with the stochastic model is feasible for all the 125 scenarios with an expected optimal objective value of 132.33 FIM over these scenarios. Second, the case study compares the performance of DICOPT, BARON, NGBD, and NGBD-PCR for the stochastic MINLP. The domains of the variables to be partitioned in NGBD-PCR are assumed to be partitioned uniformly into the same number of subdomains, denoted by K, and NGBD-PCR was implemented for three different K values to show the effect of K on the performance of NGBD-PCR. Table 1 summaries the results. It can be seen that DICOPT failed to locate a feasible solution with the given initial point, although it only took 2.8 s to terminate. BARON did not return an ε-optimal solution within 106 s, which is about 10 times the GAMS time for NGBD to locate an ε-optimal solution. This indicates that NGBD can reduce the solution time (with BARON) by an order of magnitude for the pump network configuration problem, because the GAMS time is a loose upper bound of the solution time with NGBD as explained at the beginning of this section. Table 1 also shows that most of the NGBD total solver time was spent on solving the most difficult subproblem PP-hl (which are nonconvex NLPs to be solved to global optimality). Therefore, the three NGBD-PCR methods further reduced the

solution time by solving much fewer subproblem PP-hl. As the piecewise relaxation helped to generate improved lower bounding problems that have a better chance to locate an optimal integer realization earlier, the NGBD-PCR method with more finely partitioned subdomains led to fewer subproblem PP-hl to be solved. On the other hand, integrating piecewise relaxation requires solving an additional MILP subproblem PBP-PCR-hk (and sometimes DPBP-PCR-hk as well) for each integer realization visited, and these MILPs are expensive to solve compared to subproblem PBP-hk in NGBD (although they are much easier than subproblem PP-hl). As K increases, these additional MILPs contain more integer variables and are more difficult to solve. Therefore, the solution time by NGBD-PCR depends on both the solution time for subproblem PP-hl and the solution time for the additional MILPs, and these two times are in principle negatively correlated. In this case study, NGBDPCR with K = 10 had the fastest solution because it achieved the best trade-off between the two times. This result suggests future research on a systematic approach to improve the domain partitioning for NGBD-PCR. Case 2: Energy Polygeneration Synthesis. This problem is the synthesis of an energy polygeneration plant coproducing power, liquid fuels (naphtha and diesel), and chemicals (methanol) from coal and biomass as feedstocks. The goal of the optimization is to determine the equipment capacities in the plant to achieve the best net present value over the plant lifetime while satisfying the design and operational constraints. The original version of the problem considered a static energy polygeneration plant where the operation is assumed to be fixed35 and the resulting MINLP only has one period/scenario. The nonlinear and nonconvex functions in the problem are the bilinear functions indispensable for tracking the properties of the streams in the system. Later on a multiperiod MINLP problem was proposed that accounts for the flexible operation of the plant in response to market conditions in different time periods, and the case study results showed that allowing the flexible operation could help to increase the net present value by 15−60%.36 In this case study, the multiperiod MINLP problem is revised to exhibit the structure of Problem P with an assumption of middle oil price and carbon tax and 100% operational flexibility. Two problem instances are considered here. One addresses 8 periods/scenarios which represent peak and off-peak times in 4 seasons, respectively, and it contains 70 binary variables and 4904 continuous variables. The other addresses 24 periods/ scenarios which represent peak and off-times in 12 months, respectively, and it contains 70 binary variables and 14712 7295



Article

Table 1. Results for the Pump Network Configuration Problem NGBD-PCR optimal objective GAMS time (s) total solver time (s) time for subproblems PBP-hk-* (s)b FP-hk (s) *RMP-k (s)c PP-hl (s) integers visitedd

DICOPT

BARON

NGBD

K=5

K = 10

K = 15

Infeasible 2.8

− −

132.3 115 147.6 5404.1

132.3 133 468.2 3936.5

132.3 116 698.0 4137.3

132.3 73 303.0 5037.7

6.8 9.7 693.5 4694.1 868/235

1224.5 4.5 666.3 2,041.2 822/149

3331.0 9.6 674.4 122.4 803/15

4553.4 10.8 433.5 39.9 632/4

a

a No global solution was returned within 1 000 000 s. bPBP-hk-* stands for subproblems PBP-hk, PBP-PCR-hk, and DPBP-PCR-hk. c*RMPk stands for subproblems RMP-k, FRMP-k, and ERMP-k. dTotal number of integer realizations/number of integer realizations visited by Problem PP-hl.

Table 2. Results for the Energy Polygeneration Synthesis Problem: 8 Scenarios NGBD-PCR optimal objective GAMS time (s) total solver time (s) time for subproblems PBP-hk-* (s)b FP-hk (s) *RMP-k (s)c PP-hl (s) integers visitedd

DICOPT

BARON

NGBD

K=5

K = 10

K = 15

778.7 7.5

−a −

−1 123.0 68 705.9 65 914.9

−1 123.0 23 734.5 21 493.9

−1 123.0 11 007.8 9 014.0

−1,123.0 7 421.0 6 103.7

69.9 2.0 329.1 65 514.0 464/396

1 290.6 3.5 452.3 19 747.6 471/209

2 992.8 5.6 323.8 5 691.8 477/93

4 547.6 3.6 215.6 1 337.0 382/28

No global solution was returned within 1 000 000 s. bPBP-hk-* stands for subproblems PBP-hk, PBP-PCR-hk, and DPBP-PCR-hk. c*RMPk stands for subproblems RMP-k, FRMP-k, and ERMP-k. dTotal number of integer realizations/number of integer realizations visited by PP-hl. a

Table 3. Results for the Energy Polygeneration Synthesis Problem: 24 Scenarios NGBD-PCR optimal objective GAMS time (s) total solver time (s) time for subproblems PBP-hk-* (s)b FP-hk (s) *RMP-k (s)c PP-hl (s) integers visitedd

DICOPT

BARON

NGBD

K=5

K = 10

K = 15

778.3 62.6

− −

−1 124.4 228 680.7 217 855.0

−1 124.4 64 722.6 54 818.9

−1 124.4 41 075.8 29 222.0

−1 124.4 30 809.7 23 524.8

289.5 6.4 636.5 216 922.7 613/537

3 632.5 1.1 677.2 50 508.1 590/234

9 130.8 4.0 584.5 19 502.8 631/113

20 655.6 11.3 504.5 2 353.5 537/27

a

a No global solution was returned within 30 days. bPBP-hk-* stands for subproblems PBP-hk, PBP-PCR-hk, and DPBP-PCR-hk. c*RMP-k stands for subproblems RMP-k, FRMP-k, and ERMP-k. dTotal number of integer realizations/number of integer realizations visited by PP-hl.

problems, respectively, while NGBD obtained an ε-optimal solution within 105 s and 3 days for the two problems, respectively. This indicates that NGBD can reduce the solution time (with BARON) by an order of magnitude for both problems. Since the solution time for subproblem PP-hl dominates the overall solution time in NGBD, the NGDB-PCR methods reduced the solution time by up to about another order of magnitude via solving fewer subproblem PP-hl. In turn, they solved additional MILP subproblems which offered information to reduce the number of PP-hl to be solved, and these MILPs are much more expensive than the corresponding LP subproblems PBP-hk in NGBD. Also notice that the results in the two tables indicate the scalability of NGBD and NGBD-PCR; the solution time of NGBD or

continuous variables. Details of the MINLP models can be found in Chen et al.37 Again, the case study compares the performance of DICOPT, BARON, NGBD, and NGBD-PCR for the multiperiod MINLPs, with uniform domain partitioning for K = 5, 10, 15. The relative and absolute termination criteria for optimization were set to be 10−2 and the initial value of all the integers was 0. Tables 2 and 3 summarize the results for the problems with 8 and 24 scenarios, respectively. It can be seen that for either problem, DICOPT located a suboptimal solution that is significantly worse an εoptimal solution obtained by NGBD/NGBD-PCR (although it took significantly less computation time). BARON did not return an ε-optimal solution within 106 s and 30 days for the two 7296



Article

NGBD-PCR increased moderately with the number of scenarios, which is because the number of NGBD iterations did not change significantly with the number of scenarios and the number of subproblems to be solved at each iteration only increases linearly with respect to the number of scenarios.

and according to condition eq 11, U conv(xî , x lo , x up) ≤ z r̂ ≤ U conc(xî , x lo , x up), x lo ≤ xî ≤ x up

■

Therefore, ẑr is also feasible for the constraints in eq 5. Proof of Proposition 3. Proof. Considering that Problem ERMP-k differs from Problem RMP-k only with the first group of optimality cuts, it is obvious that Problem ERMP-k is a tighter relaxation if it is a valid relaxation of Problem P. It is proved in ref 15 that Problem RMP-k is a relaxation of Problem P when Problem P is augmented with the relevant integer cuts excluding the previously examined integer realizations, so it remains to prove

CONCLUSIONS This paper develops a piecewise convex relaxation framework and combines it with NGBD for the global optimization of stochastic or multiperiod MINLPs for integrated process design and operation. With the piecewise convex relaxation technique, a tighter lower bounding problem is generated that provides improved information for NGBD to converge. The so-updated NGBD method, NGBD-PCR, solves additional MILP subproblems but solves much fewer nonconvex NLP subproblems which are much more difficult to solve. Therefore, it reduced the solution times of NGBD by up to about an order of magnitude for the case study problems, although NGBD already reduced the solution times by an order of magnitude compared to the state-of-the-art global optimization solver, BARON. The case study results also indicate that the number of partitioned subdomains impacts the efficiency of NGBD-PCR. So it will be an interesting future work to develop a systematic approach to improve the selection of this number. In addition, more flexibility can be introduced for domain partitioning, such as different numbers of subdomains for different variables, nonuniform partitioning of the domains, and ideally an adaptive partitioning scheme that allows the subdomain to be partitioned during the solution procedure according to the intermediate solution information. As the scenarios in the MINLP models are generated via a naive approach in this paper, advanced scenario generation techniques may be introduced to favor a better model for more efficient and effective optimization, such as scenario decomposition,38 scenario reduction,39 and sample average approximation.40

■

(j) objPP(y) ≥ objDPBP − PCR (y(j) , λ ̃ ) s

(j) + (∑ (whchT + (λh̃ )T Bh ))(y − y(j)), h=1 k

∀ j ∈ T̃ , ∀ y ∈ V = {y ̂ ∈ Y : Problem P is feasible for y = y }̂ (14)

Equation 13 implies objPP(y) ≥ objPBP − PCR (y),

objPBP − PCR (y) ≥ sup objDPBP − PCR (y , λ) λ≥0 s

=

∑ sup objDPBP − PCR (y , λh) h

h = 1 λh ≥ 0

(16)

where λ = (λ1, ..., λs). ̃ For all j ∈ T̃ k and h ∈ {1, ..., s}, KKT multipliers λ(j) h ⩾ 0, and

APPENDIX

(j)

objDPBP − PCR (y , λh̃ ) h

Proof of Proposition 1. Proof. For any ẑr that is feasible for the constraints in eq 4, say, corresponding to the McCormick relaxation on the ith subdomain, there exist δî ,x̂i,ŷi such that δî = 1, ẑr = ẑri , xlo ≤ xi ≤ x̂i ≤ xi+1 ≤ xup, xlo ≤ xi ≤ ŷi ≤ xi+1 ≤ xup and ẑ =

(15)

Due to weak duality, for all y ∈ V,

Proofs of Propositions

r

∀y∈V

=

zî r

=

≥ x i + 1yi ̂ + xî y up − x i + 1y up = x i + 1(yi ̂ − y up ) + xî y up ≥ x up(yi ̂ − y up ) + xî y up (noticing x i + 1 ≤ x up , y ̂ ≤ y up )

wh(chTy + inf (xh , qh , δh) ∈ Dhpw (k) + (λh̃ )T (ugpw, h(xh , qh ,

ufpw, h(xh , qh , δh)) δh) + Bh y)

wh(chTy(j) + ufpw, h(xh , qh , δh)) inf (xh , qh , δh) ∈ Dhpw (k) + (λh̃ )T (ugpw, h(xh , qh , δh) + Bh y(j)) (j) + (whchT + (λh̃ )T Bh )(y − y(j))

(j) (j) = objDPBP − PCR (y(j) , λh̃ ) + (whchT + (λh̃ )T Bh )(y − y(j)) h

Therefore, ẑr,x̂i,ŷi satisfy the first constraint of relaxation eq 1. Similarly, ẑr,x̂i,ŷi also satisfy the second to the fourth constraint of relaxation eq 1. Therefore, ẑr is also feasible for the constraints in eq 1. Proof of Proposition 2. Proof. For any ẑr that is feasible for constraints eq 8, say, corresponding to the McCormick relaxation on the ith subdomain, there exist δ̂i,x̂i,ŷi such that δ̂i = 1, ẑr = ẑri and

(17)

Equations 16 and 17 imply that (j)

objPBP − PCR (y) ≥ objDPBP − PCR (y(j) , λ ̃ ) s

+

∑ (whchT + (λh̃(j))T Bh)(y − y(j)),

k

∀ j ∈ T̃ ,

h=1

∀y∈V

U conv, i(xî , 1, x i , x i + 1) ≤ zî r ≤ U conc, i(xî , 1, x i , x i + 1), x i ≤ xî ≤ x i + 1

(18)

and eqs 15 and 18 imply eq 14. 7297



Article

Scenario Generation for Normal Distribution

Notes

The naive sampling rule used to generate scenarios for normal distribution in the case studies is introduced here. If the normal distribution has mean μ and standard deviation σ and the total number of scenarios to be sampled is s, then scenarios are only sampled in the range [−3σ + μ, 3σ + μ] and the hth scenario is

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would like to thank BP for funding this work. This project was conducted as a part of the BP-MIT conversion research program.

rh = −3σ + μ + 3σ /s + (h − 1)6σ /s

■

This is illustrated in Figure 3. Also, the probability of each scenario is

Figure 3. Scenario generation for normal distribution.

⎧ Φ( − 3σ + μ + 6σ /s) if h = 1 ⎪ ⎪ Φ( − 3σ + μ + 6σh/s) − Φ if 1 < h < s P(r = rh) = ⎨ ⎪ ( − 3σ + μ + 6σ(h − 1)/s) ⎪ ⎩1 − Φ( − 3σ + μ + 6σ(s − 1)/s) if h = s

■

where Φ denotes the cumulative distribution function of the normal distribution.

Quadratic functions and cubic functions in the equality constraints introduce nonconvexity to the optimization problem. Since all the variables in the problem are nonnegative, the piecewise relaxation of these functions leads to a set of quadratic constraints in the form of (19)

and cubic constraints in the form of

z ≥ x3

(20)

Then these functions are outer linearized for 5 points in their domains, say [xlo,xup], to yield a MILP lower bounding problem, i.e., constraint eq 19 is relaxed into z ≥ 2xî x − xî 2 , i = 1, ..., 5

(21)

and constraint eq 20 is relaxed into z ≥ 3xî 2x − 2xî 3 , i = 1, ..., 5

(22)

where x̂i = (i − 1)(x − x )/4 + x .

■

up

lo

REFERENCES

(1) Halemane, K. P.; Grossmann, I. E. Optimal process design under uncertainty. AIChE J. 1983, 29, 425−433. (2) Acevedo, J.; Pistikopoulos, E. N. Stochastic optimization based algorithms for process synthesis under uncertainty. Comput. Chem. Eng. 1998, 22, 647−671. (3) Varvarezos, D. K.; Grossmann, I. E.; Biegler, L. T. An outerapproximation method for multiperiod design optimization. Ind. Eng. Chem. Res. 1992, 31, 1466−1477. (4) van den Heever, S. A.; Grossmann, I. E. Disjunctive multiperiod optimization methods for design and planning of chemical process systems. Comput. Chem. Eng. 1999, 23, 1075−1095. (5) Tawarmalani, M.; Sahinidis, N. V. Global optimization of mixedinteger nonlinear programs: A theoretical and computational study. Math. Program. 2004, 99, 563−591. (6) Adjiman, C. S.; Androulakis, I. P.; Floudas, C. A. Global optimization of mixed-integer nonlinear problems. AIChE J. 2000, 46, 1769−1797. (7) Kesavan, P.; Allgor, R. J.; Gatzke, E. P.; Barton, P. I. Outer approximation algorithms for separable nonconvex mixed-integer nonlinear programs. Math. Program., Ser. A 2004, 100, 517−535. (8) Duran, M.; Grossmann, I. E. An outer-approximation algorithm for a class of mixed nonlinear programs. Math. Program. 1986, 66, 327−349. (9) Fletcher, R.; Leyffer, S. Solving mixed integer nonlinear programs by outer approximation. Math. Program. 1994, 66, 327−349. (10) Li, X.; Armagan, E.; Tomasgard, A.; Barton, P. I. Stochastic pooling problem for natural gas production network design and operation under uncertainty. AIChE J. 2011, 57, 2120−2135. (11) Rebennack, S.; Kallrath, J.; Pardalos, P. M. Optimal storage design for a multi-product plant: A non-convex MINLP formulation. Comput. Chem. Eng. 2011, 35, 255−271.

Outer Linearization for Pump Network Configuration Problem

z ≥ x2

LIST OF ACRONYMS BD = Benders decomposition DPBP-PCR = dual of primal bounding problem with piecewise relaxation ERMP = enhanced relaxed master problem (with piecewise relaxation) FP = feasibility problem FRMP = feasibility relaxed master problem GBD = generalized Benders decomposition LBP = lower bounding problem LBP-PCR = lower bounding problem with piecewise relaxation LP = linear program MILP = mixed-integer linear program MINLP = mixed-integer nonlinear program NGBD = nonconvex generalized Benders decomposition NGBD-PCR = nonconvex generalized Benders decomposition with piecewise relaxation PBP = primal bounding problem PBP-PCR = primal bounding problem with piecewise relaxation PP = primal problem RMP = relaxed master problem

lo

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. 7298



Article

(12) Benders, J. F. Partitioning procedures for solving mixedvariables programming problems. Numer. Math. 1962, 4, 238−252. (13) Slyke, R. M. V.; Wets, R. L-Shaped Linear Programs with Applications to Optimal Control and Stochastic Programming. SIAM J. Appl. Math. 1969, 17, 638−663. (14) Geoffrion, A. M. Generalized Benders decomposition. J. Optim. Theory Appl. 1972, 10, 237−260. (15) Li, X.; Tomasgard, A.; Barton, P. I. Nonconvex generalized Benders decomposition for Stochastic Separable Mixed-Integer Nonlinear Programs. J Optim. Theory Appl. 2011, 151, 425−454. (16) Sahinidis, N. V.; Tawarmalani, M. BARON 9.0.4: Global Optimization of Mixed-Integer Nonlinear Programs. In BARON 9.0.4 User’s Manual, 2010. (17) Karuppiah, R.; Grossmann, I. E. Global optimization for the synthesis of integrated water systems in chemical processes. Comput. Chem. Eng. 2006, 30, 650−673. (18) Meyer, C. A.; Floudas, C. A. Global optimization of a combinatorially complex generalized pooling problem. AIChE J. 2006, 52, 1027−1037. (19) Wicaksono, D. S.; Karimi, I. A. Piecewise MILP under- and overestimators for global optimization of bilinear programs. AIChE J. 2008, 54, 991−1008. (20) Gounaris, C. E.; Misener, R.; Floudas, C. Computational comparison of piecewise-linear relaxations for pooling problems. Ind. Eng. Chem. Res. 2009, 48, 5742−5766. (21) Rosen, J. B.; Pardalos, P. M. Global minimization of large-scale constrained concave quadratic problems by separable programming. Math. Program. 1986, 34, 163−174. (22) McCormick, G. P. Computability of global solutions to factorable nonconvex programs: Part I − Convex underestimating problems. Math. Program. 1976, 10, 147−175. (23) Scott, J. K.; Stuber, M. D.; Barton, P. I. Generalized McCormick relaxations. J. Global Optim. 2011, 4, 569−606. (24) Tawarmalani, M.; Sahinidis, N. Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2002. (25) Adjiman, C. S.; Dallwig, S.; Floudas, C. A.; Neumaier, A. A global optimization method, α-BB, for general twice-differentiable constrained NLPs − I. Theoretical advances. Comput. Chem. Eng. 1998, 22, 1137−1158. (26) Gatzke, E. P.; Tolsma, J. E.; Barton, P. I. Construction of convex relaxations using automated code generation technique. Optim. Eng. 2002, 3, 305−326. (27) Balas, E.; Jeroslow, R. Canonical cuts on the unit hypercube. SIAM J. Appl. Math. 1972, 23, 61−69. (28) Moore, R. E. Methods and Applications of Interval Analysis; SIAM: Philadelphia, PA, 1979 . (29) Balas, E. Disjunctive programming and a hierachy of relaxations for discrete optimization problems. SIAM J. Algebraic Discrete Methods 1985, 6, 466−486. (30) Grossmann, I. E. Review of nonlinear mixed-integer and disjunctive programming techniques. Optim. Eng. 2002, 3, 227−252. (31) Grossmann, I. E.; Raman, R.; Kalvelagen, E. DICOPT user’s manual. http://www.gams.com/dd/docs/solvers/dicopt.pdf . (32) IBM, IBM ILOG CPLEX: High-performance mathematical programming engine. http://www-01.ibm.com/software/integration/ optimization/cplex/ . (33) GAMS, General Algebraic and Modeling System. http://www. gams.com/. (34) Westerlund, T.; Pettersson, F.; Grossmann, I. E. Optimization of pump configurations as a MINLP problem. Comput. Chem. Eng. 1994, 18, 845−858. (35) Chen, Y.; Adams, T. A., II; Barton, P. I. Optimal design and operation of static energy polygeneration systems. Ind. Eng. Chem. Res. 2011, 50, 5099−5113. (36) Chen, Y.; Adams, T. A., II; Barton, P. I. Optimal design and operation of flexible energy polygeneration systems. Ind. Eng. Chem. Res. 2011, 50, 4553−4566.

(37) Chen, Y.; Li, X.; Adams, T. A.; Barton, P. I. Decomposition strategy for the Global Optimization of Flexible Energy Polygeneration System. AIChE J. 2011, DOI: 10.1002/aic.13708. (38) Higle, J. L.; Sen, S. Stochastic decomposition: An algorithm for two-stage linear programs with recourse. Math. Oper. Res. 1991, 16, 650−669. (39) Heitsch, H.; Römisch, W. Scenario reduction algorithms in stochastic programming. Comput. Optim. Appl. 2003, 24, 187−206. (40) Kleywegt, A. J.; Shapiro, A.; Homem-De-Mello, T. The sample average approximation method for stochastic discrete optimization. SIAM J. Optim. 2001, 12, 479−502.

7299


Nonconvex Generalized Benders Decomposition with Piecewise

Recommend Documents