Global Optimal Structures of Heat Exchanger Networks by Piecewise

1752

Ind. Eng. Chem. Res. 2007, 46, 1752-1763

PROCESS DESIGN AND CONTROL Global Optimal Structures of Heat Exchanger Networks by Piecewise Relaxation Maria L. Bergamini, Nicola´ s J. Scenna, and Pio A. Aguirre* INGARsInstituto de Desarrollo y Disen˜ o, AVellaneda 3657, S3002GJC, Santa Fe, Argentina

A new methodology for the global optimization of heat exchanger networks is presented, based on an outer approximation methodology, aided by physical insights. The problem is formulated as a mixed-integer nonlinear problem (MINLP). Two lower bounding convex MINLP problems are constructed, including piecewise underestimators of the nonconvex terms. A solution of the bounding problems gives an approximated optimal solution that is used as an initial point for solving the reduced NLP problem, or that proves there is no better solution. Meanwhile, the global bounding problem selects feasible structures with improved objective value. In order to reduce the number of feasible structures to be explored, rigorous constraints obtained from physical insights are included in the bounding problems. Networks with up to 9 process streams have been solved to global optimality, improving considerably the computing time required to solve them. Networks with more than 10 process streams have been solved with the global optimization strategy, giving as results a set of feasible structures, with total cost near the optimal value. This offers the possibility of alternatives to be analyzed, considering other aspects not reflected in the original MINLP model. 1. Introduction Heat exchanger network (HEN) synthesis is one of the most studied problems in chemical engineering. The mathematical programming approach to this problem, considering the simultaneous optimization of energy consumption, streams matches, and heat recovery areas, requires the solution of a mixed-integer nonlinear problem (MINLP). Usual techniques for the solution of MINLP models include branch-and-cut,1 generalized benders decomposition2 (GBD), outer approximation3,4 (OA), and the extended cutting planes.5 It is well-known that these techniques are well-suited for convex problems, but they might get trapped at suboptimal solution or even fail to obtain a feasible solution for nonconvex problems. Thermodynamically based convex underestimators were developed by Zamora and Grossmann6 and implemented in their global optimization branch-and-bound algorithm to ensure global convergence and global optimality on the HEN model. Adjiman et al.7 have also implemented the SMIN-RBB algorithm to solve the HEN synthesis problem to global optimality, assuming linear cost function for the areas. It is also worth mentioning the approach by Bjo¨rk and Westerlund.8 Their work attempts to convexify the signomial terms and to create at fpproximated convex subproblems. The HEN synthesis problem is as follows. Given are I hot process streams to be cooled and J cold process streams to be heated, with known heat capacity flowrates and known inlet and outlet temperatures. Given are also auxiliary heating and cooling utilities with their respective temperatures and their cost per unit of heat provided or removed. The problem deals with the simultaneous determination of energy consumption, number of exchangers, area requirement, and network configuration, * To whom all correspondence should be addressed. Tel.: (54-342)4534451. Fax: (54-342)-4553439. E-mail: [email protected].

which minimize annual investment cost plus the annual cost of utilities. Yee and Grossmann9 proposed a superstructure and a MINLP model for dealing with the HEN synthesis problem that takes into account appropriately the tradeoff between different costs (utility, investment, and operating costs). The proposed superstructure is fairly general, allowing both series and parallel decoupling of the heat exchangers, because of the stagewise superstructure. The superstructure for a two-cold-two-hot set of streams is shown in Figure 1. Generally, the number of stages is set equal to the maximum cardinality of the hot and cold sets of streams (although sometimes it is necessary to increase the number of stages to allow designs with minimum energy consumption). At each stage, hot and cold streams are split to allow the potential existence of a heat exchanger to match any hot-cold pair of streams. Auxiliary cooling and heating utilities are placed at the outlets of the superstructure. The model that arises with this superstructure is a nonconvex MINLP problem. In this work, a new solution strategy based on the outer approximation for global optimization (OAGO) algorithm10 is presented and applied to solve the HEN synthesis problem to global optimality. It does not rely on a branch-and-bound framework. Instead of that, the algorithm solves underestimating subproblems whose feasible solutions are approximated solutions with relaxed objective values lower than the best-known solution. The mentioned algorithm has been successfully applied to several generalized disjunctive programming (GDP) models that arise in general process synthesis problems, including process networks, wastewater treatment, reaction path, and separation network. The technique can also solve NLP problems to global optimality. In the present work, the algorithm is specialized to solve the HEN synthesis. The proposed methodology not only generates the global optimal structure and its operating conditions but also provides

10.1021/ie061288p CCC: $37.00 © 2007 American Chemical Society Published on Web 02/09/2007

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007 1753

