Computational Experience with Applications of Bilinear Cutting Planes

Article pubs.acs.org/IECR

Computational Experience with Applications of Bilinear Cutting Planes Keith Zorn and Nikolaos V. Sahinidis* Department of Chemical Engineering, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, Pennsylvania 15228, United States ABSTRACT: This work analyzes the effectiveness of cutting planes applied during the global optimization of nonconvex nonlinear and mixed-integer nonlinear programming problems arising from process synthesis and process operations. In a previous publication, we introduced algorithms for bilinear substructure detection, cutting plane identification, cut filtering, and cut selection. Our implementation embeds reformulation−linearization, lift-and-project, and advanced convex envelope construction techniques that exploit underlying bilinear substructures and tighten factorable programming reformulations at every node in the branch-and-bound tree. In this work, we utilize industrial examples to demonstrate a variety of relaxation tightening strategies afforded by these bilinear cutting plane algorithms. A detailed description and associated computational experience are presented for applied cutting planes. substructure detection, cutting plane identification, cut filtering, and cut selection associated with the exploitation of underlying bilinear substructures in generalized nonconvex programs. The computational results showed that these practical and efficient algorithms detect underlying bilinear substructures and tighten relaxations during the solution of more than 30% of the problems in the standard literature test libraries GLOBALLib9 and MINLPLib.10 In this work, we analyze the effectiveness of cutting planes applied during the solution of nonconvex nonlinear and mixed-integer nonlinear programs from process synthesis and process operations. The remainder of this article is organized as follows: In section 2, we review classes of bilinear cutting planes generated by embedded reformulation−linearization, lift-and-project, and advanced convex envelope generation techniques. Detailed descriptions of applied cutting planes and associated computational experiences for process synthesis and process operations applications are presented in sections 3−5. In section 6, we offer a discussion and conclusions about the effectiveness of applied cutting planes.

1. INTRODUCTION In the fields of engineering and operational research, mathematical models often involve nonlinear functional forms, discrete decisions, coupled interactions, and/or stochastic parametric characterizations. In recent years, industrial competition has necessitated the solution of increasingly sophisticated process and network models to maximize profit. In turn, the increased complexity of industrial nonlinear and mixed-integer nonlinear models and the desire to locate provably superior extrema has generated the need for advanced global optimization software. The result is an increasingly rich global optimization literature. Horst and Pardalos1 compiled a global optimization handbook containing chapters written by relevant experts covering many aspects of global optimization. Other published books address important global optimization topics including interval arithmetic,2 constraint propagation,3 and convex analysis.4,5 This work concerns an implementation that incorporates these important aspects into a branch-and-bound algorithm for the global optimization of industrially relevant nonlinear programs (NLPs) and mixed-integer nonlinear programs (MINLPs). The Branch-and-Reduce Optimization Navigator (BARON)6 is a state-of-the-art, general-purpose global optimization package that leverages polyhedral relaxations and outer approximations for bounding and range reduction. Polyhedral relaxations for nonconvex programs are generated by sequentially replacing factorable terms by introduced linearization variables.6,7 During factorable reformulation, each complicating factorable term is decomposed into linear, logarithmic, exponential, monomial, and bilinear subexpressions.5 The solution of each resulting linear relaxation provides a valid bound on the optimal objective. Consequently, the relaxation strength has a significant impact on BARON’s rate of convergence. BARON’s standard factorable reformulations are enhanced with cutting planes at every node in the branch-and-bound tree. In a previous work,8 we introduced algorithms for bilinear © XXXX American Chemical Society

2. CUTTING PLANE GENERATION In this section, we review classes of linear cutting planes that exploit underlying bilinear substructures. For illustration, consider the following factorable nonlinear program P Received: December 7, 2012 Revised: May 2, 2013 Accepted: May 7, 2013

A

dx.doi.org/10.1021/ie3033763 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research min s.t.

Article η

f0 (x) fi (x) ≤ 0

∑ aljwjk − blxk ≤ 0 i = 1, ..., q

(7) η

Ax ≤ b xj ∈ 

j = 1, ..., s

xj ∈ 

j = s + 1, ..., n

blxk −

∑ aljwjk ≤ 0

l = r + 1, ..., m ; xk ≤ 0 (8)

j=1

where each bilinear term wjk = xjxk is subsequently replaced by the introduced variable xj, j ∈ {n + 1, ..., η}. Because these cuts do not depend on variable bounding information, they are classified as static cuts. Static cuts are added at the root node during preprocessing and are valid for the entire branch-andbound tree. The lifted cuts in eqs 6−8 enhance range reduction and tighten relaxations by revealing relationships between existing columns xk, k ∈ {1, ..., η}, and bilinear columns xj, j ∈ {n + 1, ..., η}. When no sign restriction is available on the lifting variable xk, k ∈ {1, ..., n}, the lifted equations

where all linear constraints are embedded in Ax ≤ b, where A ∈ Rm×n; b ∈ Rm; and f i(x), i ∈ Q = {0, 1, ..., q}, represents a set of nonlinear factorable functions. Both linear and nonlinear constraints are composed of any combination of integer and continuous variables. After each factorable function f i(x) has been sequentially decomposed, a factorable programming approach introduces linearization variables xj, j ∈ {n + 1, ..., η}, and includes the following set of linear relationships in a polyhedral relaxation of program P η

∑ aljxj = bl

l = r + 1, ..., m ; xk ≥ 0

j=1

n

l = 1, ..., r

(xk − xkL) ∑ aljxj ≤ bl(xk − xkL)

(1)

j=1

l = r + 1, ..., m

j=1

(9)

η

∑ aljxj ≤ bl

l = r + 1, ..., m

n

(xkU − xk) ∑ aljxj ≤ bl(xkU − xk)

(2)

j=1

where η is the cardinality of the column space after factorable relaxation.5 In the following subsections, we detail cutting planes based on these linear relationships. These new classes of cutting planes reveal relationships between existing and intermediate linear and bilinear substructures encountered during factorable reformulation. These new relationships effectively explore range reduction, based on interval operations, and bounding, based on tightening polyhedral relaxations. 2.1. Lifted Cutting Planes. Lifted cutting planes for bilinear structures are generated by multiplying existing linear relationships with a lifting variable xk, k ∈ {1, ..., η}. These cutting planes are based on the reformulation−linearization techniques of Sherali and co-workers.11−14 In this section, we review classes of static and dynamic lifted cutting planes. In the case of the linear equality relationships in eq 1 and the linear inequality relationships in eq 2 with a sign restriction on the lifting variable xk, k ∈ {1, ..., n}, the lifted equations

(10)

are generated by applying reformulation−linearization techniques to eq 2. After extension to the complete factorable relaxation space η, algebraic expansion, and sequential relaxation, the proposed algorithm identifies the following cutting planes inspired by eqs 9 and 10 for k ∈ {1, ..., η} η

j=1

η

j=1

l = r + 1, ..., m ; xk ≥ 0

n

l = r + 1, ..., m ; xk ≤ 0 (5)

j=1

are generated by applying reformulation−linearization techniques. After extension to the complete factorable relaxation space η, algebraic expansion, and sequential relaxation, our convexification approach for the global optimization of problems with bilinear intermediates identifies the following cutting planes inspired by eqs 3−5 for k ∈ {1, ..., η} η

∑ aljwjk − blxk = 0 j=1

j=1

(12)

