5424
Ind. Eng. Chem. Res. 2010, 49, 5424–5438
Global Optimization of Large-Scale Generalized Pooling Problems: Quadratically Constrained MINLP Models Ruth Misener and Christodoulos A. Floudas* Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544-5263
The generalized pooling problem is a generic way of identifying a topology containing sources nodes, intermediate storage tanks, and sinks when monitoring specific stream qualities is imperative. One important application domain of the generalized pooling problem is wastewater treatment. Choosing among the wide array of available wastewater treatment technologies is a combinatorially complex optimization problem that requires nonconvex terms to monitor regulated stream qualities about a treatment plant. In this work, we address five instantiations of the generalized pooling problem to global optimality by introducing (i) a quadratically constrained MINLP model formulation that reduces the number of bilinear terms, (ii) novel piecewise underestimation methods for the nonconvex bilinear terms to tighten the relaxation [Gounaris et al. Ind. Eng. Chem. Res. 2009, 48, 5742-5766; Wicaksono and Karimi AIChE J. 2008, 54, 991-1008], and (iii) a branch-and-bound algorithm suited to address the combinatorial complexity of the problem. Extensive computational results are presented for small, medium, large, and two very large-scale case studies which feature 48, 300, 675, and 1260 distinct bilinear terms, respectively. We show that the small, medium, and large case studies can be solved to global optimality efficiently. The two very large case studies can be solved within 0.9% and 2.3% of the global optimum, and when the additional assumption of a limited number of plants is introduced, they can be solved to global optimality. 1. Introduction 3
The pooling problem, which we recently reviewed, arises from limited storage conditions at a chemical or petrochemical plant.4 Nonconvex bilinear terms in the problem arise from tracking the concentration of key components or qualities when multiple streams mix in intermediate storage nodes under the assumption of linear blending. Audet et al.5 and Meyer and Floudas6 introduced the generalized pooling problem, which increases the complexity of the pooling problem by transforming the network topology into a decision variable. Choosing the interconnections between process units and storage tanks, or pools, is combinatorially complex. Because the activation or deactivation of each pipe or intermediate node is a discrete decision and the linear mixing at the intermediate nodes leads to bilinear terms, the generalized pooling problem is a mixedinteger nonconvex program (nonconvex MINLP) with quadratic equalities and inequalities which exhibits multiple locally optimal solutions. The major challenge inherent in this problem is the development of rigorous global optimization methods that can address large-scale problems to global optimality. The generalized pooling problem is a fundamentally generic problem that can effectively model multiple applications in refineries, chemical plants, and water treatment facilities. These challenges are linked by the need to monitor specific stream qualities to meet environmental standards (e.g., by limiting heavy metals in wastewater or reducing sulfur in refinery endproducts) or to ensure product quality (e.g., by tracking octane in a gasoline blend). The need to address large, industrially relevant instantiations of the pooling and water network problems have led to a sustained interest in bilinear programs, which are more generally called quadratically constrained quadratic programs (QCQP).7-9 For QCQPs that arise in various * To whom all correspondence should be addressed. E-mail:
[email protected]. Tel.: (609) 258-4595. Fax: (609) 2580211).
process synthesis applications, see the work of Floudas and co-workers.10-14 The application we address in this paper is the design of wastewater treatment networks. Wastewater generated in a refinery or chemical plant is treated to remove pollutants before being released into the environment. Choosing among the wide array of available treatment technologies (e.g., physical, chemical, biological, and thermal removal techniques) and selecting the most profitable interconnections between these processing units is a combinatorially complex optimization problem.15 The industrial need to reduce water treatment costs while meeting standards such as the ones set by the United States Clean Water Act of 1977 has generated intense interest in the optimal design of water networks.15 The two mathematical programming challenges of designing water networks are (1) dealing with the combinatorial complexity of choosing a single topology from a superstructure and (2) addressing the nonconvex bilinear terms resulting from mixing several process streams in a single treatment unit. In this paper, we revisit the generalized pooling problem test cases originally posed by Meyer and Floudas.6 These test cases posit a set of sources containing regulated qualities that must be treated before release into the environment. Representing the challenges that industry faces, we consider as many as twenty treatment options and allow the possibility of interconnections between all the treatment plants. We exploit recent advances in piecewise-linear underestimation of bilinear terms within a branch-and-bound algorithm and globally optimize these test cases.1,2 We assume a command-and-control policy on environmental pollutants and place profitability in the objective function. Malcolm et al.16 recently compared the economic and environmental differences between regulatory policies (i.e., commandand-control, environmental taxes, and cap-and-trade), and Lim and Park17 have addressed the idea of minimizing environmental impact rather than cost. The mathematical complexity of the
10.1021/ie100025e 2010 American Chemical Society Published on Web 05/04/2010
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
wastewater treatment problem presented in this paper does not change between regulatory frameworks and could be easily transformed to represent different policies or objective functions. We follow Meyer and Floudas6 in assuming a linear objective function with a fixed cost associated with activating each pipeline or treatment plant and a variable cost representing the flow through each pipeline or treatment plant, but this formulation could be easily modified to admit a nonlinear investment cost objective. For cases with a concave investment cost, we have previously shown that piecewise-linear underestimators can effectively address the additional complexity.18,19 After reviewing the literature in section 2, section 3 will introduce the problem formulation and the specific test cases. Section 4 compares possible global optimization schemes, section 5 presents the optimization results, and section 6 concludes the paper. Test case parameters can be found in Appendix B. 2. Literature Review Takama et al.20 addressed the need to integrate water-using and water-treating processes within one large water network. But, despite the recognition that globally optimizing a mathematically complete superstructure leads to the greatest cost savings, the difficulty of tackling industrially sized MINLPs led researchers to develop clever approaches yielding near-optimal results. Mathematical programming-related work on the water networks problem includes that of Doyle and Smith,21 who developed a method to find a feasible point and initialize a local sequential quadratic programming (SQP) search; Alva-Arga´ez et al.,22 who recursively solved an MILP that forced problem infeasibilities to zero; Galan and Grossmann,23 who designed a series of linear programs to seed a multistart local search algorithm; Huang et al.,24 who used embedded topological constraints and current flowsheet data to guess a topology before initializing a local search; Gunaratnam et al.,25 who moved the infeasibility penalty terms of Alva-Arga´ez et al.22 into the constraint set; Bogataj and Bagajewicz,26 who worked on effectively initializing a local solver for a heat integrated water network; Chakraborty,27 who studied a special case with only sources and sinks that can be formulated as an MILP; Faria and Bagajewicz,28 who initialized a local solution of a complete water system using a linear relaxation; and Ponce-Ortega et al.,29 who reduced the number of bilinear terms by avoiding process stream mixing in intermediate nodes. The water networks problem has been alternatively addressed using graphical and heuristic approaches (e.g., Kuo and Smith30), but that work is it out of the scope of this discussion. The first twenty years of research concerning water networks is comprehensively discussed in the review article of Bagajewicz.15 In contrast to the water networks literature, studies addressing the pooling problem have focused almost exclusively on global optimization since the early 1990s. Haverly31 used a local search algorithm, Lasdon et al.32 employed successive linear programming to address pooling problem instantiations, and Floudas et al.33 and Floudas and Aggarwal34 designed an algorithm based on generalized Benders’ decomposition. Floudas and Visweswaran developed the first rigorous deterministic global approach for biconvex and bilinear problems based on duality theory.35-38 Foulds et al.39 implemented the bilinear envelopes of McCormick.40 Ben-Tal et al.41 introduced the q-formulation and applied it to pooling problem case studies. Adhya et al.42 explored Lagrangian approaches. Quesada and Grossmann43 used the reformulation-linearization technique (RLT) of Sherali and Alameddine.44 Also see Sherali and Adams45 for a comprehen-
5425
sive study of the reformulation-linearization technique. Tawarmalani and Sahinidis46 proposed the pq-formulation. Audet et al.5 applied the branch-and-cut RLT-based QCQP algorithm of Audet et al.47 Almutairi and Elhedhli48 suggested a new Lagrangian relaxation for the pooling problem and demonstrated that their relaxation is often tighter than previously developed Lagrangian relaxations. Misener et al.49 extended the pooling problem to incorporate the Environmental Protection Agency complex emissions model and associated legislation into the constraint set. Global optimization algorithms for bilinear programs include the branch-and-bound method of Al-Khayyal and Falk50 that relaxes the bilinear terms using convex envelopes, the duality-based global optimization algorithm of Floudas and Visweswaran,35,38 the RLT-based branch-and-cut algorithm of Audet et al.,47 the branch-and-bound procedure of Linderoth51 that generates tight relaxations by partitioning two-dimensional regions into triangles and rectangles, the integration of RLT and semidefinite programming relaxations by Anstreicher,52 and the cutting plane algorithm of Bao et al.53 that generates a multiterm relaxation of a QCQP. Recently, the water network and pooling literature have converged in the recognition that industrially sized bilinear programs can be globally optimized using piecewise-linear relaxations. In this paper, we focus on the ab initio piecewise application of convex envelopes which tighten the linear relaxation, thus reducing the number of nodes needed in a branch-and-bound tree. This idea was first introduced by Meyer and Floudas6 and Karuppiah and Grossmann54 in the context of tightly underestimating large-scale wastewater treatment and water networks problems, respectively. Recognizing the importance of formulating the piecewiselinear relaxations of Meyer and Floudas6 and Karuppiah and Grossmann54 in the best manner possible, Wicaksono and Karimi2 introduced 15 mathematically equivalent alternative formulations and compared their performance on several test cases. We recently proposed five additional piecewise-linear formulations and conducted a comprehensive comparative study on the computational performance of the piecewise linear formulations on a collection of benchmark pooling problems.1 Other groups who have used piecewise-linear underestimators are Bergamini et al.,55 who expedited the convergence of their outer approximation for global optimization algorithm using the piecewise envelopes; Saif et al.,56 who used piecewise underestimation of bilinear terms to globally optimize a reverse osmosis network; Pham et al.,57 who coupled the piecewise underestimators with a fast-solving algorithm to generate nearoptimal solutions; Hasan and Karimi,58 who partitioned both variables appearing in each bilinear term; and Misener et al.,49 who used piecewise envelopes to address the extended pooling problem. 3. Problem Formulation The generalized pooling problem formulation presented here, which is adapted from Meyer and Floudas,6 represents the water treatment topology illustrated in the Figure 1. The only topological restriction we place on this formulation is that the treated wastewater cannot loop between two nodes or within a single node. This generic superstructure allows the globally optimal solution to determine the topology that meets environmental restrictions in the most cost-effective way possible. The source nodes are process effluents, the tanks are treatment plants with various removal technologies, and the regulated sinks are the environment. The source and sink nodes are fixed, but planners can decide whether or not to activate each treatment
5426
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
yct,t' + yct',t e 1
∀t ∈ T; t' ∈ T\{t}
(7)
The flow rates are balanced at the source (eq 8), treatment (eq 9), and sink nodes (eq 10): f ssource )
∑a
s,e
∑d
+
∑d
s,t
∑
+
s∈S
∀s ∈ S
s,t
e∈E
(8)
t∈T
∑
ct',t )
t'∈T\{t}
∑b
ct,t' +
∀t ∈ T
t,e
t'∈T\{t}
e∈E
(9)
∑f
source s
∑a
)
s,e
s∈S
∑b
+
(10)
t,e
s∈S e∈E
t∈T e∈E
The flow rate into each of the sinks is controlled by an upper bound on the capacity:
∑a
Figure 1. Superstructure of the mid-size generalized pooling test case.
plant and pipeline. The objective is to minimize the cost of operation, that is, the sum of the fixed and variable costs associated with using a treatment plant and pipeline. In addition to network topology constraints, the concentrations of the regulated qualities are bounded by environmental standards at the sink nodes. Table 1 displays the problem nomenclature, sections 3.1-3.4 introduce the model, and section 3.5 describes the five test cases. 3.1. Objective. The objective is to minimize the costs of ya yc yd ye operation. There is a fixed cost (cs,e , cyb t,e , ct,t′, cs,t , ct ) for a activating a pipeline or treatment plant and a variable cost (cs,e , c d cbt,e, ct,t′ , cs,t , cte) for the flow rate through each pipeline and treatment plant: min zp)
∑ (c
a s,e · as,e
∑ (c
ya a + cs,e · ys,e )+
s∈S e∈E
∑
c yc c (ct,t' · ct,t' + ct,t' · yt,t' )+
(cte · [
t∈T
∑
ct,t' +
t'∈T\{t}
∑
∑ (c
bt,e] + ctye · yte)
e∈E
3.2. Intranetwork Flows. Each pipe in the network (connecting source to sink, source to treatment plant, treatment plant to another treatment plant, and treatment plant to sink) is activated or deactivated by a binary decision variable and bounded by topologically feasible flow rates: a a ys,e · a_s,e e as,e e ys,e · ajs,e
∀s ∈ S; e ∈ E
(2)
ybt,e · b_t,e e bt,e e ybt,e · bjt,e
∀t ∈ T; e ∈ E
(3)
∀t ∈ T; t' ∈ T\{t}
(4)
∀s ∈ S; t ∈ T
(5)
yct,t' · _ct,t' e ct,t' e yct,t' · jct,t' d d j ys,t · d_s,t e ds,t e ys,t · ds,t
The activation of each treatment plant is also controlled by a binary decision variable and constrained by maximum and minimum flow rates:
∑
t'∈T\{t}
ct,t' +
∑b
t,e
e yet · jet
∀t ∈ T
e∈E
Looping between two treatment plants is disallowed:
∑a
s,e
+
s∈S
(6)
∀e ∈ E
e f sink,max e
t,e
(11)
∑b
t,e)
∑q
source · as,e c,s
g
t∈T
s∈S
+
∑q
c,t · bt,e
t∈T
∀c ∈ C;e ∈ E (12)
Each monitored quality c is also tracked across the treatment nodes:
qc,t · [
yd d + cs,t · ys,t) +
t∈T s∈S
∑b t∈T
(1 - rc,t) · [
(1)
yet · _et e
qmax c,e · (
yb b + ct,e · yt,e) +
d s,t · ds,t
+
3.3. Quality Monitoring. The nonconvex bilinear terms in the generalized pooling problem arise from monitoring the regulated qualities at the sink nodes and across the treatment plants. At the sink nodes, the qualities are limited by applicable environmental laws (represented by qc,emax):
t∈T e∈E
t∈T t'∈T\{t}
∑
b t,e · bt,e
s,e
s∈S
∑
∑q
source · ds,t c,s
s∈S
ct,t' +
t'∈T\{t}
+
∑
qc,t' · ct',t] )
t'∈T\{t}
∑b
∀c ∈ C; t ∈ T
t,e]
(13)
e∈E
Equation 13 is logically equivalent to the quality balance in Meyer and Floudas,6 but it substantially reduces the number of nonconvex bilinear terms in the problem formulation. By eq 9, we can multiply qc,t on the right-hand side of eq 13 by either
∑d
+
∑
ct,t' +
s,t
s∈S
∑
ct',t
t'∈T\{t}
or
t'∈T\{t}
∑b
t,e
e∈E
But choosing the former option necessitates |C| · |T| · (2 · (|T| 1) + |S|) bilinear terms in the problem formulation (because qc,t · ct′,t and qc,t′ · ct′,t are distinct terms), while choosing the latter requires only |C| · |T| · (|T| - 1 + |E|) bilinear terms. If there are more effluent sources than environmental sinks, formulating the model as presented here leads to more than a 50% reduction in the number of distinct bilinear terms. In the case of the Meyer and Floudas6 case studies, the number of bilinear terms reduces from 156 to 48 for the four-plant test case and from 750 to 300 for the ten-plant test case. The primary advantage in reducing the number of bilinear terms is the average computation time saved in solving a termwise relaxation of the QCQP at each branch-and-bound tree node. An additional advantage we observed for the Meyer and
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5427
Table 1. Wastewater Treatment Problem Nomenclature s∈S t∈T e∈E c∈C n∈N f ssource f min f esink,max source qc,s max qc,e rc, t [e_t, jet] ya cs,e a cs,e yb ct,e b ct,e yc ct,t′ c ct,t′ yd cs,t d cs,t ye ct cte a ys,e b yt,e c yt,t′ d ys,t yte as,e ∈ 0 ∪ [a_s,e, ajs,e] bt,e ∈ 0 ∪ [b_t,e, bjt,e] ct,t′ ∈ 0 ∪ [c_t,t′, jct,t′] ds,t ∈ 0 ∪ [d_s,t, djs,t] qc,t ∈ [qc,t, qjc,t]
indices
parameters
binary variables
contin variables
Floudas6 four-plant test case was a tighter problem relaxation. We motivate this improved relaxation in Appendix A. 3.4. Hard Bounds. Topological network analysis leads to the hard bounds listed in eqs 14-19. The minimum flow rate f min represents the smallest physically feasible flow rate. The flow rates emerging from each source are bounded by the effluent in the source and the flow rates between the treatment plants are limited by the total material flowing through the system. The concentration of the monitored qualities at the exit of each treatment plant can be no greater than if the most polluted effluent stream is treated at the plant. a_s,e ) f
min
(
ajs,e ) min min
b_t,e ) f min
bjt,e )
qmax c,e ·
source qc,s
c
∑f
∑f
source s
)
; f ssource
∀s ∈ S; e ∈ E (14) ∀t ∈ T; e ∈ E (15)
source s
s∈S
_ct,t' ) f min
jct,t' )
∑f
source s
∀t ∈ T; t' ∈ T\{t}
(16)
∀s ∈ S; t ∈ T
(17)
s∈S
d_s,t ) f min _et ) f min
djs,t ) f ssource jet )
∑f
source s
∀t ∈ T
(18)
s∈S
q_c,t ) 0
source qjc,t ) (1 - rc,t) · [max qc,s ] s∈S
∀c ∈ C; t ∈ T
effluent source nodes treatment plants regulated sinks qualities subject to environmental constraints number of partitions for piecewise relaxation effluent source s flow rate the minimum flow rate through any pipeline the maximum capacity of sink e concentration of quality c in effluent s concentration limit of quality c in sink e removal ratio of quality c in plant t flow limits through activated plant t fixed pipeline cost from source s to sink e variable flow cost from source s to sink e fixed pipeline cost from plant t to sink e variable flow cost from plant t to sink e fixed pipeline cost from plant t to t′ (t * t′) variable flow cost from plant t to t′ (t * t′) fixed pipeline cost from source s to plant t variable flow cost from source s to plant t fixed cost for treatment plant t variable flow cost through treatment plant t activates pipeline from source s to sink e activates pipeline from plant t to sink e activates pipeline from plant t to plant t′ activates pipeline from source s to plant t activates plant t flow from source s to sink e flow from plant t to sink e flow from plant t to plant t′ (t * t′) flow from source s to plant t concentration of quality c in plant t effluent
distance between bilinear term x · y and its convex relaxation scales with the distance between the variable bounds: 1 · (xj - _x) · (yj - _y) 4 This result is in keeping with the commonly made observation that tight bounds lead to tight convex relaxations. Because the qualities qc,t participate in both the (qc,t · bt,e) and (qc,t · ct,t′) bilinear terms, the bounds on qc,t are especially important. When the four-plant test case of Meyer and Floudas6 is relaxed using termwise application of the convex envelopes first proposed by McCormick40 and optimized using an MILP solver, the resulting lower bound on the global optimum is 7.219 × 105. But if an unnecessarily loose bound of (qjc,t ) maxs∈S qc,ssource) is imposed on the qc,t variables, the resulting lower bound is weakened by nearly 20% (to 5.825 × 105). 3.5. Problem Size and Test Cases. The full model of the generalized pooling problem can be characterized by eqs 1-19. Representing the number of entities in each set using norm nomenclature (i.e., |C| is the number of monitored qualities, |S| is the number of effluent source nodes, |T| is the number of treatment plants, and |E| is the number of regulated sinks): |C| no. equations f 2 2 + |S| + |T| · (|T| + |E|) +
(
)
|T| · |T| + |S| + 2 (|S| + |T|) · (|T| + |E|) + |T| · no. continuous variables f (|C| - 1) no. binary variables f (|S| + |T|) · (|T| + |E|) no. bilinear terms f |C| · |T| · (|T| - 1 + |E|)
(19) Establishing tight bounds on the variables participating in bilinear terms (i.e., bt,e, ct,t′, and qc,t) is important based on the result of Androulakis et al.,59 showing that the separation
Table 2 lists the specific sizes of the five test cases discussed in this paper. The smallest test case, which is identical to the simplified test case of Meyer and Floudas,6 has five available treatment plants (in Appendix B.1, this is represented as T )
5428
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 Table 3. Results of the Meyer and Floudas6 4-Plant Test Case Using CPLEX 1260
Table 2. Sizes of the Five Test Cases topology
no. variables
R)1
|C| |S| |T| |E| no. eqs contin binary no. bilinear terms 4-plant 10-plant 15-plant 20-plant FDC 20-plant DDC
3 3 3 3 3
7 7 7 7 7
4 10 15 20 20
1 1 1 2 2
150 516 986 1663 1663
63 207 382 634 634
55 187 352 594 594
48 300 675 1260 1260
{t3, t7, t9, t10}). The 10-plant test case is equivalent to the larger case study of Meyer and Floudas,6 (T ) {t1, t2, ..., t10}). The 15-plant test case includes five additional plants (T ) {t1, t2, ..., t15}). Finally, two large 20-plant test cases with fixed disposal costs (FDC) and distance-dependent disposal costs (DDC) are presented in Appendix B.2. Appendix B records the values of all relevant parameters. Table 2 shows that this model uses more continuous variables than Meyer and Floudas,6 used in their formulation of the same test cases. Meyer and Floudas6 chose to make the following substitutions which are valid in the case of exactly one sink (|E| ) 1): source as,e r f s
bt,e r
∑
t'∈T\{t}
∑d
s,t
t∈T
ct',t -
∑
t'∈T\{t}
ct,t' +
∑d
s,t
s∈S
In this work, we are not making those substitutions. In the former case, we do not substitute for as,e because we have found no computational advantage for doing so and we would prefer to keep the model generic for instances like the 20-plant test cases (when |E| > 1). In the latter case, we do not substitute for bt,e because the replacement would force us to include the bilinear terms (qc,t · ct′,t), (qc,t · ct,t′), and (qc,t · ds,t) in the formulation rather than just (qc,t · ct,t′) and (qc,t · bt,e) terms. Including the extra (|T| · |E|) continuous variables leads to the large reduction in the number of bilinear terms we achieved in section 3.3. In sum, the 4- and 10-plant test cases use 11 and 17 more continuous variables, respectively, than the Meyer and Floudas,6 formulation but require 108 and 450 fewer bilinear terms. 3.6. Comparison to Previous Work. Meyer and Floudas,6 addressed the generalized pooling problem using what was stateof-the-art hardware and software more than five years ago (the runs were completed in 2003). We want to demonstrate the advantage of using the proposed alternative formulation (see section 3), piecewise underestimators (section 4.3), and a good branch-and-bound optimization algorithm (sections 4.4 and 4.5), but we use faster computers and more advanced software than was available in 2003 and therefore cannot directly compare our results to those of Meyer and Floudas.6 Table 3 compares the timing results that Meyer and Floudas,6 achieved using the mixed integer linear programming solver CPLEX 7.0 within the GAMS modeling language on an HP J2240 to the results we obtain running the same code on GAMS/ CPLEX 12.1 using an Intel Core 2 Quad processor with four 2.83 GHz cores.60,61 We will judge our algorithms by comparing our timing studies in Table 4 to the final column of Table 3. The data in Table 3 represents root node relaxations of the 4-plant generalized pooling problem for different partitioning levels and R values. The piecewise underestimation formulation Meyer and Floudas6 use is equivalent to the bm formulation1,2 while we use the nf4r formulation1 in section 4.3. The two formulations should produce equivalent relaxations, but the root node relaxations in Table 3 differ from the root nodes in Table 4 because
N
relaxation
4 5 6 7 8 9
1.002 × 10 1.002 × 106 1.014 × 106 1.034 × 106 1.037 × 106 1.053 × 106 6
% gap
2003 CPU s
2009 CPU s
7.79 7.79 6.58 4.84 4.57 3.02
1788 6294 10302 13248 36227 111119
101.3 401.9 382.3 933.8 1578.7 1880.7
% gap
2003 CPU s
2009 CPU s
7.79 4.15 2.00 1.22
1742 9599 22261 285449
99.6 344.3 502.2 1561.6
R ) 1/2 N
relaxation
4 5 6 7
1.002 × 10 1.041 × 106 1.064 × 106 1.073 × 106 6
(1) As described in section 3.3, we altered the Meyer and Floudas6 quality balance formulation to the equivalent eq 13. Reducing the number of bilinear terms does not affect the global optimum, but it tightens the relaxation. (2) The Meyer and Floudas6 code contains an additional integer cut that eliminates looping between three treatment plants: yct,t' + yct',t'' + yct'',t e 2
∀t ∈ T; t' ∈ T\{t}; t'' ∈ T\{t, t'} (20)
We follow the model described in the Meyer and Floudas6 manuscript rather than the model in the code and therefore omit eq 20. Although this second-order looping constraint tightens the relaxation, it does not impact the global optimum of any of the five test cases. (3) Meyer and Floudas6 propose piecewise underestimation not only of the bilinear terms, but also of RLT relaxations similar to the ones presented in section 4.2. Although augmenting each RLT constraint tightens the relaxation, the additional piecewise underestimators also increase the solution time of the MILP. Including these tightening constraints would probably lead to fewer nodes in a branch-and-bound tree, but we tradeoff tightness of each node for the time needed to solve a node and only piecewise underestimate each of the bilinear terms. (4) As we describe in Gounaris et al.,1 irregular partitioning may lead to improvements in the relaxation. Thus, Meyer and Floudas6 consider partitioning schemes using R ) 1 and R ) 1/2 as parameters that determine the partitioning pattern. Because the choice of R ) 1 is often best for pooling problems,1,58 we have not considered a variety of partitioning schemes in this work. On the basis of the computational studies reported in Table 3, we observe that the improved hardware and software alone allows for significant improvements in computational efficiency. The CPU times needed in 2009 are, on average, 4.9% of the CPU times required in 2003. 4. Global Optimization Strategies This section discusses the branch-and-bound global optimization algorithm we used to solve the generalized pooling problem. At each node in the branch-and-bound tree, we minimize a linear relaxation of the node using the solver
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 60
61
CPLEX in the modeling language GAMS and split (branch) the node to create two child nodes. Then, we fix the current values of the binary variables, initialize the continuous variables using their current values, and locally minimize the resulting NLP using CONOPT.62 Unlike the topologies studied by Bergamini et al.,55 this topology apparently does not require bound contraction of the NLP to generate good upper bounds. If the local solve is less than the current upper bound, then the new upper bound replaces the old. At each step, nodes with relaxations within (ε ) 0.001) or larger than the current upper bound are considered, that is, nodes such that (LBNODE · (1 + ε) g UB) are eliminated. The algorithm terminates with ε convergence when there are no remaining nodes in the branch-and-bound tree. Refer to the textbooks of Floudas63,64 for comprehensive coverage of branch-and-bound algorithms. The following sections, which discuss the global optimization strategies we explored, include McCormick convex/concave envelopes,40,50 reformulation-linearization technique (RLT) relaxations,45 piecewise-linear underestimators,1,2 branching strategies, and bound tightening. 4.1. McCormick Convex/Concave Envelopes. To termwiseunderestimate the bilinear terms in the model, we replace each bilinear term with an auxiliary variable:
products of the qc,t variable bound factors and eqs 3, 4, and 6. Equations RLTQB-RLTQE represent three sets of four constraints each: [bt,e - ybt,e · b_t,e ; ybt,e · bjt,e - bt,e] × [qc,t - q_c,t ; qjc,t - qc,t] g 0 ∀c ∈ C; t ∈ T; e ∈ E (RLTQB) [ct,t' - yct,t' · _ct,t' ; yct,t' · jct,t' - ct,t'] × [qc,t - q_c,t ; qjc,t - qc,t] g 0 ∀c ∈ C; t ∈ T; t' ∈ T\{t} (RLTQC) [
∑b
t,e
e∈E
∑
+
qc wc,t,t'
{
∑ qc,ssource · [ fssource - ∑ as,e] - ∑ e∈E
s∈S
t∈T t'∈T\{t}
qb wc,t,e
{
gqc,t · bjt,e + qjc,t · (bt,e - bjt,e) gq_c,t · bt,e eqjc,t · bt,e eqc,t · bjt,e + q_c,t · (bt,e - bjt,e)
(21)
∀c ∈ C; t ∈ T; e ∈ E
(22)
In eqs 21 and 22, the semicontinuous variables ct,t′ and bt,e are bounded below by 0 rather than _ct,t′ and b_t,e to allow for the possibility that the given pipeline is deactivated. We term the MILP generated by replacing all bilinear terms with wqc c,t,t′ or wqb c,t,e the McCormick-based relaxation. 4.2. RLT Relaxations. The reformulation-linearization technique (RLT) developed by Sherali and co-workers (e.g., Sherali and Alameddine44 and Sherali and Adams45) augments an optimization problem with constraints that, while redundant in the full model, tighten the problem relaxation. Meyer and Floudas6 augmented their relaxation of the generalized pooling problem by appending to the formulation the products of qc,t variable bound factors (i.e., qc,t - qc,t and qjc,t - qc,t) and activation constraints analogous to eqs 2-6. The advantage of the resulting redundant constraints is that the nonlinear terms introduced by the additional equations are bilinear products of a continuous and binary variable (e.g., c qc,t · yt,t′ ) and hence are fully characterized by their convex envelope. Because the bilinear terms endogenous to our new model are (qc,t · ct,t′) and (qc,t · bt,e), we follow Meyer and Floudas6 in considering additional tightening constraints that are
t,e
-
e∈E
∑
ct,t'] ×
t'∈T\{t}
In their integrated water systems study, Karuppiah and Grossmann54 appended a constraint which balanced the contaminants. Similarly, we consider eq RLTBALwhich balances the change in contaminant flow rate between the sources and sinks (left-hand side) with the flow rate extracted from the treatment plants (right-hand side).
∑
∀c ∈ C; t ∈ T; t' ∈ T\{t}
∑b
[qc,t - q_c,t ; qjc,t - qc,t] g 0 ∀c ∈ C; t ∈ T (RLTQE)
wqb c,t,e
qc Each auxilary variable wc,t,t′ and wqb c,t,e is constrained by the convex envelope of the bilinear term:40
gqc,t · jct,t' + qjc,t · (ct,t' - jct,t') gq_c,t · ct,t' eqjc,t · ct,t' eqc,t · jct,t' + q_c,t · (ct,t' - jct,t')
ct,t' - _et · yte ; jet · yte -
t'∈T\{t}
ct,t' · qc,t r wqc c,t,t' bt,e · qc,t r
5429
qc,t · bt,e )
t∈T e∈E
rc,t ·q ·c 1 - rc,t c,t t,t'
∀c ∈ C (RLTBAL)
The advantage in introducing eq RLTQC into the termwise relaxation of this model is to tighten the McCormick-based relaxation of the small and midsize case studies from 7.219 × 105 to 7.223 × 105 (a 0.05% improvement). Equations RLTQB, RLTQE, and RLTBAL offer no additional augmentation to the McCormick-based relaxation for this formulation of the water treatment topology. Because the introduction of the piecewise underestimators described in the following section eliminates the small advantage of eq RLTQC, we choose not to increase the model size with any redundant constraints. Additionally, we reduce the solution time of each MILP node by choosing not to include piecewise versions of eqs RLTQB-RLTBAL as Meyer and Floudas6 propose. 4.3. Piecewise-Linear Underestimators. Following the implications of recent studies by Wicaksono and Karimi2 and Gounaris et al.,1 we choose to piecewise-underestimate each of the (qc,t · ct,t′) and (qc,t · bt,e) bilinear terms in the model using the nf4r formulation.1 In the five test cases we consider, there are fewer qc,t variables than bt,e and ct,t′ variables, so we expedite the solution time by partitioning the qc,t variables.1 We uniformly partition the qc,t variables because equal partitioning often produces the tightest relaxation.1 The following nf4r formulation is taken from the work of Gounaris et al.1 To construct the grid that will be used to build the piecewiselinear underestimators, each qc,t variable is ab initio partitioned into N segments: n qc,t(n) ) q_c,t + · (qjc,t - q_c,t) N
∀c, t, n ) 0, ..., N (23)
and an SOS1 variable λc,t(n) is introduced to activate one and only one domain segment: λc,t(n) )
{
1 0
if qc,t(n - 1) e qc,t e qc,t(n) else ∀c ∈ C; t ∈ T (24)
5430
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
N
∑
N
qc,t(n - 1) · λc,t(n) e qc,t e
n)1
∑q
c,t(n) · λc,t(n)
n)1
∀c ∈ C; t ∈ T; n ) 1, ... , N (25) N
∑λ
c,t(n)
)1
∀c ∈ C; t ∈ T
(26)
n)1
Continuous variables ∆cc,t,t′(n), ∆bc,t,e(n), n ) 1, ..., N are place holders for flow rates ct,t′ and bt,e that are nonzero in all but the active interval where they are equal to ct,t′ or bt,e. N
ct,t' )
∑ ∆c
c,t,t'(n)
∀t ∈ T; t' ∈ T\{t}
(27)
∀t ∈ T; e ∈ E
(28)
n)1
N
∑ ∆b
bt,e )
c,t,e(n)
n)1
0 e ∆cc,t,t'(n) e jct,t' · λc,t(n) ∀c ∈ C; t ∈ T; t' ∈ T\{t}; n ) 1, ... , N (29) 0 e ∆bc,t,e(n) e bjt,e · λc,t(n) ∀c ∈ C; t ∈ T; e ∈ E; n ) 1, ... , N. (30) Equations 27-30 could have been partitioned based on the difference between the upper and lower bounds of ct,t′ or bt,e, but we found no advantage for determining a nonzero lower bound for these variables because of the generalized pooling problem’s combinatorial nature (see section 4.5). The final relaxation of the bilinear terms ∀ c, t, t′, e is
{ {
N
g
qc wc,t,t'
∑q
c,t(n
N
gcjt,t' · qc,t +
e
g
c,t(n) · (∆cc,t,t'(n)
- jct,t' · λc,t(n))}
∑ {q
c,t(n
- 1) · (∆cc,t,t'(n) - jct,t' · λc,t(n))}
n)1
N
∑q
c,t(n) · ∆cc,t,t'(n)
n)1
∀c, t, t' (31)
arg max c∈C t∈T
N
∑q
c,t(n
∑
c qc yˆt,t' · |wˆc,t,t' - qˆc,t · cˆt,t' | +
t'∈T\{t}
∑ yˆ
b qb ˆ c,t,e t,e · |w
- qˆc,t · bˆt,e |
t'∈E
- 1) · ∆bc,t,e(n)
(33)
n)1
N
gbjt,e · qc,t +
∑ {q
c,t(n) · (∆bc,t,e(n)
- bjt,e · λc,t(n))}
n)1 N
ebjt,e · qc,t +
e
∑ {q n)1 N
ecjt,t' · qc,t +
qb wc,t,e
- 1) · ∆cc,t,t'(n)
n)1
tioning level is intuitive because the ideal number of discretizations varies between topologies and test cases (e.g., the pooling problem test cases of Foulds et al.39 were designed to converge with only a single McCormick-based relaxation, so extra partitions would delay solution time, but Pham et al.57 showed that the larger pooling problem test cases of Adhya et al.42 may benefit from partitioning). Because real-world applications need to be solved quickly once and engineers may not have the luxury of trying multiple discretization levels, we suggest that problems with similar topologies benefit from similar discretizations. In this work, we report several partitioning levels. After observing the effect of each partitioning level on the small and midsize test cases with four and ten plants, respectively, we are able to optimize the large fifteen plant test case as well as address the two very large twenty plant case studies to 0.9% and 2.3% optimality gaps or to global optimality under an additional assumption of a limited number of treatment facilities. 4.4. Branching Strategies. After solving a relaxation of each node using termwise McCormick convex/concave envelopes or one of the piecewise-linear underestimators described in section 4.3, we select (1) a variable for branching and (2) the branching point. This section outlines the algorithm we use to make both selections. First, we choose a variable set for branching. Because we use the MILP solver CPLEX60 to optimize the relaxation of each node, we only need to branch on variables participating in the nonconvex bilinear terms. The bilinear terms in the water treatment model are products of regulated qualities (qc,t) at the exit of a treatment plant and the flow rate of the streams leaving each treatment plant (bt,e and ct,t′), so each qc,t variable participates in (|T| - 1 + |E|) bilinear terms and each bt,e and ct,t′ variable participates in (|C|) terms. In all of the generalized pooling test cases shown in Table 2, (|T| - 1 + |E| > |C|), so we increase the effect of a single branching step by splitting the qc,t variables. After generating the optimal solution (yˆa, ..., yˆe, aˆ, ..., qˆ) of a given node, we follow common practice by branching on the variable qc,t where (c, t) contributes to the greatest discrepancy between the auxiliary and original problem variables:47,65,66
∑ {q
c,t(n
- 1) · (∆bc,t,e(n) - bjt,e · λc,t(n))}
After choosing the specific variable qc,t for branching, we select the point within the bounds [qc,t, qjc,t] at which to divide the single node into two nodes. The three options we considered were branching halfway between the bounds:
n)1
N
1 (q + qjc,t) 2 _c,t
∑q
c,t(n) · ∆bc,t,e(n)
n)1
∀c, t, e (32)
Our previous work indicates a trade-off between few partitions, which solve quickly, and many partitions, which may close the gap in a single node. But it is not clear how many partitions are appropriate for a given topology. Meyer and Floudas6 reported timing statistics for several partitioning levels, but that study applies only to the root node of a global optimization algorithm. Karuppiah and Grossmann,54 Bergamini et al.,55 Pham et al.,57 and Misener et al.49 chose the number of discretizations, presumably by running each example with several partitioning levels and reporting the discretization which generated the best results. Experimentally choosing the parti-
(BCHHALF)
branching at the optimal solution qˆc,t of the node and branching 10% away from the optimal solution qˆc,t:
{
1.1qˆc,t
branch
point ) 0.9q ˆ c,t qc,t
1 s.t. 1.1qˆc,t < (q_c,t + qjc,t) 2 1 s.t. 0.9qˆc,t > (q_c,t + qjc,t) (BCHTEN) 2
1 (q + qjc,t) else 2 _c,t
Note that branching scheme (BCHTEN) assumes that qc,t g 0, a legitimate assumption for regulated qualities in a wastewater treatment application. Branching at the halfway point was the
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5431
Figure 2. Optimal topologies for the five generalized pooling test cases.
least advantageous for the water treatment problem. For example, the algorithm using eq BCHHALF required 59 nodes in the branch-and-bound tree to solve the small four-plant test case when the partitioning level was N ) 5 in comparison to the 7 nodes needed by the other two branching schemes. In our trials, branching scheme (BCHTEN) always performed at least as well as branching at the optimal solution qˆc,t and produced significantly better results at smaller partitioning levels where there were many nodes in the branch-and-bound tree (e.g., solving the small four-plant test case with N ) 3 required 43 nodes when branching ten percent away from the optimal solution, 62 nodes when branching at the optimal solution, and 102 nodes when branching at the variable halfway point). On the basis of these trials of the four-plant test case, we exclusively used branching scheme (BCHTEN) for the remaining case studies. Discussion on branching selection can be found in the books of Floudas63 and Tawarmalani and Sahinidis46 and the article of Adjiman et al.66
4.5. Bounds Tightening. As we discussed in section 3.4, tight bounds on variables participating in bilinear terms are critically important.59 Therefore, bounds on the bt,e and ct,t′ flow rates are updated according to topological analysis on the quality bounds at each node in the branch-and-bound tree:
bjt,e ) min c∈C
jct,t' ) min c∈C
(
(
qmax c,e ·
∑f
source s s∈S max[qmax _c,t] c,e , q
qjc,t' ·
)
∑f
source s
s∈S
max[qjc,t', (1 - rc,t') · q_c,t]
∀t ∈ T; e ∈ E
)
(34)
∀t ∈ T; t' ∈ T\{t}
(35)
To expedite convergence, we also use topological reasoning to eliminate binary variables at nodes in the branch-and-bound tree
5432
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Table 4. Branch and Bound Results for the 4-Plant Test Case (Serial Optimization) root node
termination
CPU time (s)
N
relaxation
CPU (s)
no. nodes
lower bnd
upper bnd
% gap
solving
total
McC RLT 3 5 7 10 15
7.219 × 105 7.223 × 105 9.607 × 105 9.896 × 105 1.005 × 106 1.041 × 106 1.086 × 106
0.1 0.6 2.9 5.9 5.6 16.6 109.3
5000a 1712 95 7 5 3 1
1.082 × 106 1.085 × 106 1.085 × 106 1.085 × 106 1.085 × 106 1.085 × 106 1.086 × 106
1.106 × 106 1.086 × 106 1.086 × 106 1.086 × 106 1.086 × 106 1.086 × 106 1.086 × 106
2.2 0.1 0.1 0.1 0.1 0.1