the required utility, and nonlinear cost for the exchanger area. It is assumed that the cost of the area is of the form f(x) ) ∑ica‚areaβi , where ca is the cost coefficient, β is a scale factor, and the summation index i runs for all the units. The m binary variables y are linearly constrained. The set of linear constraints involves energy balances, temperature difference definition, and constraints for the monotonic decrease of temperature. The nonlinear constraint hi(x) is bilinear and represents the area definition for the ith unit,

hi(x) ) areai‚LMTDi -

Figure 1. HEN superstructure proposed by Yee and Grossmann.

the optimal operating conditions of alternative structures. The algorithm selects a potentially optimal structure in an outer phase and then searches for the optimal values of the continuous variables in this structure. The structures selected before the optimal one are also globally optimized in the inner phase. Therefore, the information provided in these steps allows subsequent consideration of alternative structures, taking into account aspects not considered in the model (safety, operability, controllability, reliability, etc.). The algorithm is also capable of dealing with a less ambitious problem, i.e., the global convergence of a nonlinear problem. The methodology can always obtain a local optimal solution of a nonconvex problem, without requiring an externally supplied initial point. 2. Mathematical MINLP Model The methodology refers to the model by Yee and Grossmann,9 with the isothermal mixing assumption. It leads to a model almost linearly constrained (since it does not consider energy balance at the mixers), with nonconvexities in the equality constraints defining the required areas for each exchanger and the area costs. The inlet and outlet temperatures of the hot and cold streams are known. The heat capacity flowrates are constant values and also known. The model, denoted P(x,y), involves binary variables y representing the existence of exchangers and continuous variables x denoting stream temperatures, approach temperature differences, exchanged heat, and exchanger areas. The model P(x,y) can be outlined as follows (for a detailed description of the model, see Appendix A),

min TAC ) cf‚y + c‚x + f(x)

gi(x) e 0

where Ui is the heat transfer coefficient and qi and LMTDi denote the heat exchanged in unit i and the logarithmic mean temperature difference of the exchange, respectively. The nonlinear constraints gi(x) are convex and estimate the logarithmic mean temperature difference by the Paterson approximation,11

gi(x) ) LMTDi -

[61 (∆T + ∆T ) + 32 x∆T ∆T ] 1 i

2 i

1 i

2 i

(2)

In eq 2, ∆T1i and ∆T2i denote the temperature approach at both extremes of the unit i. Note that these constraints are convex. It is possible to express the LMTDi definition as an inequality since the objective function forces the variable LMTDi to increase. These inequalities will be active in an optimal solution. To show that, recall that NLP solvers search for a point that satisfies the optimality conditions. They say that, in an optimal point, the gradient of the objective function is a linear combination of the gradients of the equality constraints plus a linear negative combination of the active constraints of the form e. In our model, the objective function reflects the total annual cost, TAC. The component of the gradient corresponding to the variable LMTDi, for an existent exchanger i (areai * 0), is

∂(areai - qi/UiLMTDi) ∂TAC ) λi + ∂LMTDi ∂LMTDi µi

∂[LMTDi - (1/6(∆T1i + ∆T2i ) + 2/3x∆T1i .∆T2i )] ∂LMTDi

where µi e 0, and µi < 0 if and only if eq 2 is active. Calculating the derivatives, the following is obtained:

qi 0 ) λi + µi UiLMTDi2 Moreover, the component of the gradient corresponding to areai is as follows:

Calculating the derivatives,

i ) 1, ..., m i ) 1, ..., m

(1)

∂(areai - qi/UiLMTDi) ∂TAC ) λi ∂areai ∂areai

s.t. Ax + By e b hi(x) ) 0

qi Ui

ca‚β‚areaβ-1 ) λi i P(x,y)

x∈Ω y ∈ {0, 1}m The objective function reflects the total annual cost. It considers fixed costs for existing exchangers, linear cost for

Then, λi > 0, and also qi > 0 (since areai > 0). This shows that µi < 0. In conclusion, the inequality defining LMTDi is active in an optimal solution. Finally, Ω is a set defined by the bound of the continuous variables. Some authors avoid the isothermal mixing assumption, requiring that each stream participates in at most one exchanger

1754

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007