where, again, each bilinear term wjk = xjxk is subsequently replaced by the introduced variable xj, j ∈ {n + 1, ..., η}, and xLk and xUk represent the lower and upper bounds, respectively, on variable xk. Because these cutting planes require variable bounding information, they are classified as dynamic cuts and are updated after each bounding and range reduction operation at every node in the branch-and-bound tree. The lifted cuts in eqs 11 and 12 reveal relationships between existing columns xj and xk, j,k ∈ {1, ..., η}, and bilinear columns xj, j ∈ {n + 1, ..., η}, that aid in bounding and range reduction. More information about lifted cutting planes and the algorithms used by the proposed implementation to form these relationships is available in ref 8. 2.2. Disaggregation Cutting Planes. In this section, we review dynamic disaggregation cutting plane generation. In the case that a variable xj, j ∈ {n + 1, ..., η}, in one of the linear equality relationships in eq 1 or linear inequality relationships in eq 2 arises from sequential replacement of a factorable bilinear term, the proposed approach scans for a factored structure. For illustration, consider the bilinear term wjk composed of a variables xj and xk, that is

(4)

blxk ≤ xk ∑ aljxj

∑ aljxkUxj + blxk ≤ blxkU

l = r + 1, ..., m

(3)

j=1

(11)

η

−∑ aljwjk +

n

xk ∑ aljxj ≤ blxk

j=1

l = r + 1, ..., m

l = 1, ..., r

j=1

η

∑ aljwjk − ∑ aljxkLxj − blxk ≤ −blxkL

n

xk ∑ aljxj = blxk

l = r + 1, ..., m

j=1

l = 1, ..., r (6) B


Industrial & Engineering Chemistry Research

Article

Figure 1. Description of the trim loss problem.

wjk = xjxk

j , k ∈ {1, ..., η}

η

(13)

wjk ≥

i=1

If at least one of the constituent variables, for example, xk, arises from a factorable linear identity

− xiUxjU}

η

xk =

∑ aixi

wjk ≤

(14)

− xiUxjL}

j , k ∈ {1, ..., η} (15)

Alternatively, each term in eq 15 can be distributed and sequentially replaced by the introduced variable wij η

∑ aiwij

j , k ∈ {1, ..., η} (16)

i=1

during factorable reformulation. Quesada and Grossmann15 introduced a linear variable xk as in eq 14 and used the resulting formulation to bound the original bilinear term wjk. Kearfott2 described the effects of a factored formulation such as that in eq 15 on the range reduction of the final linearization variables wjk compared to the disaggregated form in terms of variable dependence. Interval arithmetic performed on the disaggregated form neglects the coupled dependence of each individual term on the common component xj because xj appears separately in the convex and concave envelopes of every introduced term wij. Consequently, range reduction is enhanced for the factored form because it inherently preserves this variable dependence. Tawarmalani and Sahinidis5 and Tawarmalani et al.,16 however, demonstrated that the disaggregated formulation produces a sharper relaxation. To promote enhanced range reduction, the proposed methodology maintains the factored formulation. To leverage the relaxation strength of the disaggregated formulation, an additional relationship is generated between the two formulations

3. TRIM LOSS PROBLEM In this section, we detail the application of bilinear cutting planes to trim loss optimization. Among the first formulations published in the context of linear programming,19 cutting and packing problems span many fields of literature including: industrial engineering, operational research, computer science, and production engineering. Since the 1960s, many variations of cutting stock problems have been identified and classified in the literature.20−27 In the following subsections, we detail a specific formulation associated with the minimization of trim loss, detail cutting planes that exploit an underlying bilinear substructure, and provide computational results demonstrating enhanced branch-and-bound performance. 3.1. Problem Statement. The trim loss problem is a production planning problem that involves cutting rolls of raw material into smaller rolls of product. During the cutting process, blades are oriented in fixed patterns to produce product rolls of specified widths. The objective is to satisfy product orders while minimizing production costs. In practice, these challenging mixed-integer nonlinear programming problems are often solved heuristically and/or after mathematical transformation because standard convexification of the original MINLPs generates large factorable relaxation gaps.

η

wjk =

∑ aiwij

j , k ∈ {1, ..., η}

i=1

(19)

for j,k ∈ {1, ..., η}, is added to the standard factorable relaxation after wjk has been replaced by the introduced variable xj, j ∈ {n + 1, ..., η}. Because these relationships require variable bounding information, they are classified as dynamic cuts and are updated after each branching and range reduction operation at every node in the branch-and-bound tree. Disaggregation cuts tighten standard relaxations by revealing relationships between columns xi, xj, j ∈ {1, ..., η}, and bilinear column xj, j ∈ {n + 1, ..., η}. More information about dynamic disaggregation cutting planes and how to form these relationships is available in ref 8.

η

i=1

