Tightening piecewise McCormick relaxations for bilinear problems

G Model CACE-4938; No. of Pages 12

ARTICLE IN PRESS Computers and Chemical Engineering xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Tightening piecewise McCormick relaxations for bilinear problems Pedro M. Castro a,b,∗ a b

Centro de Investigac¸ão Operacional, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal Laboratório Nacional de Energia e Geologia, 1649-038 Lisboa, Portugal

a r t i c l e

i n f o

Article history: Received 6 February 2014 Received in revised form 21 March 2014 Accepted 24 March 2014 Available online xxx Keywords: Optimization Mathematical modeling Nonlinear programming Generalized Disjunctive Programming Water minimization

a b s t r a c t We address nonconvex bilinear problems where the main objective is the computation of a tight lower bound for the objective function to be minimized. This can be obtained through a mixed-integer linear programming formulation relying on the concept of piecewise McCormick relaxation. It works by dividing the domain of one of the variables in each bilinear term into a given number of partitions, while considering global bounds for the other. We now propose using partition-dependent bounds for the latter so as to further improve the quality of the relaxation. While it involves solving hundreds or even thousands of linear bound contracting problems in a pre-processing step, the benefit from having a tighter formulation more than compensates the additional computational time. Results for a set of water network design problems show that the new algorithm can lead to orders of magnitude reduction in the optimality gap compared to commercial solvers. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction The simplest and perhaps most common type of constraint in Chemical Engineering is the blending equation, in which the property of a product resulting from a mix of materials is estimated as a weighted sum by flowrates of the properties of the components. The products of flowrates by properties form bilinear terms that are nonconvex and therefore give rise to multiple local solutions. Blending constraints arise in crude oil operations in refineries (Lee et al., 1996; Jia et al., 2003; Yadav and Shaik, 2012), in the blending of different distilled fractions such as gasoline and diesel (Moro et al., 1998; Jia and Ierapetritou, 2003; Kolodziej et al., 2013b), in the design of distributed wastewater treatment systems (Galan and Grossmann, 1998; Meyer and Floudas, 2006; Teles et al., 2012), integrated water (Karuppiah and Grossmann, 2006; Faria and Bagajewicz, 2012; Rubio-Castro et al., 2013) and mass and property integration networks (Nápoles-Rivera et al., 2010). Bilinear terms also arise in the trim loss problem in paper plants (Harjunkoski et al., 1998; Zorn and Sahinidis, 2013) and in the operation of hydro energy systems (Catalão et al., 2011; Castro and Grossmann, 2014). In order to find rigorous global optimal solutions to bilinear problems, which can be of the nonlinear (NLP) or mixed-integer

∗ Correspondence to: Centro de Investigac¸ão Operacional, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal. Tel.: +351 217500707. E-mail address: [email protected]

nonlinear (MINLP) type, alternative algorithms can be used that have in common the generation of linear (LP) or mixed-integer linear (MILP) relaxations of the original problem. Assuming that we are dealing with a minimization problem, a relaxation provides a lower bound, any feasible solution provides an upper bound and global convergence is achieved when bounds lie within a specified tolerance. It is thus critical to derive tight relaxations to improve the lower bound and also, indirectly, to provide initialization points leading to better solutions. In the standard McCormick (1976) relaxation, each bilinear term is replaced by a new variable and four sets of linear inequality constraints are added to the formulation. In these, the new variable is related to the two variables forming the bilinear term and their lower and upper bounds. As the tighter the lower and upper bounds, the higher the quality of the relaxation, most global optimization solvers do variable bounding to reduce the search space of the relaxed problem. For instance, the new solver GloMIQO (Misener and Floudas, 2013b) uses different techniques such as interval arithmetic, reduced cost and optimality-based bound contraction. In the latter, which is the subject of this paper, a sequence of minimization and maximization problems are solved for each nonlinearly participating variable to find the tightest possible bounds. Faria and Bagajewicz (2011) proposed a global optimization algorithm based on a more thorough optimality-based bound contraction procedure for one of the variables of each bilinear term. It consists of successive contracting steps that first involve solving a relaxation problem to compute reference values for the bilinear participating variables and then on solving multiple auxiliary linear

http://dx.doi.org/10.1016/j.compchemeng.2014.03.025 0098-1354/© 2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model CACE-4938; No. of Pages 12 2

ARTICLE IN PRESS P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

models to update the bounds. Note that direct discretization was used, which has the advantage of requiring two rather than four sets of inequality constraints per bilinear term and the disadvantage of leading to weaker relaxation problems. Further reduction in the feasible space of the relaxed problem requires domain partitioning. It can be done sequentially as in spatial branch and bound (Tawarmalani and Sahinidis, 2005; Misener and Floudas, 2013b) by, for example, bisecting the domain of one of the variables belonging to the bilinear term that most violates the feasible region of the original problem (Ruiz and Grossmann, 2011); or simultaneously. For the simultaneous approach, one can rely on the piecewise McCormick envelopes, proposed by Grossmann and co-workers for the case of univariate and uniform partitioning (Bergamini et al., 2005; Karuppiah and Grossmann, 2006) and extensively studied by, among others, Floudas and co-workers (Gounaris et al., 2009; Misener et al., 2011), Karimi and co-workers (Wicaksono and Karimi, 2008; Hasan and Karimi, 2010), and Faria and Bagajewicz (2012); or on multiparametric disaggregation, a new concept developed by Castro and co-workers (Teles et al., 2013a,b; Kolodziej et al., 2013a). While piecewise McCormick and multiparametric disaggregation are conceptually standalone approaches, proving global optimality for a typically large number of partitions, they can naturally be integrated within spatial branch and bound frameworks as is the case with GloMIQO. In addition, they can replace the standard McCormick relaxation to provide stronger bounds when performing optimality-based bound contraction, at the expense of higher computational times resulting from the solution of multiple MILPs instead of LPs (Castro and Grossmann, 2014). The comparative study between uniform piecewise McCormick and multiparametric disaggregation has shown that the former is tighter when in the presence of quadratic terms but generates considerably larger MILPs for the same number of partitions (Kolodziej et al., 2013a). More precisely, the number of binary variables in piecewise McCormick grows linearly with the number of partitions, whereas for multiparametric disaggregation the growth is only logarithmic. In practice, this is translated into the ability to reach considerably lower optimality gaps due to a computational performance that is often orders of magnitude better, as can be seen for the design of water-using networks (Castro and Teles, 2013). Yet, piecewise McCormick allows for more freedom when choosing the number of partitions, which can be important when they are just a few, more likely to occur for larger problems or when integrated within a spatial branch-and bound solver. In this respect, the piecewise relaxation option in GloMIQO (Misener and Floudas, 2013b), although turned off by default since convergence is generally faster that way, increases the probability that the bilinear problem will solve to ␧-global optimality within a given time limit. Experiments were conducted for 2, 4 and 8 partitions and for a maximum number of partitioned variables equal to 30 (Misener and Floudas, 2012). The novelty of this work is to propose tighter constraints for the piecewise McCormick relaxation of a bilinear term that feature improved lower and upper bounds for the partitioned variable as well as for the non-partitioned variable. Partition-dependent bounds for the non-partitioned variable are determined through optimality-based bound contraction, which involves solving two linear problems per partition. Once the solution of the tightened relaxation approach is found, it is used as a starting point for the solution of the original problem so as to find a near optimal solution with a fast local solver and to compute a rigorous optimality gap. To illustrate the performance of the new algorithm, we will consider small test problems and larger water-using network design problems from the literature. While such problems are of the nonlinear type, it should be noted that mixed-integer bilinear problems can also be tackled by the proposed approach simply by solving MILPs