in each stage. This is known as the no stream splitting assumption, and it is modeled by using constraints on the binary variables representing the existence of units. Of course, the resulting model is much restricted. This assumption is implemented in some examples in Section 7. Removing the isothermal mixing assumption makes the model harder to solve, introducing bilinear terms in constraints representing energy balances in mixing points. Although the model thus generated is more general and could have optimal solutions with lower cost, in practice the improvement is rarely greater than 5%. Therefore, the general model is not considered here. The described model is a nonconvex MINLP. Both the objective function and the feasible region are nonconvex, which can cause the local solver to miss the global solution or even to be unable to find a feasible solution. The algorithm proposed here is based on piecewise relaxation that replaces the nonconvex terms, so that a rigorous bounding problem is generated. This problem has a solution, as long as the nonconvex problem P(x,y) has an improved optimal solution, and this property is used as the stop criterion in the algorithm. 3. Piecewise Relaxation Convex MINLP underestimating subproblems are set up, replacing the bilinear constraints of the form qi/Ui ) areai‚ LMTDi and the concave terms like areaβi in the objective function by piecewise relaxations. The relaxations are constructed over a partition of the domain for areai. Piecewise linear function have been used widely for approximating separable functions.12 With the significant progress in linear (continuous, integer, or mixed) optimization codes, piecewise approximations have become a more attractive alternative for handling nonlinear functions, even at the cost of increasing the size of the model in the continuous and integer domains.13 Several formulations of piecewise functions have been studied, both considering binary and continuous variables and including only continuous variables.14 up Consider a partition K of the interval [arealo i , areai ], con* * sisting on N + 1 grid point, K ) {areai 1, areai 2, ..., areai*N+1}, up * where areai*1 ) arealo i and areai N+1 ) areai . The piecewise linear formulation requires the addition of new continuous variables δin that denote the involvement of the / ] in the current value of areai. In subinterval [area/in, areain+1 addition, binary variables win are necessary for denoting the subinterval where the value of areai belongs. Consider the concave term CAi ) areaβi . The piecewise linear relaxation is formulated with the following restrictions, N

/ areai ) areai1 +

∑ δin

C ˜ Ai )

+

∑

/ (areain+1 )β - (area/in)β / areain+1 - area/in

As is known, piecewise linear relaxation of concave terms holds the following properties: •C ˜ Ai ) CAi ) areaβi when areai ) areai*n for some n ) 1, ..., N + 1. • C ˜ Ai e CAi: the piecewise linear relaxation is a rigorous lower bound for the concave function (this is not true for nonconcave functions). Moreover, the approximation gap (CAi - C ˜ Ai) decreases when the size of the subintervals [area/in, / areain+1 ] is reduced. Consider now the bilinear term relating area, logarithmic mean temperature difference, and exchanged heat. Multiplying eq 3 by LMTDi and using eq 1,

qi Ui

N

/ ) areai1 LMTDi +

∑ γin

(6)

n)1

where the new continuous variables γin represent the bilinear term δinLMTDi. In order to avoid bilinear terms δinLMTDi ) γin, the following linear cuts defining the convex and concave envelopes are added,

}

lo lo lo γin g δinLMTDlo i + δinLMTDi - δinLMTDi up up up γin g δinLMTDup i + δin LMTDi - δin LMTDi n ) 1, ..., N lo lo up γin e δinLMTDup i + δinLMTDi - δinLMTDi up up lo γin e δinLMTDlo i + δin LMTDi - δin LMTDi (7)

Equation 7 defines the convex and concave envelope of the up lo up bilinear term in the subregion [δlo in, δin ] × [LMTDi , LMTDi ]. As is well-known, these inequalities make γin exact at the boundaries of the region, and the relaxation gap is proportional lo up lo to the area of the subregion, (δup in - δin)(LMTDi - LMTDi ). The lower and upper bound for the new continuous variables up / / δin are δlo in ) 0 and δin ) areain+1 - areain. Note that, if win ) / / / 1, then δin ) areain+1 - areain and γin ) (areain+1 / areain)LMTDi. Moreover, if win ) 1, then wim ) 1 and γim ) / (areaim+1 - area/im)LMTDi for all m < n. On the other hand, if win-1 ) 0, δin ) 0 and also γin ) 0. Moreover, wim ) 0 for all m > n - 1 and γim ) 0. Thus, the bilinear variable γin is exact if either win ) 1 or win-1 ) 0, and there is an approximation gap only in one subinterval.

(3)

n)1

/ β (areai1 )

From the previous restrictions, it follows that 1 g wi1 g wi2 g ... g wi N-1 g 0.

4. OAGO Algorithm

δin

(4)

/ / / / (areai2 - areai1 )wi1 e δi1 e areai2 - areai1

The OAGO algorithm consists of two phases: an outer optimization and an inner optimization. In the outer optimization, a potentially optimal structure is selected. This structure is globally optimized in the inner phase, after performing a bound contraction on the complicating variables, as is explained below.

/ / (areain+1 - area/in)win e δin e (areain+1 - area/in)win-1 n ) 2, ..., N - 1 (5)

4.1. Subproblems Definition. The algorithm relies on three subproblems:

/ δiN e (areaiN+1 - area/iN)wiN-1

GR(x,y,w)K. This is a convex MINLP subproblem that globally relaxes the problem P(x,y). It has the form,


min TAC ) cf‚y + c‚x + bψ s.t.TAC e TACup Ax + By e b DK(x,ψ,dx,w,y) ) 0 g(x) e 0

GR(x,y,w)K