∑ ai min{xixjU + xiLxj − xiLxjU , xixjL + xiUxj i=1

then wjk contains a factored structure

wjk =

(18)

η

k ∈ {1, ..., η}

i=1

wjk = xj ∑ aixi

∑ ai max{xixjL + xiLxj − xiLxjL , xixjU + xiUxj

(17)

A polyhedral outer approximation for eq 17, given by17,18 C



Article I

Figure 1 describes a cutting process associated with a trim loss problem. The width of a roll of raw stock material is labeled Bmax, and each variable bi represents the specified width of a product i ∈ 1, ..., I. Variables nij indicate the number of times product i is produced by each cutting pattern j ∈ 1, ..., J. Finally, mj is the number of times that each cutting pattern is repeated during the process. Shaded areas in Figure 1 indicate waste, or “trim loss”. Harjunkoski et al.28 detailed the following MINLP trim loss formulation that arises naturally in a paper production facility

∑ wij − Nmaxmj ≤ 0

(20a)

j=1 I

s.t.

∑ binij − Bmax ≤ 0

j = 1, ..., J (20b)

i=1 I

Bmax −

∑ binij − Δ ≤ 0

j = 1, ..., J (20c)

i=1

I

∑ nij − Nmax ≤ 0

j = 1, ..., J

yj − mj ≤ 0

Table 1. Instance Statistics for Trim Loss Problems

(20d)

i=1

j = 1, ..., J

mj − Myj ≤ 0

(20e)

j = 1, ..., J

(20f)

J

ni ,order −

∑ mjnij ≤ 0

i = 1, ..., I (20g)

j=1

nij , mj ∈ 

+

yj ∈ {0, 1}

i = 1, ..., I ; j = 1, ..., J

j = 1, ..., J

(20h) (20i)

j = 1, ..., J

∑ biwij ≤ 0 i=1

no. of variables

tln2 tln4 tln5 tln6 tln7

13 25 31 37 43

9 25 36 49 64

8 24 35 48 63

no. of nonzeroes (NZ)

no. of nonlinear nonzeroes (NNZ)

33 105 156 217 288

8 32 50 72 98

model name

Mstat

tbil (s)

tnobil (s)

tnobil/tbil

tln2 tln4 tln5 tln6 tln7

6 12 15 18 21

0.1 1.0 2.8 4.8 >500.0

0.1 6.3 >500.0 >500.0 >500.0

1.0 6.4 >178.0 >104.0 −

generated using a GAMS 23.8.229 interface to a development version of BARON 11.56,30 on a 64-bit Intel Xeon X5650 2.66 Ghz processor. With the exception of enabling and disabling the bilinear cutting plane algorithms and setting the absolute and relative termination tolerances to 10−6, all GAMS/BARON options were set to their default values. In Table 2, Mstat indicates the number of static cutting planes in eqs 21−23 applied at the root node during preprocessing. tbil and tnobil represent the times required to solve each instance to global optimality with and without the use of bilinear cutting planes, respectively. The final column of Table 2 contains a relative comparison of tnobil and tbil signifying an acceleration achieved through the application of bilinear cutting planes. On the smallest instance, tln2, the effect of bilinear cutting planes is obfuscated by a rapid solution time. tln4, however, exhibits an acceleration factor of greater than 6 after the application of 12 static cutting planes. BARON is unable to solve the midsized instances, tln5 and tln6, to global optimality

(21)

I

(Bmax − Δ)mj −

no. of equations

Table 2. Computational Experience for Trim Loss Problems

I i=1

model name

no. of discrete variables (D)

Table 2 describes our computational experience while solving these instances to global optimality. All reported results were

Here, binary variables yj indicate the activity of each cutting pattern in the production process. The objective in eq 20a minimizes the individual costs of raw materials for each individual cutting pattern cj and a setup cost for each cutting pattern Cj. The constraints in eqs 20b and 20c ensure that each cutting pattern fits on the roll of raw material within a specified distance tolerance Δ, and the number of product repeats in a single cutting pattern is limited by eq 20d. Inactive cutting patterns are disabled by the constraints in eqs 20e and 20f, and eq 20g ensures that all product demands ni,order are satisfied. 3.2. Cutting Planes. Standard factorable relaxation of the trim loss MINLP formulation results in a weak linear relaxation. To tighten relaxations and aid in range reduction, static lifted cutting planes are added. Here, we detail cuts added by the proposed methodology to exploit the bilinear structure of the original problem. The original demand constraint in eq 20g contains bilinear terms mjnij. Standard factorable reformulation introduces the variable wij = mjnij for each bilinear combination. The linear relationships

∑ biwij − Bmax mj ≤ 0

(23)

are generated by lifting the linear constraints in eqs 20b−20d with the number of pattern repeats mj. These new cutting planes in eqs 21−23 represent lifted material and cutting pattern usage limits that reveal relationships between the existing bilinear combinations of pattern repeats and product repeats in the constraint in eq 20g and the number of cutting pattern repeats. More about the filtering and application of lifting cutting planes is available in ref 8. 3.3. Computational Results. To assess the effectiveness of the cutting planes applied by our proposed implementation during the solution of trim loss problem formulations, we include computational results for trim loss instances from MINLPLib.10 These specific problems have been shown to suffer from weak linear programming relaxations and are typically solved only after transformation.28 Model statistics, including the numbers of equations, variables, discrete variables (D), nonzeroes (NZ), and nonlinear nonzeroes (NNZ), are reported in Table 1.

J

min ∑ (cjmj + Cjyj )

j = 1, ..., J

i=1

j = 1, ..., J (22) D



Article

factorable programming relaxation “nobil”. Tightened relaxations result in sharper bounds and enhanced range reduction and are responsible for the improved computational performance demonstrated in Table 2.

within 500 s without the use of bilinear cutting planes. After the application of a modest number of cutting planes, these problems are solved in only a few seconds. For the largest and most challenging instance, tln7, BARON is unable to converge to a global optimum in a time limit of 500 s even with the aid of bilinear cutting planes. For all problems, as the number of products and the number of cutting patterns increase, an increasing number of relaxation tightening static lifting cuts, representing lifted material and usage limits, are applied. For the trim loss application, the number of integer combinations increases with problem size. As combinatorics begin to dominate the degree of difficulty imposed by nonlinearity, increased computational resources are required to certify global optimality. tln7 is an excellent example of this, and a combination of integer programming preprocessing and bilinear cutting planes might be necessary to enable BARON to solve larger trim loss instances. Static cutting planes applied by the proposed bilinear cutting plane algorithms are applied during preprocessing and aid in bounding and range reduction throughout the branch-andbound tree. The effect on required computational resources is demonstrated in Table 2. Because root node optimality gaps are obtained before beginning the branch-and-bound procedure, they are useful metrics for assessing relaxation fitness independently from partitioning and branching decisions, bound propagation and range reduction, and other algorithmic effects. To examine the effect of lifted static cutting planes on the linear relaxations of the trim loss problem, optimality gaps at the root node before and after the application of bilinear cutting planes are compared in Figure 2.

4. HEAT EXCHANGER NETWORK SYNTHESIS In this section, we investigate the application of bilinear cutting planes to heat exchanger network synthesis problems. Heat integration is one of the oldest and most studied process synthesis topics and is the focus of a rich literature including several books31−34 and review articles.35−38 In the following subsections, we examine a simultaneous optimization model for heat integration, detail cutting planes applied during branchand-bound optimization, and present computational results demonstrating the effectiveness of cutting plane application. 4.1. Problem Statement. Heat integration is a design strategy that considers energy flows throughout an entire process instead of optimizing single process subunits individually. The result is a reduction of wasted energy and increased process efficiency. As the price of energy increases, process integration results in increased cost savings. Figure 3 depicts a two-stage heat integration superstructure. Simultaneous optimization of this superstructure selects and sizes heat exchangers and connections between process streams to minimize fixed and operational costs. Mixed-integer nonlinear models representing these superstructures have been studied extensively.39−41 Yee and Grossman40 detailed a widely cited formulation with the objective of minimizing the sum of fixed and unit costs I

J

min c unit(∑ ∑

K−1

I

J

J

U ∑ zijk + ∑ ziCu + ∑ zjHu) + ∑ qjHuhcost

i=1 j=1 k=1

i=1

j=1

j=1

I

+

U + acoeff ∑ qiCuccost i=1

⎫aexp ⎧ ⎛1 1 ⎞ ⎪ ⎪ ⎜ H + C ⎟qijk h h ⎪ ⎪ ⎝ i j ⎠ ⎬ + ε ∑∑∑⎨ ⎤1/3 ⎪ i = 1 j = 1 k = 1 ⎪ ⎡ ΔTijk ΔTijk + 1(ΔTijk + ΔTijk + 1) + ε + ε ⎪ ⎪ ⎣⎢ 2 ⎦⎥ ⎭ ⎩ I

J

K

H + ucoeff

⎫aexp ⎧ ⎛1 1 ⎞ Hu ⎪ ⎪ ⎜ C + Hu ⎟qj J ⎪ ⎪ h ⎠ ⎝ hj + ε⎬ ∑⎨ 1/3 ⎤ ⎪ ⎛ T Huin − T Cout + ΔT jHu ⎞ j=1 ⎪ ⎡ j ⎟ + ε⎥ ⎪ ⎪ ⎢(T Huin − T jCout)ΔT jHu⎜ 2 ⎝ ⎠ ⎣ ⎦ ⎩ ⎭ ⎧ ⎫aexp ⎛1 1 ⎞ Cu ⎜ ⎟q ⎪ ⎪ + H Cu ⎪ ⎪ ⎝h h ⎠ i ∑ ⎨ Hout Cuin Cui Hout Cuin Cu + ε⎬ 1/3 ⎤ ⎪ )ΔTi (Ti −T −T + ΔTi ) i = 1 ⎪ ⎡ (Ti + ε⎥ ⎪ ⎢⎣ ⎪ 2 ⎦ ⎩ ⎭ I

C + ucoeff

Figure 2. Root node gap comparison for trim loss problems.

(25a)

Figure 2 presents percentage root node optimality gaps ⎛ obj* − LBroot ⎞ gap = ⎜ ⎟ × 100 |obj*| ⎝ ⎠

Hu Here, binary variables zijk, zCu i , and zj represent the activity of an exchange between hot stream i and cold stream j at location k, hot stream i and the cold utility, and cold stream j and the Hu hot utility, respectively. Similarly, variables qijk, qCu i , and qj are the energies exchanged by streams i and j at location k, hot stream i and the cold utility, and cold stream j and the hot Hu utility, respectively. ΔTijk, ΔTCu i , and ΔTj signify the approach temperatures between streams i and j at location k, hot stream i and the cold utility, and cold stream j and the hot utility, respectively. Parameter cunit is the fixed cost of an exchanger, whereas hUcost and cUcost are costs associated with use of the heating and cooling utilities, respectively. Similarly, acoeff is the

(24)

where obj* represents the globally optimal solution to the original nonlinear problem and LBroot represents optimal objective function of the linear programming relaxation at the root node. For instance tln7, the incumbent solution is used for obj* because the problem does not converge to global optimality within the allotted time. For all instances, the optimality gap after the application of static lifted cuts “bil” is considerably reduced compared to a corresponding standard E



Article

Figure 3. Description of the heat exchanger network synthesis problem.

area cost coefficient for heat exchangers, whereas uHcoeff and uCcoeff are area cost coefficients for heaters and coolers, respectively. Parameters hCi and hHj are stream-individual film coefficients for hot stream i and cold stream j, respectively, and hHu and hCu are stream-individual film coefficients for the hot and cold utilities, respectively. THout and TCout are the target temperatures of hot i j stream i and cold stream j, respectively, and THuin and TCuin are the inlet temperatures of the hot and cold utilities, respectively. Finally, aexp is the cost exponent for heat exchangers, and ε is a small positive number. The constraints J

(TiHin − TiHout)fiH =

balance the energy exchanged by hot stream i and cold stream j, respectively, in stage k. Variables THik and TCjk represent the temperature of hot stream i as it enters stage k and the temperature of cold stream j as it leaves stage k, respectively. Similarly, the equations

i ∈ {1, ..., I }

j=1 k=1

(25b) I

(T jCout − T jCin)f jC =

K−1

∑ ∑ qijk + qjHu

j ∈ {1, ..., J } (25c)

balance the total energy exchanged by hot stream i and cold stream j, respectively. Variables THin and TCin are the supply i j temperatures of hot stream i and cold stream j, respectively. The heat capacity flow rate of hot stream i and cold stream j are given by f Hi and f Cj , respectively. The equations J

−

TikH+ 1)

=

∑ qijk (25d)

∑ qijk

(25g)

TiHin = TiH1

i ∈ {1, ..., I }

(25h)

C T jCin = T jK

j ∈ {1, ..., J }

(25i)

TikH ≥ TikH+ 1

i ∈ {1, ..., I }, k ∈ {1, ..., K − 1}

(25j)

T jkC ≥ T jkC + 1

j ∈ {1, ..., J }, k ∈ {1, ..., K − 1}

(25k)

TiKH ≥ TiHout

i ∈ {1, ..., I }

(25l)

C TiCout ≥ T jK

j ∈ {1, ..., J }

(25m)

i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K − 1}