rather than LPs when performing optimality-based bound contraction and by fixing the binary variables before using a local NLP solver to find a feasible solution (see Castro and Grossmann, 2014, who consider standard McCormick and multiparametric disaggregation relaxation approaches). Relaxation techniques for bilinear terms with binary variables and nonlinear integer problems with signomial terms can also be found in Rodriguez and Vecchietti (2013) and Tsai and Lin (2013). 2. Problem definition We consider the class of nonconvex, nonlinear problems where all nonlinear terms are of the bilinear type xi xj and all variables are continuous. In (P), x is an m-dimensional vector of non-negative variables that lie between given lower xL and upper xU bounds. Set Q includes all functions fq , including the objective function f0 and all the constraints. The function hq (x) is linear in x, BL is an (i,j)index set that defines the bilinear terms xi xj present in the problem and aijq is a scalar. Note that i = j can be allowed to accommodate quadratic problems. min x = f0 (x)

(P)

subject to fq (x) ≤ 0 ∀q ∈ Q \{0} fq (x) = xL ≤ x



aijq xi xj + hq (x) ∀q ∈ Q

(i,j) ∈ BL ≤ xU

x ∈ Rm 3. Lower bounding formulation from piecewise McCormick envelopes McCormick (1976) envelopes can be used to provide different types of relaxations to problems of type (P). In the standard approach, a linear programming (LP) relaxation is derived by replacing each bilinear term involving variables xi and xj with a new variable wij = xi xj and adding four sets of constraints. A tighter mixed-integer linear programming (MILP) relaxation can be constructed by partitioning the domain of one of the variables (xj ) of the bilinear term into n disjoint regions, with new binary variables being added to the formulation to select the optimal partition for xj . In the seminal work by Grossmann and co-workers (Bergamini et al., 2005; Karuppiah and Grossmann, 2006), partitions are generated uniformly and there is a linear increase in problem size with the number of partitions. This is the lower bounding approach considered in this work. A different uniform partitioning scheme featuring a logarithmic growth in the number of binary variables was proposed by Misener et al. (2011) but the results failed to show major benefits in their use. Other extensive studies involving a variety of piecewise under- and overestimators, some of them featuring non-identical segment lengths, have been performed by Wicaksono and Karimi (2008) and Gounaris et al. (2009). Hasan and Karimi (2010) report on the advantage of uniform placement of grid points. L and xU represent respectively the lower and upper Let xjn jn bounds of variable xj for partition n. If the value of xj does belong to such partition, then binary variable yjn = 1 and the McCormick envelopes hold. The piecewise McCormick relaxation (PR-GDP) can be formulated as a Generalized Disjunctive Program (Raman and Grossmann, 1994) and is tighter due to the use of the L and xU in the four constraints partition-dependent parameters xjn jn inside the disjunction, instead of the global bounds xjL and xjU .

min z R = f0 (x) =



aij0 wij + h0 (x)

(PR-GDP)

(i,j) ∈ BL

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model CACE-4938; No. of Pages 12

ARTICLE IN PRESS P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

3

subject to

Remark 1. The quality of the relaxation is influenced by the choice of variables xj selected for partitioning. Hasan and Karimi (2010) propose a very simple integer program to determine the minimum cardinality set of partitioned variables but as can be seen in Section 5.1 for test problem P3, selecting the fewest partitioning variables is not necessarily the best option. In general, there may exist an optimal set of variables that would result in the best relaxation but it is beyond the scope of this paper to find such set. This is because in most engineering problems, as the one considered in section 6, bilinear terms involve two different sets of variables (e.g. flowrates, qualities, split fractions) and so there are just a couple of obvious choices. For instance, Hasan and Karimi (2010) chose to partition: flow variables in problems involving column sequencing for nonsharp distillation and integrated water use and treatment systems; qualities in the generalized pooling problem on wastewater treatment networks; and areas and split fractions when synthesizing heat exchanger networks. Remark 2. In this paper, we consider the case of univariate partitioning. In bivariate partitioning the range of both variables forming the bilinear term is known a priori, leading to a relaxation that is usually tighter (Hasan and Karimi, 2010). In the Appendix, we provide the piecewise McCormick relaxation for the case of bivariate partitioning. The results for the motivating example show that the approach to be presented next is more efficient.

4. Improving the quality of the piecewise relaxation To further improve the quality of the piecewise McCormick relaxation, we now propose partition-dependent lower and upper bounds also for variable xi . The difference compared to xj is that the tighter bounds are not easily computed using information from the global bounds and the number of partitions N as can be seen in (PR-GDP). Instead, they will be determined through optimalitybased bound contraction, a more time consuming procedure to be explained in the next section. For now, let us assume that paramL and xU , representing the lower and upper bounds of eters xijn ijn variable xi when variable xj is constrained to partition n, are known. The tightened piecewise relaxation approach (PRT-GDP) is given below. Notice that the changes occur inside the disjunction, which now features the bounding constraints of variable xi .

The two most common alternatives to reformulate a linear GDP into a MILP are big-M and convex hull relaxations (Balas, 1985; Grossmann and Trespalacios, 2013). It is well known that the big-M reformulation can exhibit poor relaxation quality and for the specific case of piecewise partitioning relaxation schemes it is not competitive (Wicaksono and Karimi, 2008). We thus apply a convex hull reformulation, along the lines of Karuppiah and Grossmann (2006), Meyer and Floudas (2006) and Castro and Teles (2013). The full lower bounding MILP formulation (PRT-MILP) is given below. Notice the appearance of new sets of disaggregated variables xˆ ijn and xˆ jn .



min z R = f0 (x) =

aij0 wij + h0 (x)

(PRT-MILP)

(i,j) ∈ BL subject to



fq (x) =

aijq wij + hq (x) ≤ 0 ∀q ∈ Q \{0}

⎫  ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ n=1 ⎪ ⎪ N ⎪   ⎪ ⎪ U U U U ⎪ xˆ ijn · x + x · xˆ jn − x · x · yjn wij ≥ ⎪ jn ijn ijn jn ⎪ ⎪ ⎪ n=1 ⎪ ⎪ N   ⎬ (i,j) ∈ BL N



wij ≥

xˆ ijn · xL + xL · xˆ jn − xL · xL · yjn jn ijn ijn jn

xˆ ijn · xL + xU · xˆ jn − xU · xL · yjn

wij ≤

⎪ ⎪ ⎪ ⎪ ⎪  ⎪  ⎪ xˆ ijn · xU + xL · xˆ jn − xL · xU · yjn ⎪ wij ≤ ⎪ ijn ijn jn jn ⎪ ⎪ ⎪ n=1 ⎪ ⎪ N ⎪  ⎪ ⎪ ⎪ xi = xˆ ijn ⎪ ⎭ n=1 ⎫ N  ⎪ ⎪ xˆ jn ⎪ xj = ⎪ ⎬ jn

ijn

ijn

jn

∀ (i, j) ∈ BL

n=1 N

N 

n=1

yjn = 1



n=1

xL jn

⎪ ⎪ ⎪ ⎪ ⎭

=

xL j

+

∀ {j|(i, j) ∈ BL}

xU − xL



j

j

n (n − 1)

N

xU = xL +

(xU − xL )nn j

j

⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭

∀ {j|(i, j) ∈ BL}, n ∈ {1, . . ., N}

j jn N xL · yjn ≤ xˆ jn ≤ xU · yjn jn jn xL · yjn ≤ xˆ ijn ≤ xU · yjn ∀ (i, j) ∈ BL, n ∈ {1, . . ., N} ijn

ijn

xL ≤ x ≤ xU x ∈ Rm yjn ∈





0, 1

∀

j| (i, j) ∈ BL



, n ∈ {1, . . ., N}

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model

ARTICLE IN PRESS

CACE-4938; No. of Pages 12

P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

4