Φ(y) ) φ x ∈ Ω, dx ∈ Rn1K y ∈ {0,1}m, w ∈ {0,1}n2K This problem has binary variables y related to the existence of the exchangers and binary variables w related to piecewise underestimators. DK(x,ψ,dx,w,y) represents the formulation of the piecewise relaxation, where dx are the new continuous variables. The superscript K means that this relaxation is constructed over a specific partition K. The set of constraints Φ(y) ) φ defines integer cuts. An upper bound to the total annual cost is imposed, TACup ) (1 - )TACopt, where TACopt is the best-known solution (or infinite, while no solution is at hand) and is a relative optimality tolerance. Because of this bound constraint, the infeasibility of this problem indicates that the global solution of the original nonconvex problem P(x,y) has been found (the best-known solution) or that P(x,y) is infeasible (if no solution has been found). P(x,yl). This is a reduced nonconvex NLP subproblem representing a fixed structure in the network. This problem is obtained from P(x,y) by fixing the structural binary variables y in values yl. The variables in this problem are continuous and correspond to stream temperatures, approach temperature differences, exchanged heat, and exchanger areas for existing units. LR(x,yl,w)K. It is a convex local underestimating MINLP subproblem, obtained from GR(x,y,w)K by fixing the structural variables y (and, therefore, without considering the integer cuts). This problem is also a relaxation of P(x,yl). It has binary variables w defining the piecewise underestimators. Infeasibility of this problem implies that either the P(x,yl) problem corresponding to the set of values yl is infeasible or that the global solution of the current structure is not better than the best-known solution. The algorithm requires solving the underestimating problems GR(x,y,w)K and LR(x,yl,w)K to feasibility, i.e., a feasible solution of the subproblems is sufficient for the algorithm goal. 4.2. Procedures. 4.2.a. Outer Approximation for Feasibility (OAF) Procedure. An ad hoc procedure was implemented for solving the GR(x,y,w)K to feasibility, modifying the outer approximation algorithm.3 The procedure relies on two problems: a linear mixed-integer problem GL(x,y,w)K,L and a convex nonlinear subproblem GR(x,yo,wo)K. The linear problem is obtained by replacing the convex constraints of the form g(x) e 0 by valid linearizations at points xi. The convex problem GR(x,yo,wo)K is obtained by fixing the binary variables at the values yo and wo obtained as solutions of the GL(x,y,w)K,L subproblem. Given a gridpoint set K, and an upper bound for the TAC, the procedure is as follows. (1) Initialization: Initialize the set L of linearization points. (2) Solve MILP: Solve GL(x,y,w)K,L with a MILP solver to get a feasible solution. If a feasible solution is found, fix the values of the binary variables y and w and go to Solve NLP. Otherwise, go to Stop Infeasible.

Figure 2. Structure with several binary variables representation.

Figure 3. Forbidden structure by the integer cut when stream splitting occurs.

(3) Solve NLP: Solve GR(x,yo,wo)K with a NLP solver. If it is infeasible, add the values of x to L, construct GL(x,y,w)K,L by adding an integer cut to avoid repeating the values for the binary variables, and go to Solve MILP. Otherwise, go to Stop feasible. (4) Stop feasible: A new feasible structure with better TAC was found. Exit. (5) Stop infeasible: Report infeasibility for GR(x,y,w)K. Exit. An equivalent procedure is applied for solving the local relaxed problem LR(x,yl,w)K, holding the structural binary variables fixed. 4.2.b. Bound Contraction Procedure. This procedure eliminates parts of the feasible region where there are no feasible points or the objective function is worse than the best solution found. Moreover, since the relaxation of nonconvex terms depends on the variable bounds, tight bounds adjust the relaxation. The bound contraction procedure, based on the work of Zamora and Grossmann,15 is performed after a feasible structure is fixed; it calculates tight bounds valid for the current structure. The procedure relies on solving the problems

min/max xi s.t. constraints in LR(x,yl,w)K

(8)

for each complicating variable xi. In the HEN synthesis problem, the complicating variables are those involved in nonconvex terms (area and temperature difference for each potential exchanger). Although the structural binary variables are fixed, the bound contraction problems may involve binary variables related to the piecewise estimators. However, using a grid with two points up ({arealo i , areai }), these are LP subproblems. 4.2.c. Grid Update Procedure. In order to improve the approximation of the piecewise relaxation of bilinear and concave terms, the grid is refined in the interval where the solution of the last bounding problems GR(x,y,w)K or LR(x,yl,w)K lie, when the relaxation gap is greater than a certain tolerance. Two strategies for updating the grid are used.

1756