i=1

j ∈ {1, ..., J }, k ∈ {1, ..., K − 1}

j ∈ {1, ..., J }

qijk − min(eiHc , ejCc)zijk ≤ 0

I

f jC (T jkH − T jkC + 1) =

f jC (T jCout − T jC1) = qjHu

ensure monotonicity. Logical restrictions are provided in the constraints

j=1

i ∈ {1, ..., I }, k ∈ {1, ..., K − 1}

(25f)

set the supply temperatures of hot streams and cold streams, respectively, whereas the constraints

i=1 k=1

fiH (TikH

i ∈ {1, ..., I }

balance the energy exchanged between hot stream i and the cold utility and between cold stream j and the hot utility, respectively. The constraints

K−1

∑ ∑ qijk + qiCu

fiH (TiKH − TiHout) = qiCu

qiCu − eiHcziCu ≤ 0

(25e) F

i ∈ {1, ..., I }

(25n) (25o)


Industrial & Engineering Chemistry Research qjHu − ejCcz jHu ≤ 0

j ∈ {1, ..., J }

Article L L ijk U t3ijk ≥ max{t1ijk Δtijk + t1ijk LΔtijk − t1ijk LΔtijk ,t1 Δtijk + t1ijk UΔtijk

(25p)

U L ijk L − t1ijk UΔtijk } + max{t1ijk Δtijk + 1 + t1 Δtijk + 1

Cu Hc Cc for qijk, qHu i , and qj . Here, ei and ej are the heat contents of hot stream i and cold stream j, respectively. The additional logical constraints

L U U ijk ijk U ijk U − t1ijk LΔtijk + 1 , t1 Δtijk + 1 + t1 Δt1jk + 1 − t1 Δtijk + 1} U U ijk L t3ijk ≤ min{t1ijk Δtijk + t1ijk LΔtijk − t1ijk LΔtijk ,t1 Δtijk + t1ijk UΔtijk

Hu ΔTijk ≤ TikH − T jkC + γij(1 − zijk)

i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K − 1} ΔTijk + 1 ≤

H Tijk +1

−

T jkC + 1

L U ijk L − t1ijk UΔtijk } + min{t1ijk Δtijk + 1 + t1 Δtijk + 1

(25q)

U L L ijk ijk U ijk U − t1ijk LΔtijk + 1 ,t1 Δtijk + 1 + t1 Δtijk + 1 − t1 Δtijk + 1}

(25r)

are applied at the hot and cold ends. Here, γij represents an upper bound on the driving force for energy exchange between streams i and j. Finally, the logical constraints ΔT jHu

T jC1

j ∈ {1, ..., J }

(25s)

ΔTiCu ≤ TikH − T Cuout

i ∈ {1, ..., I }

(25t)

≤T

Huout

−

t4ijk =

(34)

t5ijk = (t4ijk)1/3 i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

Cu

(35)

t6ijk = t5ijk + ε i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

ΔTijk ≥ Tmapp i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

1 ijk t3 + ε 2

i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

are added to restrict ΔT and ΔT , respectively. The Tmapp parameter restricts approach temperatures between streams. The additional constraints Hu

(33)

for i ∈ {1, ..., I}, j ∈ {1, ..., J}, and k ∈ {1, ..., K} are identified as potential cutting planes after disaggregation of eq 31. The additional linear variables

+ γij(1 − zijk)

i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K − 1}

(32)

t 7ijk =

(26)

qijk t6ijk

(36)

i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K } (37)

ΔT jHu ≥ Tmapp

j ∈ {1, ..., J }

(27)

ΔTiCu ≥ Tmapp

i ∈ {1, ..., I }

(28)