Remark 3. For N = 1, the added disaggregated and binary variables can be removed from the formulation so that (PRT-MILP) is reduced to the standard McCormick linear relaxation of problem (P). Remark 4. In (PRT-MILP), it is assumed that the number of partitions N is the same for all partitioned variables xj but it is straightforward to consider otherwise. In fact, while starting with the same number, the bound contraction algorithm to be described next will typically eliminate a different number of partitions for each partitioned variable. To take advantage of this, the number of partitions in the actual implementation of (PRT-MILP) is variable dependent. The convex hull reformulation of (PR-GDP) that uses global bounds for variables xi is named (PR-MILP) and will be used for comparative purposes. 4.1. Generating tighter bounds for the non-partitioned variables Replacing the global bounds xiL and xiU of the non-partitioned L and xU , will result in a tighter variables xi with better bounds xijn ijn

relaxation, leading to a zR value that is closer to the global optimal solution. The improved bounds can be obtained through optimality-based bound contraction. LP formulation (BC) assumes than an upper bound z’ on z is known, which can be obtained by solving (P) with a local NLP solver. By constraining in (BC) the objective function of (P), we limit the domain of xi∗ to regions that can actually improve the quality of the solution. To determine the lower bound xiL∗ j∗ n∗ through (BC) we need first to pick a partitioned vari-

able j∗ |(i∗ , j∗ ) ∈ BL and a partition n∗ ∈ {1, . . ., N}. The bounds for xj∗ are set to the partition limits, xjL∗ = xjL∗ n∗ and xjU∗ = xjU∗ n∗ , whereas global bounds are used for all other bilinear term variables. These bounds appear in the standard McCormick envelopes that are used to provide the relaxation of the bilinear terms xi · xj . To determine the upper bound xiU∗ j∗ n∗ , the constraints remain the same but the objective function changes from minimizing to maximizing variable xi∗ . xiL∗ j∗ n∗ := min xi∗

(xiU∗ j∗ n∗ := max xi∗ )

subject to f0 (x) =





aij0 wij + h0 (x) ≤ z 

(i,j) ∈ BL0

fq (x) =

aijq wij + hq (x) ≤ 0 ∀q ∈ Q \{0}

(i,j) ∈ BLq wij ≥xi · xjL + xiL

· xj − xiL · xjL

wij ≥xi · xjU + xiU · xj − xiU · xjU wij ≤ xi · xjL + xiU · xj − xiU · xjL wij ≤ xi · xjU + xiL · xj − xiL · xjU xjL∗ = xjL∗ n∗ ≤ xj∗ ≤ xjU∗ n∗ = xjU∗

⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭

(BC)

∀ (i, j) ∈ BL

xL ≤ x ≤ xU x ∈ Rm Remark 5. A particular problem (BC) is not necessarily feasible. In such case, partition n* of variable xj∗ can be eliminated from the domain of some summations and constraints in (PRT-MILP) leading to a smaller problem size for the latter. More specifically, binary variables yj∗ n∗ and continuous variables xˆ j∗ n∗ and xˆ ij∗ n∗ , will no longer be required. Since the feasible region is the same for the minimization and maximization problems, if the former is found infeasible, there is no need to solve the latter. Remark 6. Problems (BC) are in essence independent and so can be solved in parallel (Quesada and Grossmann, 1995). However, in the optimality-based bound contraction algorithm described in

Fig. 1. Optimization algorithm featuring new bound contraction approach.

the next section, they will be solved sequentially, with the global bounds being updated whenever possible. 4.2. Bound contraction algorithm In Fig. 1, we show the proposed optimization algorithm featuring the tightened piecewise McCormick bound contraction

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model CACE-4938; No. of Pages 12

ARTICLE IN PRESS P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

procedure. Assuming that problem (P) is feasible and that the local solver can converge following initialization with the values from (PRT-MILP), the algorithm is capable of finding a local optimal solution to the problem and compute an optimality gap, given the bilinear terms BL, the bounds on variables xi and xj and the number of partitions N. The algorithm starts by solving (P) in order to find the upper bound z required by the bound contraction formulations. In case (P) is found infeasible due to a poor initialization, one can continue with z  = +∞ and still improve the variable bounds. The next step is to solve one minimization and one maximization bound contraction problem for each variable appearing in a bilinear term. In this stage, we use the standard McCormick relaxation leading to linear programming (LP) problems identical to (BC) without constraint xjL∗ = xjL∗ n∗ ≤ xj∗ ≤ xjU∗ n∗ = xjU∗ . Doing bound contraction as a pre-processing step takes time but can in general lead to a major reduction in the global domain of the variables. As a consequence, in subsequent steps, the domain of the partitioned variables in every partition n will be shorter, and the relaxation will be tighter, meaning a lower optimality gap. The initialization of the tightened bound contraction procedure starts by picking the first partitioned variable j∗ = First(J). We then select the first xi variable that forms a bilinear term with xj∗ : i∗ = First(i ∈ I|(i∗ , j∗ ) ∈ BL); and choose the first partition n∗ = 1. After constraining xj∗ to the partition bounds, two LP problems (BC) are solved. If feasible, we have the tightened bounds for variable xi∗ when constrained to partition n∗ of variable xj∗ . In case the problem is infeasible, such partition can be discarded. The algorithm proceeds to the next iteration by making n∗ = n∗ + 1 and continues solving (BC) problems up to n∗ = N. It is then time to update the global bounds of variable xi∗ . The lower bound is set to the minimum over active partitions n of lower bound xiL∗ j∗ n , while xiU∗ = max n

xiU∗ j∗ n . This can be done since the optimal value of xj∗

must belong to one such partition. Once we reach the last variable xi participating in a bilinear term with xj∗ , i.e. i∗ = Last(I), we can update the bounds on xj∗ . Tighter bounds will be generated whenever one of the extreme partitions is found to be infeasible when solving (BC). Once all partitioned variables xj are handled, the tightened piecewise McCormick formulation (PRT-MILP) can be solved to generate a lower bound zR to the solution of (P). In case the problem is not solved to optimality, the lower bound is set to the best possible solution at the maximum computational time limit termination. The algorithm terminates with the second run of (P) with a local solver, which should provide a better feasible solution z’ due to a much better initialization (solution from the relaxed problem (PRT-MILP), solved immediately before), although not necessarily the global optimal solution.

5

5. Motivating examples We first consider the very simple example problem reported by Al-Khayyal and Falk (1983) and solved to global optimality by Quesada and Grossmann (1995) and Kolodziej et al. (2013a). In particular, the latter reference provides the full set of constrains for (PR-MILP). Variable x1 is selected as the partitioned variable and so BL = {(2,1)}. min f0 (x) = −x1 + x1 x2 − x2 subject to −6x1 + 8x2 ≤ 3 3x1 − x2 ≤ 3

(P1)

0 ≤ x1 , x2 ≤ 1.5 While the variables domain belongs to [0,1.5], tighter bounds can be inferred after applying the standard optimality-based bound contraction. Solving NLP (P1) with a local solver from starting point (0,0) generates what is actually a suboptimal solution z  = f0 (x) = −1. This is information required by the bound contraction subproblems that are able to reduce the domain to 0.357143 ≤ x1 ≤ 1.375 and 0 ≤ x2 ≤ 1.26. Let us now assume N = 3 partitions for x1 , corresponding to L = 0.357143, xU = xL = 0.696429, xU = xL = 1.035714 and x11 11 12 12 13 U = 1.375. Solving (BC) for j∗ = 1, i∗ = 2 and n∗ = 1 (check Fig. 1), x13 one realizes that the LP is infeasible and so the first partition can be eliminated. The other two partitions define regions that can L = 0, actually improve the current upper bound, leading to x212 U U L x212 = 1.137609, x213 = 0.107143 and x213 = 1.195150. The subsequent solution of MILP model (PRT-MILP) generates a lower bound z R = −1.169619, which is tighter than the value (−1.185207) resulting from (PR-MILP) that uses global bounds for x2 . The optimal solution returned by the NLP solver is now z = −1.083333, obtained at (1.166667, 0.5), which is the global optimal solution. The optimality gaps are thus 7.38% and 8.60% respectively. While the difference in optimality gap may not seem significant, it tends to increase with the number of partitions as can be seen in Table 1. More precisely, we get reductions of one and two orders of magnitude for 150 and 1500 partitions. The drawback from the (PRT-MILP) approach results from the need to solve multiple (BC) problems, up to 2N LPs if all partitions were found to be active, which makes the total computational time increase by respectively one and two orders of magnitude. The interesting part is that (PRTMILP) is an overall better approach, since solving (PR-MILP) for 7500 partitions takes more or less the same time than solving (PRTMILP) for 1500 partitions but the optimality gap is 19 times higher. Finally, in order to identify the prevailing effect in the reduction of