(4) Inner optimization: Find the global solution of P(x,yl) or detect if the current structure is suboptimal, through the following procedure: (4.a) Exact solution: solve P(x,yl) with a local optimizer. If a feasible solution is found, denote TAC* the optimal objective value and y* the selected structure. If TACup > TAC*, set TACup ) (1 - )TAC* and yopt) y*. Go to the Lower approximation step. (4.b) Lower approximation: Update the gridpoint set K. Solve the problem LR(x,yl,w)K to feasibility with the OAF procedure. If it is feasible, go to the Exact solution step. Otherwise go to the New structure step. (5) New structure: Generate and add integer cuts to make the examined structure infeasible. Unfix the binary variables y and go to the Outer optimization step. (6) Stop: The global optimal solution is TACopt ) TACup/(1 - ), and the optimal configuration is yopt. Figure 4. Amin vs Q bounding curve with solution points from P(x,yl) (dots) and GR(x,y,w)K (stars) without imposing the bounding curve, corresponding to example 3.

5. Integer Cuts As was detailed in the previous section, the OAGO algorithm implements a search in two phases. An outer optimization phase selects the network structure, and the inner optimization phase globally optimizes the selected structure. In each outer optimization, a new structure is desirable, and then integer cuts are added in order to avoid the selection of the same structure in subsequent iterations, as it is implemented in the classical outerapproximation algorithm by Duran and Grossmann.3 If the analyzed structure is defined by the binary values yli, the added cut in the outer approximation algorithm is the following,

∑

i:yli ) 1

Figure 5. Amin vs Q bounding curve and solution points of problems P(x,yl) (dots) and (stars) imposing the curve as a restriction in the problem GR(x,y,w)K; example 3.

(1) Strategy A: add to the gridpoint set the middle point of the active subinterval. (2) Strategy B: add to the gridpoint set the solution point of the bounding problem 4.3. Algorithm. The steps of the OAGO algorithm are detailed below: (1) Initialization: Set the upper bound TACup ) ∞ or a finite value if it is known (obtained as a heuristic solution of P(x,y)); set global convergence tolerance and the initial gridpoint set K. (2) Outer optimization: Set global bounds for the complicating variables. Find a feasible solution of the problem GR(x,y,w)K with the OAF procedure. If the procedure reports infeasibility, go to the Stop step. Otherwise denote yl the values of the binary variables, fix y ) yl, update the grid, and go to the Bound contraction step. (3) Bound contraction: Apply Bound contraction procedure for the variables areai, ∆T1i , and ∆T2i , associated to existing units i, to update their bounds and the bounds on LMTDi. Fix to zero the variables related to nonexistent units. If one of these problems is infeasible, go to the New structure step.

yi -

∑

i:yli ) 0

yi e

∑i yli - 1

(9)

In the HEN problem, different binary vectors may frequently represent equivalent structures. This is due to the staged location of the exchangers. Then, in our algorithm, it is necessary to include cuts for all the equivalent structures. An attempt was made to model the equivalent structures. However, this presents the drawback of requiring a combinatorial number of additional variables and constraints, making the bounding problems considerably hard to solve. Therefore, other strategies that prevent exploring equivalent structures were studied. The first one comes up from the fact that the location of a match in a particular stage is fairly free. For example, in the network depicted in Figure 2, the two existent units are dashed. The structure can be represented by assigning the binary variables values as y6 ) 1 and y7 ) 1, with all the others equal to zero. However, the same configuration is obtained with y3 ) 1 and y6 ) 1, or y7 ) 1 and y10 ) 1, or y3 ) 1 and y10 ) 1. In order to prevent this, integer cuts forcing a match to be placed most on the left in the superstructure are added. The cut can be expressed in logical language as follows: if an exchanger matching hot stream i and cold stream j is placed in stage k, then there is another exchanger in the stage k - 1, involving either i or j, or both. It should be mentioned that these cuts are valid only with the assumption of no stream splitting, since in other cases, they forbid feasible structures, as is shown in the example in Figure 3. Another approach of cutting all the y-representation of a structure is to generate at the end of each outer iteration, among

Ind. Eng. Chem. Res., Vol. 46, No. 6, 2007 1757 Table 1. Data of Examples

stream

Tin (°C)

Tout (°C)

Table 5. Solution Time for Example 1, with Strategy A

heat capacity flowrate (kW/°C)