are introduced into the objective function as factorable reformulation continues. The bilinear relationship in eq 37 exhibits a factored form. The linear relationships qijk ≥ t 7ijkε + max{t 7ijkt5ijk L + t 7ijk Lt5ijk − t 7ijk Lt5ijk L , t 7ijkt5ijk U + t 7ijk Ut5ijk

impose limits on the temperature differences between hot and cold streams i and j, cold stream j and the hot utility, and hot stream i and the cold utility, respectively. In addition to these restrictions, scalar upper and lower bounds are set for all problem variables. 4.2. Cutting Planes. The heat exchanger network synthesis formulation produces a weak standard factorable relaxation. Dynamic disaggregation cuts enhance bounding and range reduction during branch-and-bound optimization. Here, we analyze cuts selected by the proposed methodology to exploit the factored bilinear structure of the original problem. The linearization variables

− t 7ijk Ut5ijk U}

qijk ≤ t 7ijkε + min{t 7ijkt5ijk U + t 7ijk Lt5ijk − t 7ijk Lt5ijk U , t 7ijkt5ijk L + t 7ijk Ut5ijk − t 7ijk Ut5ijk L}

(29)

t 2ijk = ΔTijk + ΔTijk + 1 i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

(30)

t3ijk = t1ijkt 2ijk i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

(39)

for i ∈ {1, ..., I}, j ∈ {1, ..., J}, and k ∈ {1, ..., K} are identified by disaggregating eq 37. In addition to eqs 32, 33, 38, and 39, four similar families of cutting planes are revealed after disaggregation of the remaining factored bilinear and fractional forms in the objective function. Throughout the branch-and-bound algorithm, rounds of cuts are tested, selected, and applied dynamically to exploit an underlying factored structure and tighten linear programming relaxations at every node. More information about the filtering, selection, and application of dynamic cutting planes is available in ref 8. 4.3. Computational Results. Here, we examine the effects of dynamic disaggregation cutting planes applied by our implementation during the solution of heat exchanger network synthesis instances. The tested problem formulation, synheat, was developed by Yee and Grossmann40 and is available in MINLPLib.10 All instances were solved in the same test environment a sdescribed in section 3.3, and with the exception of enabling and disabling the bilinear cutting plane algorithms and setting the absolute and relative termination tolerances to 10−6, all GAMS/BARON options were set to their default values.

t1ijk = ΔTijk ΔTijk + 1 i ∈ {1, ..., I }, j ∈ {1, ..., J }, k ∈ {1, ..., K }

(38)

(31)

are introduced into the objective function in eq 25a and added to the sequential factorable relaxation. The intermediate relationship in eq 31 involves a factored functional form similar to that of 15. The linear relationships G



Article

5. SUPPLY CHAIN AND INVENTORY MANAGEMENT In this section, we analyze the application of bilinear cutting planes to simultaneous supply chain design and stochastic inventory management. Although supply chain management and inventory management are often treated separately,42−46 a simultaneous optimization approach has been shown to yield considerable benefits.47,48 Recently, stochastic inventory models have been proposed to deal with demand uncertainty at the enterprise level.49−52 In the following subsections, we assess the effectiveness of bilinear cutting planes applied to an MINLP model for large-scale supply chain design with stochastic inventory management. 5.1. Problem Statement. Inventory management decisions are crucial when designing and operating a supply chain. After production or purchase at suppliers, goods are transported and stored at distribution centers before being transported and sold at retailers. The supply chain must be designed to optimize a balance between transportation costs and inventory storage costs while ensuring that a sufficient amount of product is delivered to meet retailer demands. Figure 4 depicts a three-echelon supply chain network structure incorporating J distribution centers and I service retailers. In practice, precise retailer demands are not available during the design of a supply chain network. Often, product demands are modeled as scenarios or selected from continuous distributions. The resulting stochastic product demand specifications are used to obtain an optimal minimum inventory level, or safety stock, required by each distribution center to ensure an acceptable demand service level with the overall goal of minimizing total cost. Shen et al.49 detailed a stochastic inventory management system coupled with a supply chain network design problem. You and Grossmann50 showed that a similar MINLP formulation for combined supply chain network design and stochastic inventory management, namely

The minimum approach temperature parameter, Tmapp, acts as a restriction on the feasible region of the original MINLP. Consequently, this parameter has a significant impact on the rate of convergence to a global optimum. Table 3 compares Table 3. Computational Experience for Heat Integration Problems ID no.

Tmapp (K)

Mdyn (no./N)

tbil (s)

tnobil (s)

tnobil/tbil

1 2 3 4 5 6 7 8 9

1 3 5 7 9 10 11 13 15

39.16 35.47 40.29 33.49 36.99 44.51 49.35 45.79 35.56

>500.0 147.0 254.1 88.5 226.2 121.7 111.9 161.1 109.2

>500.0 377.9 >500.0 309.6 267.2 181.7 223.6 221.1 183.6

− 2.6 >2.0 3.5 1.2 1.5 2.0 1.4 1.7

solution times across nine parameter values over the range Tmapp ∈ [1,15], including the default value Tmapp = 10. Lower Tmapp values allow for a wide range of heat exchanger approach temperatures. When the feasible region of the original MINLP is expanded through the use of a small Tmapp parameter value, as in instance 1, the resulting increase in computational complexity causes an increase in required solution time. Each parameter value in Table 3 is assigned an ID number. Here, Mdyn represents the number of dynamic cutting planes applied per node. Solution times with and without the application of bilinear cuts are labeled tbil and tnobil, respectively. A relative comparison of solution times, constituting acceleration, is included in the final column. The application of dynamic cutting planes resulted in accelerated convergence for all combinations of parameter values. Further evidence of the beneficial effects of bilinear cutting planes during the solution of heat exchanger network synthesis problems is provided in Table 4.

J

Tmapp (K)

Nbil

Nnobil

Nnobil/ Nbil

2 4 5 6 7 8 9

3 7 9 10 11 13 15

4840 3242 7519 4405 3596 4905 5212

10952 9661 9197 5596 7764 7027 6268

2.3 3.0 1.2 1.3 2.2 1.4 1.2

S

bil

310 260 560 407 348 533 303

S

nobil

849 809 520 383 485 653 610

J

J

min ∑ f j xj + β ∑ ∑ dijχμi yij + θhzα ∑ j=1

Table 4. Cut Generation Statistics for Heat Integration Problems ID no.

I

i=1 j=1

j=1

I

∑ Lσi 2yij i=1

⎧ ⎡ ⎛ θh Dj ⎞⎫ ⎛ Dj ⎞⎤ ⎪ ⎪ ⎟⎟⎬ + ∑ xj⎨Fjnj + β ⎢gj + ⎜⎜aj ⎟⎟⎥nj + ⎜⎜ ⎪ ⎢⎣ ⎝ 2 nj ⎠⎪ ⎝ nj ⎠⎥⎦ j=1 ⎩ ⎭ J

Snobil/ Sbil 2.7 3.1 0.9 0.9 1.4 1.2 2.0

(40a)

J

s.t.

∑ yij = 1

i ∈ {1, ..., I } (40b)

j=1

yij ≤ xj

i ∈ {1, ..., I }, j ∈ {1, ..., J }

(40c)

I

Dj = χ ∑ μi yij

j ∈ {1, ..., J } (40d)

i=1

Instances 1 and 3 were omitted from this table because one or more solvers failed to converge within the 500-s time limit. The number of nodes required to converge and prove global optimality with and without the application of bilinear cutting planes are labeled Nbil and Nnobil, respectively. Sbil and Snobil represent the maximum numbers of nodes stored simultaneously during the branch-and-bound algorithm. The ratios Nnobil/Nbil and Snobil/Sbil indicate relative iteration and storage requirements. Although instances 5 and 6 exhibited a relatively small increase in maximum storage requirements after the addition of bilinear cutting planes, the use of cutting planes resulted in fewer required iterations in all cases.