Table 1 Key computational statistics for (P1) for different number of partitions. Relaxation approach

Partitions (N)

15

150

1500

7500

(PR-MILP)

Binary variables Total variables Equations Total CPUs Lower bound Optimality gap

15 49 70 0.93 −1.10173 1.67%

150 454 610 1.60 −1.085376 0.188%

1500 4504 6010 3.63 −1.083538 0.0189%

7500 22,504 3010 302 −1.083374 0.0038%

(PRT-MILP)

Optimality gap Lower bound Total CPUs Binary variables Total variables Equations

0.77% −1.09175 3.63 9 31 46

0.010% −1.083446 27.8 75 229 310

0.0002% −1.083335 267 739 2221 2966

(PR-MILP) with active partitions only

Lower bound Optimality gap

−1.10047 1.56%

−1.08522 0.174%

−1.083522 0.0174%

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model

ARTICLE IN PRESS

CACE-4938; No. of Pages 12

P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

6 Table 2 Key computational statistics for (P2). Relaxation approach

Partitions (N)

9

90

900

2700

(PR-MILP)

Binary variables Total variables Equations Total CPUs Lower bound Optimality gap

18 51 71 0.83 53.85311 8.41%

180 456 557 1.2 57.89790 0.84%

1800 4506 5417 7.76 58.33503 0.083%

5400 13,506 16,217 478 58.36747 0.0278%

(PRT-MILP)

Optimality gap Lower bound Total CPUs Binary variables Total variables Equations

0.84% 57.89737 3.33 5 18 31

0.02% 58.37400 24.0 32 81 103

0.0002% 58.38353 223 297 694 799

optimality gap, we solved (PR-MILP) considering only the active partitions (problem size equal to that from (PRT-MILP)). As can be seen in Table 1, the benefit from fewer partitions is minor, leading to the conclusion that the effect from the partition dependent bounds L and xU is the dominant one. xijn ijn

and (c) J = {4, 5, 6, 7, 8}. The global optimal solution is f0 (x) = 7049.2480.

5.1. Other test problems

0.0025 (−x4 + x5 + x7 ) − 1 ≤ 0 0.01 (−x5 + x8 ) − 1 ≤ 0

We consider two other small test problems. Problem (P2) is taken from Floudas et al. (1999) and has been thoroughly studied by Teles et al. (2013a,b). Since it involves a bilinear and two quadratic terms, both x1 and x2 need to be partitioned. In this case, applying the standard bound contraction results in 1.212711 ≤ x1 ≤ 6.220287 and 1.636883 ≤ x2 ≤ 6.180705, when considering as upper bound z the global optimal solution f0 (x) = 58.383672.

min f0 (x) = x1 + x2 + x3 subject to 0.0025 (x4 + x6 ) − 1 ≤ 0

(P3)

100x1 − x1 x6 + 833.33252x4 − 83333.333 ≤ 0 x2 x4 − x2 x7 − 1250x4 + 1250x5 ≤ 0 x3 x5 − x3 x8 − 2500x5 + 1, 250, 000 ≤ 0 100 ≤ x1 ≤ 10, 000, 1000 ≤ x2 , x3 ≤ 10, 000

min f0 (x) = 6x12 + 4x22 − 2.5x1 x2

10 ≤ x4 , x5 , x6 , x7 , x8 ≤ 1000

subject to

The results in Table 2 for (P2) show a similar trend to those obtained for (P1). Nevertheless, the reduction in problem size is more pronounced, with (PRT-MILP) requiring one sixth of the binary variables due to a more extent elimination of partitions (the number of partitions is directly related to the number of binary variables). For (P3), the impact of the new approach in optimality gap is not as significant, in part because fewer partitions can be discarded. The most interesting result in Table 3 is that the total computational time for (P3a) and (P3b) for N = 200 and (P3c) for N = 100, is lower for (PRT-MILP) than (PR-MILP). It tells us that the

x1 x2 − 8≥0

(P2)

1 ≤ x1 , x2 ≤ 10

Problem (P3) can be found in Hock and Schittowski (1981), Teles et al. (2013b), and is particularly challenging for the piecewise McCormick approach (Kolodziej et al., 2013a). There are now 5 bilinear terms: x1 x6 , x2 x4 , x2 x7 , x3 x5 and x3 x8 ; corresponding to three partitioning alternatives: (a) J = {1, 2, 3}; (b) J = {2, 5, 6, 8}; Table 3 Key computational statistics for (P3). Alternative

Relaxation approach

Partitions (N)

10

100

200

(P3a)

(PR-MILP)

Binary variables Total CPUs Optimality gap Optimality gap Total CPUs Binary variables

30 1.80 53.1% 27.1% 10.8 30

300 43.7 3.8% 1.7% 110 285

600 1021 1.96% 0.83% 853 568

Binary variables Total CPUs Optimality gap Optimality gap Total CPUs Binary variables

40 1.78 40.7% 35.1% 10.2 36

400 71.1 3.0% 2.2% 144 325

800 1095 1.30% 0.94% 792 648

Binary variables Total CPUs Optimality gap Optimality gap Total CPUs Binary variables

50 1.87 26.4% 23.6% 10.2 46

500 253 1.7% 1.5% 194 422

1000 2100 0.72% 15.6% 3771a 842

(PRT-MILP)

(P3b)

(PR-MILP)

(PRT-MILP)

(P3c)

(PR-MILP)

(PRT-MILP)

a

Maximum resource limit of 3600 CPUs reached when solving (PRT-MILP), after bound contracting stages.

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model

ARTICLE IN PRESS

CACE-4938; No. of Pages 12

P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

7

Table 4 Effect of improved bounds on key computational statistics for water network problems (suboptimal solutions in italic, lowest gap in bold, termination due to maximum computational limit in bold italic). Non-partitioned variables (|I|)

Partitioned variables (|J|)

Partitions N

Ex2

20

16

10 50

Ex3

20

16

10 100

Ex6

30

25

Ex7

30

Ex8

(PR-MILP)

(PRT-MILP)

Solution

Total CPUs

Gap

Gap

Total CPUs

Solution

74.4699

18.8 3609

0.1250% 0.1245%

0.0484% 0.0103%

141 1536

74.4699

143.4126

10.2 453

0.0037% 0.0003%

0.0020% 0.0000%

10 50

142.0816

44.7 3614

0.0435% 0.1379%

0.0110% 0.0041%

205 1090

142.0816

30

10 50

280.7712

30.3 1848

0.0799% 0.0198%

0.0371% 0.0076%

287 1671

280.7712

30

20

10 100

164.4898

18.8 3613

0.7134% 1.9508%

0.2132% 0.0130%

207 1803

164.4898

Ex9

42

36

10 20

312.922

142 1552

0.1989% 0.1148%

0.1253% 0.0537%

460 1029

312.9215

Ex10

15

9

10 100

169.1173

0.0160% 0.0023%

0.0008% 0.0000%

45.3 323

169.1173

Ex11

15

12

10 50

104.8861

0.1946% 0.0322%

0.1698% 0.0237%

122 824

104.8861

Ex12

20

8

10 200

165.1953

0.9872% 0.0438%

0.4827% 0.0116%

58.9 892

165.1953

Ex13

28

16

10 40

178.2629

31.0 3612

0.2559% 0.1645%