film transfer coefficient (kW//(°C m2))

CPU time (s) cost ($/(kW year))

Example 1 CU HU H1 H2 H3 C1 C2

20 300 159 267 343 26 118

60 300 77 80 90 127 265

Cu HU H1 H2 C1 C2 C3 C4

25 325 180 240 40 120 40 80

40 325 75 60 230 260 130 190

CU HU H1 H2 H3 H4 C1 C2 C3 C4 C5

15 330 327 220 220 160 100 35 85 60 140

30 250 40 160 60 45 300 164 138 170 300

2.285 0.204 0.538 0.933 1.961

0.20 0.05 0.10 0.04 0.50 0.01 0.50

10 110

1.0 2.0 2.0 2.0 1.5 1.5 2.0 2.0

20 120

0.5 0.5 0.5 0.4 0.14 0.60 0.35 0.70 0.50 0.14 0.60

6 60

38 236 160 249 271 226 199 82 93 38 60 116

82 236 93 138 149 66 66 177 205 221 160 206

GR(x,yl,w)K

contraction

total

15

51.22

0.82

1.98

41.34

95.36

iteration

binary variables

continuous variables

equations

1 5 10 15

64 91 105 119

317 398 440 482

578 775 895 1017

exchanger

area

H1-CU H2-CU H3-C2 H3-C1 HU-C2

36.942 9.634 0.963 99.728 65.626

the topology of the network. The topology is preserved if the change does not superimpose units, neither add nor eliminate splits, and does not change the order in which the streams match each other. If it is possible, a cut for the resulting structure is generated and added to the global bounding problem. 6. Relaxation Tightening by Additional Constraints

Example 4 CU HU H1 H2 H3 H4 H5 C1 C2 C3 C4 C5

P(x,yl)

Table 7. Areas for Example 1, Calculated with the Paterson Approximation

Example 3 100 160 60 400 100 70 350 60 200

GR(x,y,w)K

Table 6. Global Bounding Problem Size (Strategy A)

Example 2 30 40 20 15 25 20

outer iterations

18.12 37.64 8.79 10.55 12.56 14.77 17.73 17.28 13.9 8.44 7.62 6.08

Table 2. Optimal Values in Example 1 exchanger

∆T1

∆T2

area (m2)

Q (kW)

H1-CU H2-CU H3-C2 H3-C1 HU-C2

99.00 207.00 203.64 138.15 35.00

57. 00 60.00 147.15 64.00 160.64

36.033 8.57 0.955 95.093 55.412

187.37 38.14 41.88 94.23 246.38

Table 3. Solution Time for Example 1, with Strategy B CPU time (s) outer iterations

GR(x,y,w)K

P(x,yl)

LR(x,yl,w)K

contraction

total

20

53.20

1.01

2.76

60.41

117.38

Table 4. Global Bounding Problem Size (Strategy B) iteration

binary variables

continuous variables

equations

1 5 10 15 20

64 102 118 137 153

317 431 478 536 584

578 853 988 1133 1261

the classical cuts (eq 9), other valid cuts for the equivalent structures. At the end of an iteration, the procedure inspects the optimized structure, looking for the possibility of changing the placement of an existent unit in other stages, without altering

The hardest step in the proposed algorithm is the solution of MILP subproblems. In order to help the MILP solver, additional constraints can be included in GR(x,y,w)K and LR(x,yl,w)K such that the continuous relaxation becomes tighter. These constraints are based on thermodynamic or graph theory knowledge, and they are, in general, redundant to the original model P. 6.1. Area Targeting. From pinch theory, it is possible to calculate the thermodynamical minimum utility requirement, given ∆Tmin. Area targeting methods predict the minimum area for heat transfer, given the minimum energy recovery.16 The area is calculated assuming equal heat transfer coefficients for all streams and utilities. However, according to Linnhoff and Ahmad,17 the deviation of the target area calculated with this assumption is not more than 10% of the entire area, even when the heat transfer coefficients differ by 1 order of magnitude. The minimum area calculated with the maximum heat transfer coefficient is a rigorous lower bound. The area target depends on the minimum utility consumption, or the ∆Tmin. Since, in the model of Yee and Grossmann, ∆Tmin does not have to be specified, any utility consumption may be feasible. Thus, the minimum area can be calculated as a function of utility consumption, resulting in a convex function. Since an analytical expression for this relation is not achievable, it is calculated for discrete values of utility consumption. This information is included in the relaxed problems through supporting linear functions in several points. Figure 4 shows the bounding curve Amin vs Q calculated for example 3, where Q is the cooling utility and Amin is the minimal heat recovering area. Figure 4 also shows the values of cooling energy and total heat recovery area for several structures, as well as the corresponding values taken from the bounding problem solutions, when the bounding curve is not imposed in GR(x,y,w)K. As can be seen, several bounding solutions lie below the curve, with infeasible values of Q and total area. Figure 5, as its own, shows the values corresponding to

1758


Figure 6. Optimal solution of example 2.

bounding problem solutions where the bounding curve was indeed imposed. Since the values of area are present in the objective function to be minimized and these variables are relaxed in GR(x,y,w)K and LR(x,yl,w)K, bounding these variables from below is crucial in getting good approximated solutions. 6.2. Overall Minimum Number of Units. The minimum overall units target is given by Nmin ) Ns - 1, where Ns is the number of both process and utility streams.16 The restriction bounding the total number of units with this value is redundant in model P. However, it is probably not redundant in the relaxed problem. Additionally, this restriction makes the continuous relaxation of the MILP problems tighter, and thus, the MILP solver is helped in finding a feasible integer solution. 6.3. Minimum Exchanger Area for Each Stream. A lower bound for exchange area for each stream can be derived considering the best case, which is the case when every stream transfers its heat load at maximum ∆T, with maximum film heat transfer coefficient. This nonzero lower bound allows disaggregating the minimum area calculated in Section 6.1 between all streams, avoiding the area accumulation in a few streams. 7. Examples The algorithm was implemented in GAMS, CPLEX 9.0 was used for solving the MILP problems, and CONOPT 3 was used for solving the NLP subproblems. Several examples were solved, and this section reports the results of four of them. A Pentium 4 PC was used, with a 1.5 GHz processor and 256 kB RAM memory. Previous works in global optimization of HEN synthesis problems consider up to 5 process streams in the network, or 6 streams in the case of the work of Bjo¨rk and Westerlund.8 Because of the combinatory complexity of the model, industrialsize problems are unmanageable within the current state of the art of global optimization. Table 1 shows the data for the problems reported in this section. The optimality tolerance used is ) 1%, unless another value is explicitly mentioned. 7.1. Example 1. The first example refers to a process with 3 hot streams and 2 cold streams. The superstructure consists of 3 stages, with 18 potential exchangers between process streams, 3 coolers, and 2 heaters. The area cost is 7400 + 80‚area (then,

Figure 7. T vs ∆H diagram of the optimal solution of example 2.

there are no concave terms in the objective function). The MINLP model P(x,y) has 23 binary variables and 103 continuous variables. In order to compare the results reported by other authors, in the first place, the example was solved by assuming the arithmetic mean temperature difference as an approximation to the logarithmic mean temperature difference, instead of eq 2. In this case, the bounding problems, GR(x,y,w)K and LR(x,yl,w)K, are linear. The optimal solution is a network with 2 exchangers between process streams, 2 coolers, and 1 heater, with an annual cost of $82 042.35. The units and their areas and exchanged heat are reported in Table 2. Zamora and Grossmann6 find the same solution and report a computation time of ∼6 h. Adjiman et al.7 solve the same problem with their branch-and-bound algorithm, requiring 577 CPU s for convergence to the global optimal solution. The example was solved with the proposed algorithm, using strategy B for updating the grid. The initial grid for each variable areai consists of 2 subintervals. For this example, the Amin vs Q curve was not included, since this is a small problem and does not present convergence difficulties. The first outer iteration gives as the solution a network with optimal cost $85 479.20, found and globally optimized in the