xj , yij ∈ {0, 1}

Dj , nj ≥ 0

i ∈ {1, ..., I }, j ∈ {1, ..., J }

j ∈ {1, ..., J }

(40e) (40f)

suffers from a weak standard factorable relaxation. Here, decision variables xj and yij indicate the activity of each distribution center j and service activity between each distribution center j and each retailer i, respectively. Variable Dj represents the expected annual demand needs of each distribution center j, and nj represents the number of inventory replenishment orders fulfilled each year by a supplier. H



Article

Figure 4. Description of the supply chain with inventory management problem.

Parameter f j signifies the fixed cost of constructing a distribution center at location j, and Fj represents the fixed cost of placing a replenishment order to each distribution center. Parameters gj and aj indicate the fixed and unit transportation costs from a supplier to each distribution center, respectively. Similarly, dij is the unit transportation cost from each distribution center j to each supplier i. L is the lead time on an order placed with a supplier, and χ is the number of business days per year. Stochastic demand at each retailer i is represented by a probability distribution with a mean μi and variance σi2. Parameter h represents inventory holding costs, and zα is a service level satisfying Pr(z ≤ zα) = α.53 Finally, β and θ are weight factors assigned to transportation and inventory costs, respectively. The objective in eq 40a minimizes the total fixed and transportation costs of the supply chain network. The constraint in eq 40b ensures that each retailer i is serviced by a single distribution center j. Equation 40c prevents the transportation of products from inactive distribution centers. Finally, eq 40d translates daily retailer demands to calculate the total annual needs of each distribution center. 5.2. Cutting Planes. Standard factorable relaxation of the supply chain and stochastic inventory management MINLP formulation results in a weak linear relaxation. The addition of both static lifted cutting planes and dynamic disaggregation cutting planes tightens the standard relaxation. Here, we detail cuts added by the proposed methodology that exploit the underlying bilinear structure of the original mixed-integer nonlinear problem. The linearization variables t1j =

Dj nj

This new lifted static cutting plane 43 reveals a relationship between the intermediate bilinear and linear relationships (eqs 41 and 42, respectively) in the original nonlinear objective and the number of annual replenishments. The additional linearization variables t3j = t 2jnj

t4j = Fjnj + βt j3 +

t5j = xjt4j

j ∈ {1, ..., J }

j ∈ {1, ..., J }

(45) (46)

t3j ≥ gjnj + aj max{njt1j L + njLt1j − njLt1j L , njt1jU + njUt1j − njUt1jU} (47)

t3j ≤ gjnj + aj min{njt1jU + njLt1j − njLt1jU , njt1j L + njUt1j − njUt1j L} (48)

t5j ≥ Fj max{xjnjL + xjLnj − xjLnjL , xjnjU + xjUnj − xjUnjU} + β max{xjt3j L + xjLt3j − xjLt3j L , xjt3jU + xjUt3j − xjUt3jU} +

θh max{xjt1j L + xjLt1j − xjLt1j L , xjt1jU + xjUt1j − xjUt1jU} 2 (49)

t j5 ≤ Fj min{xjnjU + xjLnj − xjLnjU , xjnjL + xjUnj − xjUnjL} + β min{xjt3jU + xjLt3j − xjLt3jU , xjt3j L + xjUt3j − xjUt3j L}

(41)

+

j ∈ {1, ..., J }

(42)

j ∈ {1, ..., J }

θh min{xjt1jU + xjLt1j − xjLt1jU , xjt1j L + xjUt1j − xjUt1j L} 2 (50)

for j ∈ {1, ..., J} are available as disaggregation cutting planes to tighten relaxations dynamically throughout the branch-andbound tree. These new cutting planes in eqs 47−50 reveal relationships between the introduced variables tj1, tj3, and tj5 and the original problem variables xj and nj. Selected dynamic cutting planes are applied during bounding and range reduction at each node in the branch-and-bound tree.8 5.3. Computational Results. Here, we analyze the effectiveness of static and dynamic cutting planes applied by

are introduced into the objective function in eq 40a during sequential factorable reformulation. The resulting bilinear identity 41 and linear identity 42 are added as intermediate constraints in the standard relaxation. The linear relationship gjnj + ajw1j − w2j = 0

θh 1 tj 2

(44)

are introduced into the objective in eq 40a and added to the relaxation as factorable reformulation continues. The intermediate relationships in eqs 44 and 46 include a factored functional form similar to 15. The linear relationships

j ∈ {1, ..., J }

t 2j = gj + ajt1j

j ∈ {1, ..., J }

(43)

is generated by lifting the intermediate constraint in eq 42 with the number of replenishments nj, where w1j = tj1nj and w2j = tj2nj. I



Article

Table 5. Computational Experience for Supply Chain and Inventory Management ID no.

β

θ

Mstat

Mdyn (no./N)

tbil (s)

tnobil (s)

tP2 (s)

tnobil/tbil

tP2/tbil

1 2 3 4 5 6 7 8 9 10 11 12

0.01 0.01 0.01 0.01 0.01 0.01 0.10 0.10 0.10 0.10 0.10 0.10

0.001 0.005 0.010 0.050 0.100 0.500 0.001 0.005 0.010 0.050 0.100 0.500

3 3 3 3 3 3 3 3 3 3 3 3

15.33 39.51 30.76 21.66 19.99 17.82 16.00 12.00 32.00 35.00 28.25 35.66

0.3 1.5 3.2 8.9 13.2 4.0 0.1 0.1 0.1 0.1 0.5 3.3

12.9 29.6 98.7 55.9 112.2 9.5 71.1 24.0 9.5 3.8 7.7 24.4

5.0 5.1 4.3 0.8 0.7 0.9 5.1 5.0 5.1 0.7 3.8 0.9

49.7 20.4 31.0 6.3 8.5 2.4 889.0 300.5 105.2 41.8 16.3 7.4

19.2 3.5 1.4 0.1 0.1 0.2 63.8 62.5 56.7 7.8 8.1 0.3

Table 6. Cut Generation Statistics for Supply Chain and Inventory Management ID no.

β

θ

Nbil

Nnobil

Nnobil/Nbil

1 2 3 4 5 6 7 8 9 10 11 12

0.01 0.01 0.01 0.01 0.01 0.01 0.10 0.10 0.10 0.10 0.10 0.10

0.001 0.005 0.010 0.050 0.100 0.500 0.001 0.005 0.010 0.050 0.100 0.500

21 129 379 1298 1938 559 1 1 1 1 51 351

3515 7649 23135 13783 21836 1497 45337 15511 5143 1701 2373 5203

167.4 59.3 61.0 10.6 11.3 2.7 45337.0 15511.0 5143.0 1701.0 46.5 14.8

S

bil

4 20 40 148 170 50 1 1 1 1 6 34

S

nobil

206 611 1817 1207 1759 141 2071 590 208 56 147 374

Snobil/Sbil 51.5 30.6 45.4 8.2 10.3 2.8 2071.0 590.0 208.0 56.0 24.5 11.0