0.2197% 0.1152%

201 1165

178.2629

Ex14

32

24

10 20

329.9689 329.9566

282 3615

1.2259% 2.3368%

1.0499% 0.9206%

447 4270

329.5698

Ex15

40

30

10 20

361.5177 361.6856

231 3618

0.7963% 0.9361%

0.7427% 0.7779%

723 4486

361.5177 361.6856

Ex16

24

16

10 50

285.9343

17.4 3610

1.1105% 1.5755%

0.8861% 0.0526%

212 2042

285.9343

16

0.0000%

0.0000%

149

157.0944

0.0372% 0.0035%

0.0372% 0.0035%

108 998

238.7333

0.9546% 0.9691%

0.9507% 0.7029%

625 4615

403.1960

Ex17

24

12

10

157.0944

Ex18

15

12

10 100

238.7333

Ex20

40

20

10 20

403.1960

former MILP is significantly easier to solve, more than compensating the additional time taken by the thorough bound contracting procedure executed before, which is exclusive to (PRT-MILP). Thus, there is the double advantage of lower computational time and smaller optimality gap. On the other hand, the MILP solver had more difficulties with (P3c) for N = 200, being unable to close the integrality gap within the hour. It explains the higher optimality gap in italic in Table 3, since given enough time the gap should at least be the same as the one for (PR-MILP). 6. Computational results for larger test problems To illustrate the performance of the proposed approach on a set of larger NLPs, we consider the 32 most challenging problems involving the optimal design of water-using networks with respect to freshwater minimization. Missing from the original test set (data provided in the supplementary material of Teles et al., 2009) are problems Ex1, Ex4 and Ex5, which can be solved quite easily by a global optimization algorithm relying on piecewise McCormick relaxation (Castro and Teles, 2013), whereas Ex19 is indistinguishable from Ex2. These problems have been solved by heuristic decomposition strategies (Teles et al., 2009) and other global optimization solvers (Teles et al., 2012; Kolodziej et al., 2013a; Misener and Floudas, 2013a), and for some of them, the global optimal solution is yet to be proven. It should be

6.69 809 12.5 1947 7.11 707

8.18 190 88.8 3617

93.5 938

143.4126

highlighted that bilinear terms arise from the product of flowrates by outlet concentrations and so there are two natural choices for selecting the partitioning variables, either the set of flowrate or concentration variables. Results in Castro and Teles (2013) have shown that it is preferable to select the concentration variables for partitioning, and so this is the approach followed in this paper. The algorithm and mathematical formulations have been implemented in GAMS 24.1.3. The resulting LP and MILP problems were solved using CPLEX 12.5.1, whereas the local solver for the NLPs was CONOPT 3. Global solvers GloMIQO 2.3 and BARON 12.3.3 have also been used for comparison. For the results in Sections 5 and 6.1, the hardware consisted on an Intel Core i7 950 (3.07 GHz), whereas for Section 6.2, we have used an Intel Core i7-3770 (3.40 GHz) processor. In both cases, there was 8 GB of RAM, the operating system was Windows 7 Professional and we have used a single thread and default options except for the 10−6 relative optimality tolerance. The other termination criterion was a maximum computational time for problems generated by (P), (PR-MILP) and (PRT-MILP) equal to 3600 CPUs. Such cases are identified in italic. It is thus possible for the total computational time listed in the tables to exceed 3600 CPUs, even by a wide margin. This occurs when solving hundreds or even thousands of (BC) subproblems and despite the fact that in some cases (not in italic) problem (PRT-MILP) can be solved to optimality within the hour.

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model

ARTICLE IN PRESS

CACE-4938; No. of Pages 12

P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

8

Table 5 Effect of improved bounds on key computational statistics for water network problems. Non-partitioned variables (|I|)

Partitioned variables (|J|)

Partitions N

(PR-MILP)

(PRT-MILP)

Solution

Total CPUs

Gap

Gap

Total CPUs

Solution

Ex21

28

24

10 30

216.3702

49.1 3613

0.1823% 0.2146%

0.1613% 0.0515%

267 1357

216.3702

Ex22

32

24

10 20

323.6816 323.5051

683 3614

0.3210% 0.3330%

0.2106% 0.5799%

592 4466

323.5051

Ex23

40

30

10 30

366.4735

97.8 3619

0.2500% 1.5827%

0.0992% 0.0462%

371 1631

366.4735

Ex24

32

24

10 30

511.9642

30.2 1660

0.1866% 0.0770%

0.0773% 0.0305%

311 1075

511.9642

Ex25

40

15

10 15

410.6354

448 3613

0.6926% 1.6928%

0.6825% 2.3181%

798 3984

410.6354

Ex26

130

60

5 10

797.9772

461 3694

0.9877% 0.9177%

0.9547% 0.5362%

2024 7111

797.9772

Ex27

136

32

5 10

556.6752

313 3665

2.9393% 2.2458%

2.9386% 2.3105%

1377 5272

556.6752

Ex28

420

120

2 5

1812.170

2809 4946

2.2508% 3.2003%

2.2363% 2.8960%

17,314 38,933

1812.1703

Ex29

294

84

2 5

1285.006 1285.090

1652 4301

0.8487% 1.2948%

0.8414% 1.2808%

8003 19,079

1285.006

Ex30 Ex31 Ex32

315 420 420

90 120 80

2 2 2

751.4707 694.3594 639.2615

4358 5486 4441

11.312% 9.0220% 10.282%

11.962% 7.4028% 11.993%

11,230 23,356 9800

753.6011 683.7706 639.2615

Ex33

187

44

2 5

572.9515 574.2588

198 3713

14.758% 7.4202%

14.753% 7.5227%

879 5338

572.9515 567.3166

Ex34

221

52

2 5

595.2940 587.0710

268 3738

9.8073% 3.7920%

9.6669% 6.4956%

962 5446

595.2940 589.1558

Ex35

221

65

5 10

783.9510 783.9520

1915 3867

0.6415% 1.3879%

0.6251% 1.4957%

6330 11,615

783.9510 783.9520

Ex36

126

45

5 10

662.8070

1039 3671

2.1435% 3.7501%

2.1422% 2.0689%

2695 5991

662.8070

6.1. Tightened vs. conventional piecewise relaxation approach The new algorithm, featuring (PR-MILP) or (PRT-MILP), is run for two different N values, unless global optimality can be proven or the problem is already too complex when considering just two partitions. With respect to the selection of the number of partitions, we have used N = 10 by default since above this value, piecewise McCormick relaxation becomes less efficient than multiparametric disaggregation (Kolodziej et al., 2013a). A ten-fold increase typically leads to an order of magnitude reduction in optimality gap (Castro and Teles, 2013), making N = 100 an obvious choice for a second run. However, this was often too much and so the value needed to be tuned. The results have been divided in two tables, Tables 4 and 5, with the former featuring problems that are on average easier to solve. Complexity is very much related to the number of bilinear terms, which for this particular problem is equal to the product of the nonpartitioned variables |I| and partitioned variables |J|, as well as to the number of partitions. Problems with more bilinear terms typically feature fewer partitions so that the resulting MILP can still provide useful information in one hour. Increasing the number of partitions always leads to a reduction in the optimality gap, provided that the solver does not terminate due to a maximum time limit. In such cases, the lower bound is computed as the best possible solution at termination, which may be worse than the lower bound obtained from using fewer partitions (it explains for instance why the optimality gap for Ex15 is higher for N = 20 than for N = 10). Lower optimality gaps provide better initialization points, making it more likely for the local NLP solver to find the global optimum solution.