Figure 8. Optimal solution of example 3.

inner optimization phase. The second outer iteration finds a network with larger cost, and the global optimal structure is found in the third iteration. Following this, the algorithm performs 17 outer iterations, exploring other structures that are suboptimal. The 20th global bounding problem, including 60 integer cuts (representing the 19 explored structures, and other equivalents), is infeasible, and the algorithm stops. Table 3 shows the consumed time for solving each kind of problem. Note that the solution time is mostly consumed in the solution of the global bounding problem and the bound contraction procedure. In almost every iteration after finding the global optimal solution (third iteration), the algorithm does not enter to the inner phase, since the bound contraction procedure detects that the current structure is suboptimal (one of the bound contraction problems is infeasible). Table 4 shows the size of the global bounding problems GR(x,y,w)K in some representative iterations. If strategy A is used for updating the grid, the optimal solution is found in the 6th iteration, but the algorithm stops in the 15th iteration. Tables 5 and 6 show the solution time and the size of the problems, respectively. When the Paterson approximation (eq 2) is used for the logarithmic mean temperature difference, the bounding problems are convex, and they are solved with the OAF procedure. With this model, the optimal solution corresponds to the same structure, but the total annual cost is $83 389.13, since the exchanger areas are larger (see Table 7 for the optimal area values). The algorithm performance is similar in this case. After exploring 15 structures, the search stops (15 global bounding problems) and the optimal structure is found in iteration 9. The OAF procedure solves the 15 global bounding problems in 89.97 CPU s (83.42 s for the problems GL(x,y,w)K,L and 3.55 s for the convex problems GR(x,yo,wo)K). As in the previous case, the total times for solving the inner-phase problems are ∼1 CPU s. 7.2. Example 2. This example corresponds to a process with 6 streams (2 hot and 4 cold), taken from Bjo¨rk and Westerlund.8 The area cost is 8000 + 50‚area0.85. The cited authors propose a superstructure with 3 stages, allowing stream splits. The same conditions are specified in order to compare the results. The MINLP problem P(x,y) has 30 binary variables and 207 continuous variables.

Figure 9. T vs ∆H diagram of the optimal solution of example 3.

Figure 10. Progress of the best-known solution in example 3.

First, the algorithm was implemented for this example with strategy A for grid updating, with the cut generation procedure for equivalent structures at each outer iteration and the curve

1760


Figure 11. Best solution found for example 4.