combination of weighting parameters. The final two columns represent the relative solution times, or accelerations, achieved. For all combinations of parameter values, Model_P2 is solved more rapidly than Model_P1 without the aid of the proposed bilinear cutting plane algorithms. For parameter combinations 4, 5, 6, and 12, formulation Model_P2 outperforms Model_P1 after the addition of bilinear cutting planes. For the remaining eight parameter combinations, however, Model_P1 dominates Model_P2 when the proposed cutting plane algorithms are enabled. Interestingly, the Model_P2 formulation exhibits enhanced performance when inventory costs outweigh transportation cots, whereas bilinear cutting planes are more effective when transportation costs are at least as important as inventory costs. Across all parameter combinations, bilinear cutting planes result in an expedited convergence of Model_P1 to a globally optimal solution. Additional cutting plane statistics are available in Table 6. In this table, Nbil and Nnobil identify the numbers of nodes required to converge to a globally optimal solution with and without the application of bilinear cutting planes, respectively. The relative number of nodes Nnobil/Nbil indicates whether the application of cutting planes results in significantly fewer BARON iterations. Sbil and Snobil signify the maximum numbers of nodes stored simultaneously throughout the branch-andbound algorithm. The final column displays the relative storage requirements Snobil/Sbil. The application of static and dynamic cutting planes results in a significant reduction in iterations, as well as decreased storage requirements. Further evidence of relaxation tightening by the proposed cutting planes is displayed in Figure 5. This figure compares root node relaxation gaps as calculated by eq 24. A modest

our implementation during the solution of supply chain design problems incorporating stochastic inventory management. The formulation tested, Model_P1, is available at MINLP.org54 and has been shown to produce a weak standard factorable relaxation.50 Typically, these challenging problems are solved by decomposition methods after transformation.49,50 In particular, a tailored reformulation Model_P2 available at MINLP.org54 was shown to outperform Model_P1 with a previous version of the BARON solver.50 For completeness, we compare the performance of the proposed cutting plane algorithms to this tailored reformulation. The weighting parameters β and θ in the objective in eq 40a were shown by You and Grossmann50 to have a significant effect on computational performance. Table 5 compares solution times across 12 combinations of parameter values, including the default parameter values β = θ = 0.01. All instances were solved in the same test environment as described in section 3.3, and with the exception of enabling and disabling the bilinear cutting plane algorithms and setting the absolute and relative termination tolerances to 10−6, all GAMS/BARON options were set to their default values. Each parameter combination in Table 5 is assigned an ID number. In this table, Mstat represents the number of static cutting planes applied during preprocessing, and M dyn represents the average number of dynamic cutting planes applied at each node of the branch-and-bound tree. tbil and tnobil are solution times with and without the application of bilinear cutting planes, respectively. Because no bilinear cutting planes are available during the soluiton of Model_P2, tbil = tnobil, and the notation tP2 is used to represent the solution time for each J



Article

relative termination tolerances to 10−6, all GAMS and solver options were set to their default values. The BARON solver demonstrates a clear performance advantage in instances from the application areas of trim loss, heat integration, and supply chain design. Figure 6 also shows that the proposed bilinear cutting plane algorithms have a significant impact on BARON’s performance in these application areas. Problem formulations for each of these applications exhibit underlying bilinear substructures. The discussions in sections 3.2, 4.2, and 5.2 reveal these exploitable substructures and detail families of cutting planes applied by BARON to enhance bounding and range reduction. Although shown to result in weak standard factorable programming relaxations, the numerical results in sections 3.3, 4.3, and 5.3 show that standard relaxations are enhanced by the addition of bilinear cutting planes. The computational results demonstrate the reduction of iterations, solution times, and storage requirements required to converge to global optimality. Accompanied by detailed descriptions of the families of cuts responsible for enhanced performance, these results demonstrate the impact of bilinear cutting planes on process synthesis and process operations applications.

Figure 5. Root node gap comparison for supply chain problems.

number of static lifted cuts is responsible for the gap reductions presented here. In particular, instances 7−10 exhibit a 100% gap reduction after the application of bilinear cutting planes and are solved at the root node before beginning the branch-andbound algorithm. For all instances, dynamic disaggregation cuts speed the solution process deeper in the branch-and-bound tree.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

6. SUMMARY In this work, we have examined the effectiveness of cutting planes that exploit underlying bilinear structures arising in instances from process synthesis and process operations. In particular, we highlight performance enhancements achieved as a result of the proposed cuts by the BARON solver in the areas of trim loss, heat integration, and supply chain optimization. For completeness, a performance comparison of other local and global optimization solvers across all instances form this work is included in Figure 6.

Notes

The authors declare the following competing financial interest(s): The second author has ownership in the company that provides BARON software.

■

REFERENCES

(1) Horst, R., Pardalos, P. M., Eds. Handbook of Global Optimization; Nonconvex Optimization and Its Applications Series; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1995; Vol. 2. (2) Kearfott, R. B. Rigorous Global Search: Continuous Problems; Nonconvex Optimization and Its Applications Series; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; Vol. 13. (3) Van Hentenryck, P.; Michel, L.; Deville, Y. Numerica: A Modeling Language for Global Optimization; The MIT Press: Cambridge, MA, 1997. (4) Floudas, C. A. Deterministic Global Optimization: Theory, Algorithms and Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1999. (5) Tawarmalani, M.; Sahinidis, N. V. Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming: Theory, Algorithms, Software, and Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2002. (6) Tawarmalani, M.; Sahinidis, N. V. A polyhedral branch-and-cut approach to global optimization. Math. Program. 2005, 103, 225−249. (7) Tawarmalani, M.; Sahinidis, N. V. Global optimization of mixedinteger nonlinear programs: A theoretical and computational study. Math. Program. 2004, 99, 563−591. (8) Zorn, K.; Sahinidis, N. V. Global optimization of general nonconvex problems with intermediate bilinear structures. Optim. Methods Software, published online Apr 23, 2013, 10.1080/ 10556788.2013.783032. (9) GLOBAL Library. Available at http://www.gamsworld.org/ global/globallib.htm (accessed May 2013). (10) Bussieck, M. R.; Drud, A. S.; Meeraus, A. MINLPLibA collection of test models for mixed-integer nonlinear programming. INFORMS J. Comput. 2003, 15, 114−119. (11) Sherali, H. D.; Adams, W. P. A hierarchy of relaxations between the continuous and convex hull representations for zero-one programming problems. SIAM J. Discrete Math. 1990, 3, 411−430.

Figure 6. Comparison of local and global solver performances.

Figure 6 contains performance profiles for local and global mixed-integer nonlinear solvers including AlphaECP,55 BARON,6 DICOPT,56 Couenne,57 LINDOGlobal,58 SCIP,59 and SBB.60 All of the results in this figure were generated using GAMS 23.8.229 on a 64-bit Intel Xenon X5650 2.66 Ghz processor. With the exception of reducing the absolute and K



Article