From the results in Table 4, it can be seen that the new tightened approach (PRT-MILP) always provides lower or equal optimality gaps than its conventional piecewise McCormick (PR-MILP) counterpart. Only for Ex17 and Ex18 are the gaps the same, being better to use (PR-MILP) since it does not need to solve multiple bound contracting (BC) problems. From the perspective of total computational time, (PRT-MILP) is faster in 56% of the cases (considering the finer partition level). Thus, the time taken on solving LP problems to obtain a better relaxation is well spent. The algorithm based on the new approach has a single failure in terms of finding the optimal solution (Ex15 for N = 20), while (PR-MILP) fails three times. There are also 5 more cases for which problems resulting from (PR-MILP) cannot be solved in one hour (and in Ex6, Ex8 and Ex16 this means that the optimality gap could not be lowered when increasing the number of partitions N), whereas those from (PRT-MILP) not only can but may lead to an order of magnitude reduction in the optimality gap (Ex8 and Ex16). Overall, Ex3, Ex10 and Ex17 can be solved to global optimality (relative tolerance = 10−6 ) and the highest gap is obtained for Ex14 at 0.9206%. For the more complex problems in Table 5, similar trends are observed. With respect to finding the global optimal solution, (PR-MILP) fails 10 times, compared to 8 failures for (PRT-MILP). Notice that for Ex30, Ex31 and Ex32, the minimum number of partitions (N = 2) already leads to MILP problems that cannot be solved to optimality in one hour. Still, considerable lower gaps are returned compared to commercial solvers (see Table 6). In problems with over 60 partitioned variables, it can be argued that the new tightened approach leads to an unreasonable number of (BC) problems to solve, increasing the total

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model

ARTICLE IN PRESS

CACE-4938; No. of Pages 12

P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

9

Table 6 Optimality gap and total computational effort for piecewise relaxation approaches and commercial solvers (gap within 1% of the lowest gap in bold, termination due to maximum computational limit in italic). Optimum

EX2 EX3 EX6 EX7 EX8 EX9 EX10 EX11 EX12 EX13 EX14 EX15 EX16 EX17 EX18 EX20 EX21 EX22 EX23 EX24 EX25 EX26 EX27 EX28 EX29 EX30 EX31 EX32 EX33 EX34 EX35 EX36 a b c d

74.470 143.413 142.082 280.771 164.490 312.922 169.117 104.886 165.195 178.263 329.570 361.518 285.934 157.094 238.733 403.196 216.370 323.505 366.473 511.964 410.635 797.977 556.675 1812.170 1285.006 744.112b 681.583b 638.717b 565.343 583.522 783.951 662.807

Partitions N

Optimality gapa (PR-MILP)

(PRT-MILP)

GloMIQO 2.3

BARON 12.3.3

58 58 25 21 31 13 135 101 114 42 25 16 48 64 101 24 28 25 16 25 31 4 6 2 2 2 2 2 4 3 3 5

0.027% 0.0007% 0.010% 0.086% 0.265% 0.247% 0.001% 2.678% 0.065% 0.117% 1.577% 0.657% 0.165% 0.0000% 0.563% 1.311% 0.099% 1.458% 0.162% 0.081% 6.570% 1.261% 2.389% 2.251% 0.849% 9.821% 7.016% 9.958% 6.725% 4.662% 1.598% 2.143%

0.050% 0.0001% 0.004% 0.040% 0.119% 0.137% 0.0000% 0.098% 0.023% 0.066% 1.358% 0.339% 0.099% 0.0000% 0.001% 0.608% 0.045% 0.104% 0.082% 0.030% 4.766% 1.228% 2.387% 2.236% 0.841% 10.31% 7.059% 11.48% 6.715% 4.598% 1.584% 2.142%

0.237% 0.383% 0.500% 0.0001% 0.0001% 0.0001% 0.0001% 1.302% 0.873% 0.0001% 0.620% 0.0001% 0.0001% 3.927% 1.676% 3.460% 0.052% 0.033% 0.010% 0.0001% 5.843% 3.385% 11.24% 7.908% 3.753% 113.1% 267.2% 137.4% 40.24% 45.14% 3.770% 13.31%

0.308% 0.247% 0.285% 0.011% 0.0001% 0.094% 0.019% 0.120% 0.704% 0.530% 1.577% 2.185% 0.001% 2.274% 0.997% 2.982% 0.039% 0.225% 0.038% 0.014% 6.435% 3.089% 7.898% 6.644% 3.425% 120.0% 276.9% 99.75% 31.41% 25.94% 2.824% 11.76%

Total CPUs (PR-MILP) 1024 80.9 1829 150 52.0 213 708 3607 76.7 3529 3616c 3618 1290 59.3 3607 3615 3614 3614c 240 336 3613 160 742 2118 1198 4130c 4775c 4177c 911c 488c 305 676

(PRT-MILP) 4357 453 476 579 525 593 395 4656 449 1323 4448 4214 1123 788 965 4722 760 3421 570 757 4376 1336 1850 13,005 6610 9631c 1762 c 8674c 2510c 1871c 2401 2268

GloMIQO 2.3 3600 3600 3600 24.2 14.1 132 15.8 3600 3600 2989 3600 3282 39.4 3600 3600 3600 3600 3600 3600 49.7 3600 3600 3600 3600 3600 3600c 3600c 3600c 3600 3600 3600 3600

BARON 12.3.3 3600 3600 3600 3600 84.8 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600c 3600 3600 3600d 3600d 3600d 3600c 3600d 3600 3600

Optimality gap calculated with respect to best known solution (optimum) rather than upper bound returned by the algorithms. Solution taken from MILP-based heuristic procedures (Teles et al., 2009). Suboptimal solution returned. No feasible solution found.

computational time to a few hours (almost 11 h for Ex28 with N = 5) without major benefits in optimality gap. In some cases, Ex27, Ex30, Ex32 and Ex33, the branch-and-bound procedure is even less efficient for (PRT-MILP), leading to higher gaps at termination.

6.2. Comparison to commercial solvers The new tightened piecewise relaxation algorithm is now compared to global optimization solvers GloMIQO (Misener and Floudas, 2013b) and BARON (Tawarmalani and Sahinidis, 2005). These are based on spatial branch-and-bound and involve multiple solves of (P) with a local solver following decisions regarding domain partitioning. Thus, better feasible solutions z’ are expected compared to the new algorithm that involves just two solves of (P) (the last one starting from solution zR , which although not feasible in (P) provides a very good initialization point). This is confirmed by the results in Table 6 that show that GloMIQO is the best at finding the optimal solution, returning suboptimal solutions only for Ex31, Ex32 and Ex33. In contrast, BARON fails to find even a feasible solution for Ex30, Ex31, Ex32 and Ex34, while the solutions for Ex27 and Ex33 are suboptimal. Models (PR-MILP) and (PRT-MILP) lead to suboptimal solutions in seven and five cases, respectively. The main focus of this paper is however on the quality of lower bound zR and on showing that better values can be obtained when using the original (Bergamini et al., 2005) or tightened

Fig. 2. Performance profile for optimality gap (water network problems).