Amin vs Q. Since, in this example, the model allows stream splits, the cut proposed in Section 5 is not valid. In the fist outer iteration, a network with a total annual cost of $143 602.72 is obtained. The global optimal (within the tolerance of 1%) configuration was obtained in the second outer iteration, with a cost of $ 140 367.07. The optimal values are shown in Figure 6, and Figure 7 shows the T vs ∆H diagram for all the exchangers in the optimal solution, between process streams. The algorithm explores one more structure in the third iteration, and in the fourth one, the global bounding problem is infeasible. The first GR(x,y,w)K has 30 binary variables, 306 continuous variables, and 662 constraints. The fourth (and last) problem GR(x,y,w)K has 128 binary variables, 574 continuous variables, and 1257 constraints. The four global bounding problems consume 146.89 CPU s (145.77 s for solving the MILP problems GL(x,w,y),K,L and 1.12 s for solving the convex NLP problems GR(x,yo,wo)K). Moreover, 14.74 s are consumed in the local bounding problem LR(x,yl,w)K, and 25 s are consumed in the bound contraction procedure. The algorithm by Bjo¨rk and Westerlund8 finds the same optimal cost, and the reported computational time is ∼8000 CPU s in a Celeron PC, with a 466 MHz processor. Using strategy B, 5 outer iterations are required for stopping the search, consuming 565.18 CPU s for solving the bounding problems. The last problem GR(x,y,w)K has 138 binary variables and 604 continuous variables. Last, the model is modified imposing the no stream splitting condition. Hence, the previous solution is no longer feasible, since it supposes the split of stream H2 in stages 1 and 2. If no stream split is allowed, the optimal solution has a total cost of $243 089.31, involving the exchange between stream H1 with C2 and C1 (in a decreasing order in temperature) and the stream H2 with C1, C4, and C3, in the same order. Moreover, streams H1, H2, and C2 use cooling and heating utilities. This modified problem has less feasible structures. The algorithm requires 5 iterations to converge, using only 20.86 CPU s in solving the bounding problems GR(x,y,w)K. 7.3. Example 3. This example deals with a system of 4 hot streams and 5 cold streams, widely studied in the literature, known as 9SP.16 As far as we know, this problem, and none other of comparable size, has been treated with global optimiza-

tion methodologies. The fixed cost of units is $2000, the area cost coefficient is $70, and the scale factor is β ) 1. Pettersson18 reports an optimal solution of 2905 M$/year. This solution corresponds to an approximated MILP model where the temperature is divided into subintervals. The superstructure considered by Pettersson has more degrees of freedom than the one considered here, because of the possibility of stream splitting, mixing of streams at different temperatures, and placement of coolers and heaters in intermediate stages. The reported network has a complicated topology, involving 17 units, 7 stream splitting, and intermediate heaters. Lewin19 also solved this example, finding a network with cost 2936 M$/year, allowing stream splitting, and a simpler network without stream splitting, with a cost of 2946 M$/year. The superstructure with 5 stages generates a model with 109 binary variables, 611 continuous variables, and 1126 constraints. The no stream splitting condition is imposed. Then, the integer cuts proposed in Section 5 are included in the bounding problem, as well as the Amin vs Q curve. To have a rough measure of how this curve improves the underestimation, note that the optimal solution of the first GR(x,y,w)K problem (using two gridpoints for each variable) without imposing the curve is 1 814.60 M$/year, whereas the optimal solution of the problem with this curve is 2 338.92 M$/year. The proposed algorithm finds the global optimal solution of the MINLP model (within a tolerance of 1%) in the outer iteration number 18, and it stops in iteration 22, with the infeasibility of the bounding problem. The optimal solution has a cost of 2 935.02 M$/year. The optimal structure consists of 10 exchangers, 3 coolers, and 2 heaters, and it is depicted in Figure 8. Note that the solution has a cost near the cost obtained with the approximated model by Pettersson, but with a much simpler structure. Figure 9 shows the T vs ∆H diagram for the exchangers between process streams in the optimal solution of example 3. The algorithm consumes ∼104.86 min, distributed as follows: 84.68% is spent in solving the bounding problems GR(x,y,w)K, 8.9% is spent in bound contraction, and 6.5% is spent in solving the bounding problems LR(x,yl,w)K. The time consumed in solving the nonconvex NLP problems P(x,yl) is negligible. Strategy B was used for updating the grid, when the relative approximation error in each nonconvex term is >0.1. The initial


Figure 12. Progress of the best-known solution in example 4.

grid considers two points (lower and upper bound of areas). The set of linearization points for the first iteration in the OAF procedure are the four corners and the center point of the domain up lo up [∆Tlo 1 , ∆T1 ] × [∆T2 , ∆T2 ]. The lower bound for ∆Ti is EMAT ) 10 °C. CPLEX runs with the default options, and the search for a feasible solution is finished when a feasible solution is within an absolute tolerance of 800 (generally, the first feasible solution satisfies that tolerance). The optional choices of the algorithm enable it to converge, but when another combination of optional values was used, the algorithm could not converge in

Global Optimal Structures of Heat Exchanger Networks by Piecewise

Recommend Documents