(12) Adams, W. P.; Sherali, H. D. Linearization strategies for a class of zero-one mixed integer programming problems. Oper. Res. 1990, 38, 217−226. (13) Sherali, H. D.; Tuncbilek, C. H. A reformulation− convexification approach for solving nonconvex quadratic programming problems. J. Global Optim. 1995, 7, 1−31. (14) Sherali, H. D.; Alameddine, A. A new reformulation− linearization technique for bilinear programming problems. J. Global Optim. 1992, 2, 379−410. (15) Quesada, I.; Grossmann, I. E. Global optimization of bilinear process networks and multicomponent flows. Comput. Chem. Eng. 1995, 19, 1219−1242. (16) Tawarmalani, M.; Ahmed, S.; Sahinidis, N. V. Product disaggregation and relaxations of mixed-integer rational programs. Optim. Eng. 2002, 3, 281−303. (17) McCormick, G. P. Computability of global solutions to factorable nonconvex programs: Part IConvex underestimating problems. Math. Program. 1976, 10, 147−175. (18) Al-Khayyal, F. A.; Falk, J. E. Jointly constrained biconvex programming. Math. Oper. Res. 1983, 28, 273−286. (19) Kantorovich, L. V. Mathematical methods of organizing and planning production. Manage. Sci. 1960, 6, 366−422. (20) Gilmore, P. C.; Gomory, R. E. A linear programming approach to the cutting stock problem. Oper. Res. 1961, 9, 849−859. (21) Gilmore, P. C.; Gomory, R. E. A linear programming approach to the cutting stock problemPart II. Oper. Res. 1963, 11, 863−888. (22) Gilmore, P. C.; Gomory, R. E. Multistage cutting stock problems of two and more dimensions. Oper. Res. 1965, 13, 94−120. (23) Golden, B. L. Approaches to the cutting stock problem. AIIE Trans. 1976, 8, 265−274. (24) Hinxman, A. The trim loss and assortment problems: A survey. Eur. J. Oper. Res. 1980, 5, 8−18. (25) Dyckhoff, H. A typology of cutting and packing problems. Eur. J. Oper. Res. 1990, 44, 145−159. (26) Lodi, A.; Martello, S.; Monaci, M. Two-dimensional packing problems: A survey. Eur. J. Oper. Res. 2002, 141, 241−252. (27) Wäscher, G.; Haußner, H.; Schumann, H. An improved typology of cutting and packing problems. Eur. J. Oper. Res. 2007, 183, 1109−1130. (28) Harjunkoski, I.; Westerlund, T.; Pörn, R.; Skrifvars, H. Different transformations for solving non-convex trim loss problems by MINLP. Eur. J. Oper. Res. 1998, 105, 594−603. (29) Brooke, A.; Kendrick, D.; Meeraus, A. GAMSA User’s Guide; Scientific Press: Redwood City, CA, 1988. (30) Sahinidis, N. V. BARON: A general purpose global optimization software package. J. Global Optim. 1996, 8, 201−205. (31) Shenoy, U. V. Heat Exchanger Network Synthesis: Process Optimization by Energy and Resource Analysis; Gulf Publishing Company: Houston, TX, 1995. (32) Smith, R. Chemical Process Design; McGraw−Hill: New York, 1995. (33) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall: Upper Saddle River, NJ, 1997. (34) Seider, W. D.; Seader, J. D.; Lewin, D. R. Process Design Principles: Synthesis Analysis and Evaluation; John Wiley & Sons: New York, 1999. (35) Gundersen, T.; Naess, L. The synthesis of cost optimal heat exchanger networks. Comput. Chem. Eng. 1988, 12, 503−530. (36) Jezowski, J. Heat exchanger network grassroot and retrofit design. The review of the state-of-the-art: Part I. Heat exchanger network targeting and insight based methods of synthesis. Hung. J. Ind. Chem. 1994, 22, 279−294. (37) Jezowski, J. Heat exchanger network grassroot and retrofit design. The review of the state-of-the-art: Part II. Heat exchanger network synthesis by mathematical methods and approaches for retrofit design. Hung. J. Ind. Chem. 1994, 22, 295−308.

(38) Furman, K. C.; Sahinidis, N. V. A critical review and annotated bibliography for heat exchanger network synthesis in the 20th century. Ind. Eng. Chem. Res. 2002, 41. (39) Yee, T. F.; Grossmann, I. E.; Kravanja, Z. Simultaneous optimization models for heat integrationI. Area and energy targeting and modeling of multi-stream exchangers. Comput. Chem. Eng. 1990, 14, 1151−1164. (40) Yee, T. F.; Grossmann, I. E. Simultaneous optimization models for heat integrationII. Heat exchanger network synthesis. Comput. Chem. Eng. 1990, 14, 1165−1184. (41) Daichendt, M. M.; Grossmann, I. E. Preliminary screening procedure for the MINLP synthesis of process systemsII. Heat exchanger network synthesis. Comput. Chem. Eng. 1994, 18, 679−709. (42) Daskin, M. S. Network and Discrete Location: Models, Algorithms, and Applications; Wiley: New York, 1995. (43) Owen, H.; Daskin, M. S. Strategic facility location: A review. Eur. J. Oper. Res. 1998, 111, 423−447. (44) Cachon, G. P.; Fisher, M. Supply chain inventory management and the value of shared information. Manage. Sci. 2000, 6, 1032−1048. (45) Tsiakis, P.; Shah, N.; Pantelides, C. C. Design of multi-echelon supply chain networks under demand uncertainty. Ind. Eng. Chem. Res. 2000, 38, 1279−1290. (46) Kok, A. G.; Graves, S. C. Supply Chain Management: Design, Coordination and Operation; Elsevier: New York, 2003. (47) Grossmann, I. E. Enterprise-wide optimization: A new frontier in process systems engineering. AIChE J. 2005, 51, 1846−1857. (48) Chopra, S.; Meindl, P. Supply Chain Management: Strategy, Planning and Operation; Prentice Hall: Upper Saddle River, NJ, 2003. (49) Shen, Z.-J. M.; Coullard, C.; Daskin, M. S. A joint location− inventory model. Transp. Sci. 2003, 37, 40−55. (50) You, F.; Grossmann, I. E. Mixed-integer nonlinear programming models and algorithms for large-scale supply chain design with stochastic inventory management. Ind. Eng. Chem. Res. 2008, 47, 7802−7817. (51) You, F.; Grossmann, I. E. Integrated multi-echelon supply chain design with inventories under uncertainty: MINLP models, computational strategies. AIChE J. 2010, 56, 419−440. (52) You, F.; Grossmann, I. E. Balancing responsiveness and economics in the design of process supply chain with multi-echelon stochastic inventory. AIChE J. 2011, 57, 178−192. (53) Zipkin, P. H. Foundations of Inventory Management; McGrawHill: Boston, MA, 2000. (54) CMU-IBM Cyber-Infrastructure for MINLP collaborative site. Available at http://www.minlp.org (accessed May 2013). (55) Westerlund, T.; Lastusilta, T. AlphaECP Solver Manual, 2012; available at http://www.gams.com/dd/docs/solvers/alphaecp.pdf (accessed May 2013). (56) Grossmann, I. E.; Viswanathan, J.; Vecchietti, A. DICOPT Solver Manual, 2012; available at http://www.gams.com/dd/docs/solvers/ dicopt.pdf (accessed May 2013). (57) Belotti, P. COUENNE: A User’s Manual, 2009; available at https://projects.coin-or.org/Couenne/browser/trunk/Couenne/doc/ couenne-user-manual.pdf?format=raw (accessed May 2013). (58) LINDOGlobal Solver Manual, 2012; available at http://www. gams.com/dd/docs/solvers/lindo.pdf (accessed May 2013). (59) Vigerske, S. SCIP Solver Manual, 2012; available at http://www. gams.com/dd/docs/solvers/scip.pdf (accessed May 2013). (60) SBB Solver Manual, 2012. Available at http://www.gams.com/ dd/docs/solvers/sbb.pdf (accessed May 2013).

L


Computational Experience with Applications of Bilinear Cutting Planes

Recommend Documents