piecewise relaxation approach, preceded by the standard McCormick optimization bound contraction. The results in Table 6 show very significant differences in optimality gap, leading to different performance profiles (Dolan and Moré, 2002) in Fig. 2. Optimality gap can be directly related to the quality of the underlying relaxation method since it was computed using the best-known (optimal) solution from the second column as the upper bound. The improved upper bound explains why the gap for (PR-MILP) in Ex31 is now lower than that from (PRT-MILP) (compare to Table 5). Overall, the new tightened approach (PRT-MILP) has the highest probability (53%) of being the best performer (results for  = 0

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model CACE-4938; No. of Pages 12

ARTICLE IN PRESS P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

10

in Fig. 2) and remains above the other solvers up to  = 3.8, meaning that it is the most likely to have a gap at most 23.8 = 14 times higher than the lowest gap of the group. GloMIQO has the second highest probability (34%) of being the best performer but quickly falls below (PR-MILP), catching up at  = 2.6 and actually becoming the best in the interval  ∈ [3.8, 10[. BARON is the worst performer, leading to the lowest gap merely in Ex8 and Ex21. Furthermore, the highest gap for (PRT-MILP) is 11.5%, obtained for Ex32, whereas GloMIQO and BARON terminate three times with gaps over 100%. 6.2.1. Remarks While the quality of the relaxation increases with the number of partitions N, one should be cautious with the extent of piecewise relaxation schemes in order to avoid generating intractable MILP problems (Gounaris et al., 2009). Several works analyze the influence of N in the reduction of optimality gap (see for instance Gounaris et al., 2009, Misener et al., 2011, WittmannHohlbein and Pistikopoulos, 2013, 2014), sometimes identifying the most appropriate partitioning level and scheme for a particular problem. However, to the best of our knowledge, none has proposed a method to estimate the number of partitions as a function of the problem complexity. To provide a fair comparison with the commercial solvers GloMIQO and BARON, which ran with default options, we have used in Table 6 the following simple formula: N =1+

 ˛ |I| · |J|

to compute the number of partitions as a function of the number of non-partitioned variables |I| and partitioned variables |J|. The formula meets three important criteria: (i) returns an integer value; (ii) ensures a minimum of two partitions per problem so as to benefit from the piecewise relaxation scheme; (iii) reduces the number of partitions with the increase in problem size. It was derived after noting that the number of partitions is roughly inversely proportional to the problem size (r2 = 0.76 for a linear regression with zero intercept), when picking the best N values from Tables 4 and 5 for (PRT-MILP). In this case, we have used ˛ = 1.8E4. It should also be noted that the proposed algorithm is not a global optimization algorithm. It is capable of providing a highquality feasible solution and an optimality gap for a given N but cannot in general benefit from augmenting the computational resources. In contrast, the optimality gap for the commercial solvers decreases monotonically with the computational time.

piecewise relaxation schemes and optimality-based bound contraction should be used to a greater extent. Future work will thus look at the development of a global optimization algorithm with such features.

Acknowledgments The author gratefully acknowledges financial support from the Luso-American Foundation under the 2013 Portugal-U.S. Research Networks Program and from Fundac¸ão para a Ciência e Tecnologia through the Investigador FCT 2013 program.

Appendix. Piecewise McCormick relaxation with bivariate partitioning In bivariate partitioning, there are no partition-dependent variable bounds to be tightened since the domain of the two variables of a bilinear term is known a priori for each twodimensional partition. Let binary variable yinjn indicate the active L and xU ) and x (xL partition for variables xi (bounded by xin j in jn U ). The piecewise McCormick relaxation corresponding to and xjn  the bivariate partitioning case can be formulated as the following Generalized Disjunctive Program (PR-BV-GDP). While we consider that the domain of both variables is divided into an equivalent number of regions N, it is straightforward to consider otherwise.

7. Conclusions This paper has shown that the quality of the relaxation resulting from piecewise McCormick envelopes for bilinear terms can be improved significantly once partition-dependent lower and upper bounds are defined for all bilinear variables. Such bounds can easily be found through optimality-based bound contraction, which involves solving a number of linear problems that increases linearly with the number of partitions. The new approach has been implemented in the form of an optimization algorithm that can find a high quality solution to a bilinear problem and compute an optimality gap. Through the solution of test problems from the literature, it has been shown that orders of magnitude reduction in optimality gap can be achieved compared to the conventional piecewise McCormick approach and that the total computational time is sometimes even lower. It has also been shown that the proposed algorithm outperforms state-of-the-art commercial global optimization solvers for the specific class of waterusing network design problems. This is an indication that

The convex hull reformulation is (PR-BV-MILP) and features disaggregated variables with five indices (e.g. xˆ iinjn ).

min z R = f0 (x) =



aij0 wij + h0 (x)

(PR-BV-MILP)

(i,j) ∈ BL

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model

ARTICLE IN PRESS

CACE-4938; No. of Pages 12

P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

11

Table 7 Univariate vs. bivariate partitioning for (P1). Relaxation approach

Partitions (N)

16

144

2500

10,000

(PR-MILP)

Total CPUs Lower bound Optimality gap

0.88 −1.102097 1.70%

1.42 −1.085467 0.197%

9.11 −1.083456 0.0114%

429 −1.083364 0.0029%

N

4

12

50

100

125

(PR-BV-MILP)

Optimality gap Lower bound Total CPUs

1.28% −1.09743 0.93

0.188% −1.085369 0.91

0.0078% −1.083418 7.55

0.0029% −1.083364 217

0.0018% −1.083353 469

(PR-MILP) & (PR-BV-MILP)

Binary variables Total variables Equations

16 52 74

144 436 586

2500 7504 10,010

10,000 30,004 40,010

15,625 46,879 62,510

subject to fq (x) =



aijq wij + hq (x) ≤ 0 ∀q ∈ Q \{0}

(i,j) ∈ BL N N

⎫  ⎪ ⎪ ⎪ wij ≥ ⎪ ⎪ ⎪  ⎪ n=1 n =1 ⎪ ⎪ N N ⎪    ⎪ ⎪ U U U U ⎪ wij ≥ xˆ iinjn · xjn + xin · xˆ jinjn − xin · xjn · yinjn ⎪ ⎪ ⎪ ⎪ n=1 n =1 ⎪ ⎪ N N  ⎪  ⎪ ⎪ U U L L ⎪ wij ≤ xˆ iinjn · xjn + xin · xˆ jinjn − xin · xjn · yinjn ⎪ ⎪ ⎪  ⎪ n=1 n =1 ⎪ N N ⎬  ⎪   

L L L L ˆ jinjn − xin · xjn xˆ iinjn · xjn  + xin · x  · yinjn

U + xL · x L · xU · y ˆ jinjn − xin xˆ iinjn · xjn  injn jn in

wij ≤

n=1 n =1

xi =

N N  

xˆ iinjn

n=1 n =1 N N

xj =



 N

xˆ jinjn

n=1 n =1 N

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

∀(i, j) ∈ BL

Fig. 3. Partitions resulting from univariate and bivariate partitioning schemes for (P1) when considering 16 binary variables to identify the optimal region of the relaxation problem (PR-MILP), on the left, and (PR-BV-MILP), on the right.

faster. Compared to the tightened univariate approach (PRT-MILP), the optimality gap of 0.0018% for 1252 = 15,625 partitions is still one order of magnitude larger, despite the higher computational time (compare to Table 1). References

Al-Khayyal FA, Falk JE. Jointly constrained biconvex programming. Math Oper Res 1983;8(2):273–86. Balas E. Disjunctive programming and a hierarchy of relaxations for discrete optin=1 n =1  U ·y L ·y mization problems. SIAM J Algebraic Discrete Math 1985;6:466–86.



ˆ iinjn ≤ xin xin injn ≤ x injn ML, Aguirre P, Grossmann IE. Logic-based outer approximation ∀ (i, j) ∈ BL, n ∈ 1, . . ., N , n ∈ 1, . .Bergamini ., N U ·y L ·y ˆ xjn ≤ x ≤ x     injn injn jinjn for globally optimal synthesis of process networks. Comput Chem Eng jn ⎫ U − xL ) · (n − 1) 2005;29:1914–33. (x ⎪ i i ⎬ L = xL + Castro PM, Grossmann IE. Optimality-based bound contraction with multipaxin



i N rametric disaggregation for the global optimization of mixed-integer bilinear ∀ i| (i, j) ∈ BL , n ∈ 1, . . ., N (xiU − xiL ) · n problems. J Global Optim 2014., http://dx.doi.org/10.1007/s10898-014-0162-6. ⎪ U = xL + ⎭ xin i Castro PM, Teles JP. Comparison of global optimization algorithms for the design of N ⎫ water-using networks. Comput Chem Eng 2013;52:249–61. (xjU − xjL ) · (n − 1) ⎪ ⎬ L = xL + Catalão JPS, Pousinho HMI, Mendes VMF. Hydro energy systems management in xjn

 j Portugal: profit-based evaluation of a mixed-integer nonlinear approach. Energy N ∀ j| (i, j) ∈ BL , n ∈ {1, . . ., N} (xjU − xjL ) · n 2011;36:500–7. ⎪ ⎭ U L xjn = xj + Dolan ED, Moré JJ. Benchmarking optimization software with performance profiles. N Math Program Ser A 2002;91:201–13. xL ≤ x ≤ xU Rubio-Castro E, Ponce-Ortega JM, Serna-González M, El-Halwagi MM, Pham V. m x∈R Global optimization in property-based inter-plant water integration. AIChE J





2013;59(3):813–33. yinjn ∈ 0, 1 ∀ (i, j) ∈ BL, n ∈ 1, . . ., N , n ∈ 1, . . ., N Faria DC, Bagajewicz MJ. Novel bound contracting procedure for global optimization of bilinear MINLP problems with applications to water management problems. The motivating example (P1), with its single bilinear term, Comput Chem Eng 2011;35:446–55. Faria DC, Bagajewicz MJ. A new approach for global optimization of a class of MINLP provides a good setting for comparing univariate and bivariate problems with applications to water management and pooling problems. AIChE partitioning schemes. Partitions were selected so as to generate J 2012;58(8):2320–35. the same number of binary variables, with Fig. 3 showing the Floudas CA, Pardalos PM, Adjiman CS, Esposito WR, Gumus ZH, Harding ST, Klepeis JL, Meyer CA, Schweiger CA. Handbook of test problems in local and global regions associated to the sixteen different partitions. While the optimization. Boston: Kluwer Academic Publishers; 1999. optimal values from the lower bounding problems are very simGalan B, Grossmann IE. Optimal design of distributed wastewater treatment ilar, (1.183867, 0.551602) for (PR-MILP) and (1.159928, 0.479784) networks. Ind Eng Chem Res 1998;37:4036. Gounaris CE, Misener R, Floudas CA. Computational comparison of piecewise-linear for (PR-BV-MILP), the latter is considerably tighter, leading to an relaxations for pooling problems. Ind Eng Chem Res 2009;48:5742. optimality gap of 1.28% (vs. 1.70%), see Table 7. Bivariate partitionGrossmann IE, Trespalacios F. Systematic modeling of discrete-continuous optiing continues to provide better quality lower bounds up to 10,000 mization models through generalized disjunctive programming. AIChE J 2013;59:3276–95. partitions, when the values become the same. It is also generally

yinjn = 1

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025

G Model CACE-4938; No. of Pages 12 12

ARTICLE IN PRESS P.M. Castro / Computers and Chemical Engineering xxx (2014) xxx–xxx

Harjunkoski I, Westerlund T, Pörn R, Skrifvars H. Different transformations for solving non-convex trim loss problems by MINLP. Eur J Oper Res 1998;105: 594–603. Hasan MMF, Karimi IA. Piecewise linear relaxation of bilinear programs using bivariate partitioning. AIChE J 2010;56:1880. Hock W, Schittkowski K. Test examples for nonlinear programming codes, vol. 187 of Lecture Notes in Economics and Mathematical Systems. Berlin Heidelberg: Springer Verlag; 1981. Jia Z, Ierapetritou M. Mixed-integer linear programming model for gasoline blending and distribution scheduling. Ind Eng Chem Res 2003;42:825–35. Jia Z, Ierapetritou M, Kelly JD. Refinery short-term scheduling using continuous time formulation: crude-oil operations. Ind Eng Chem Res 2003;42:3085–97. Karuppiah R, Grossmann IE. Global optimization for the synthesis of integrated water systems in chemical processes. Comput Chem Eng 2006;30:650–73. Kolodziej S, Castro PM, Grossmann IE. Global optimization of bilinear programs with a multiparametric disaggregation technique. J Global Optim 2013a;57:1039–63. Kolodziej S, Grossmann IE, Furman KC, Sawaya NW. A discretization-based approach for the optimization of the multiperiod blend scheduling problem. Comput Chem Eng 2013b;53:122–42. Lee H, Pinto JM, Grossmann IE, Park S. Mixed-integer linear programming model for refinery short-term scheduling of crude oil unloading with inventory management. Ind Eng Chem Res 1996;35:1630–41. McCormick GP. Computability of global solutions to factorable nonconvex programs. Part I. Convex underestimating problems. Math Progr 1976;10:146. Meyer CA, Floudas CA. Global optimization of a combinatorially complex generalized pooling problem. AIChE J 2006;52(3):1027. Misener R, Floudas CA. Global optimization of mixed-integer quadraticallyconstrained quadratic programs (MIQCQP) through piecewise-linear and edge-concave relaxations. Math Prog Ser B 2012;136(1):155–82. Misener R, Floudas CA. Mixed-Integer Quadratically-Constrained Quadratic Programs: GloMIQO 2.2 Test Suite; 2013a, http://helios.princeton.edu/GloMIQO/ MisenerFloudas GloMIQO TestSet.pdf Misener R, Floudas CA. GloMIQO: global mixed-integer quadratic optimizer. J Global Optim 2013b;57:3–50. Misener R, Thompson JP, Floudas CA. APOGEE: global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes. Comput Chem Eng 2011;35:876. Moro LFL, Zanin AC, Pinto JM. A planning model for refinery diesel production. Comput Chem Eng 1998;22:S1039–42.

Nápoles-Rivera F, Ponce-Ortega JM, El-Halwagi MM, Jiménez-Gutiérrez A. Global optimization of mass and property integration networks with in-plant property interceptors. Chem Eng Sci 2010;65:4363. Quesada I, Grossmann IE. A global optimization algorithm for linear fractional and bilinear programs. J Global Optim 1995;6:39–76. Raman R, Grossmann IE. Modeling and computational techniques for logic based integer programming. Comput Chem Eng 1994;18:563–78. Rodriguez MA, Vecchietti A. A comparative assessment of linearization methods for bilinear models. Comput Chem Eng 2013;48:218–33. Ruiz JP, Grossmann IE. Exploiting vector space properties to strengthen the relaxation of bilinear programs arising in the global optimization of process networks. Optim Lett 2011;5:1–11. Tawarmalani M, Sahinidis NV. A polyhedral branch-and-cut approach to global optimization. Math Program 2005;103(2):225–49. Teles JP, Castro PM, Matos HA. Global optimization of water networks design using multiparametric disaggregation. Comput Chem Eng 2012;40:132–47. Teles JP, Castro PM, Matos HA. Multi-parametric disaggregation technique for global optimization of polynomial programming problems. J Global Optim 2013a;55:227–51. Teles JP, Castro PM, Matos HA. Univariate parameterization for global optimization of mixed-integer polynomial problems. Eur J Oper Res 2013b;229: 613–25. Teles JP, Castro PM, Novais AQ. MILP-based initialization strategies for the optimal design of water-using networks. Chem Eng Sci 2009;64:3736–52. Tsai J-F, Lin M-H. An improved framework for Solving NLIPs with signomial terms in the objective or constraints to global optimality. Comput Chem Eng 2013;53:44–54. Wicaksono DN, Karimi IA. Piecewise MILP under- and overestimators for global optimization of bilinear programs. AIChE J 2008;54:991. Wittmann-Hohlbein M, Pistikopoulos EN. On the global solution of multi-parametric mixed integr linear programming problems. J Global Optim 2013;57:51–73. Wittmann-Hohlbein M, Pistikopoulos EN. Approximate solution of mp-MILP problems using piecewise affine relaxation of bilinear terms. Comput Chem Eng 2014;61:136–55. Yadav S, Shaik MA. Short-term scheduling of refinery crude oil operations. Ind Eng Chem Res 2012;51:9287–99. Zorn K, Sahinidis NV. Computational experience with applications of bilinear cutting planes. Ind Eng Chem Res 2013;52:7514–25.

Please cite this article in press as: Castro PM. Tightening piecewise McCormick relaxations for bilinear problems. Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.03.025