Improved sequential strategy for the synthesis of near-optimal heat exchanger networks

Computers chem. Engng, Vol. 21, Suppl., pp. S59-S64, 1997

Pergamon

© 1997 Elsevier Science Ltd All rights reserved Printed in Great Britain

PII:S0098-1354(97)00026-4

0098-1354/97 $17.00+0.00

Improved Sequential Strategy for the Synthesis of Near-Optimal Heat Exchanger Networks by Truls Gundersen t~, Petter Traedal 2 and Abdolreza Hashemi-Ahmady 2 Norwegian University of Science and Technology, Trondheim, Norway Telemark Institute of Technology, Porsgrunn, Norway 2 (~:): Corresp. to this author: Dept. of Thermal Energy and Hydro Power, NTNU, N-7034 Trondheim, Norway Phone: + 47-73-59-2700, Fax: + 47-73-59-8390, E-mail: [email protected]

Abstract A sequentialfraniework (LP, MILP, NLP) for Heat Exchanger Network Synthesis (HENS) is described that reduces significantly the numerical problems caused by non-linearities, discontinuities and combinatorial explosion in simultaneous MINLP models. The multiple trade-offs in HENS related to level of heat recovery, network topology and area distribution as well as the trade-off between total annual cost, network complexity and operability is dealt with by iteration, user interaction and a new MILP transportation model. This model, which is based on the same idea as our vertical MILP transshipment model, is crucial in the sequential approach, since it provides sets of matches that give designs with close to minimum area and cost when supplied to an NLP model for network generation and optimization.

In addition to the framework itself, this paper will focus on specific issues related to the sequential approach such as: ( 1) improving the selection of matches in the MILP model by accounting for temperature driving forces, (2) reducing the combinatorial problem of the MILP model by a tighter formulation, and (3) improving the robustness and efficiency of the NLP model by increased use of information from the MILP model and the addition of convex estimators. Results will be presented for some of these problem areas.

Introduction

Previous Work

Heat Exchanger Networks of industrial size can currently be synthesized using Pinch Technology. The tedious nature of the Pinch Design Method and its inability to properly address the multiple trade-offs in such design activities, however, ask for the use of simultaneous optimization techniques, either as a stand-alone tool or used in conjunction with Pinch Analysis based methods.

The first breakthrough in using Mathematical Programming to solve HENS problems was the use of transportation models (Cerda et al., 1983, and Cerda and Westerberg, 1983) and transshipment models (Papoulias and Grossmann, 1983) for minimum utility requirements, fewest number of units and the corresponding Heat Load Distributions (HLDs), even with forbidden matches. An HLD is a set of the matches (i,j) and their duties (Qi.j). The network generation and optimization was realized by a stream superstructure (Floudas et al., 1986). The result was a sequential approach (LP, MILP, NLP) to HENS. An important part of these developments was the use of insights from Pinch Analysis, such as temperature intervals and pinch decomposition.

Simultaneous MINLP models have been published that in theory can address and solve the trade-offs in Heat Exchanger Networks. Due to severe numerical problems related to the non-linear (non-convex) and discrete (combinatorial) nature of these models, however, local optima are often encountered and the size of problems that can be solved does not meet industrial needs.

Later, models and numerical solution strategies for simultaneous MINLP models were developed (Yee and Grossmann, 1990, and Ciric and Floudas, 1991). The motivation behind these models was limitations in the sequential approach related to improper account for the trade-offs involved, and most important, the fact that HLDs were selected without considering heat transfer area in the corresponding networks.

The iterative and sequential framework proposed in this paper provides a compromise between Pinch Technology and MINLP models with the objective to solve problems of industrial size and find near-optimal networks. Related work is carded out by Zhu et al., 1995. An important feature of our approach is a decomposition into more manageable subtasks. Level of heat recovery and number of matches are adjusted iteratively by the user after proper and automated initialization. For each pair of these two key variables, the selection of matches and generation as well as optimization of the corresponding network is obtained by the use of a new MILP model followed by an NLP optimization.

There are, however, severe numerical problems related to these MINLP models caused by non-convexities (local optima) and combinatorial explosion (computer resources). At present, it is virtually impossible to solve problems with more than 20 streams. In fact, even S59

$60

PSE '97-ESCAPE-7 Joint Conference

problems with 10-15 streams are sometimes hard to solve without resorting to supercomputers. In current research, the combinatorial problem is addressed by the use of hard and soft logic (Raman and Grossmann, 1991), while the problem with local optima is addressed by the use of convex under-estimators (Quesada and Grossmann, 1993 and 1995, and Ryoo and Sahinidis, 1995).

The Sequential Framework The motivation behind the sequential framework described in this paper is the fact that little progress has been made during the last years to improve robustness and size of the problems that can be solved by the simultaneous MINLP models, except for certain classes of rather simplified problem definitions. The objective of the research on the sequential approach is to be able to solve problems with at least 30 streams, which would make it possible to address most industrial cases. The sequential framework is based on the recognition that the selection of matches (HLDs) is of crucial importance, not only for the cost of the network, but also for qualitative aspects of complexity, operability and controllability. The vertical MILP model for selection of matches and the subsequent NLP model for generating and optimizing the corresponding network are the inner elements of the sequential approach (see figure 1). The level of heat recovery (represented by HRAT, the Heat Recovery Approach Temperature,) and the number of units (U) are adjusted iteratively by the user after proper and automated initialization. HRAT is initialized by the pre-optimization procedure (SuperTargeting - ST) by Ahmad et al., 1990, and Linnhoff and Ahmad, 1990, while U is initialized for the corresponding value of HRAT by a simple MILP transshipment model allowing the Exchanger Minimum Approach Temperature (EMAT) to take values less than HRAT to allow more feasible solutions (for example EMAT = 0).

[-~ AdjustHRAT

Matches and Topologies The topological aspects of Heat Exchanger Networks can be put into layers of a hierarchy of decision making. Related to these issues is the concept of topology traps discussed by Sagli et al., 1990, and Trivedi et al., 1990. Selecting the matches (i,j) and corresponding duties (Qi.j), is the top level decision. For each such HLD, however, there are sub-decisions that have to be made regarding the sequence of matches, the choice between series or parallel heat exchange, and the possible use of bypasses. For HLDs (and corresponding networks) with more than minimum number of matches (units), there are also degrees of freedom with respect to the duties of the matches. Thus, there are no longer well defined duties (Qij) corresponding to a set of matches (i,j). This is pointed out to explain that the NLP optimization following the MILP model in the sequential framework is not an easy task and becomes harder the higher the number of units is above the minimum. Nevertheless, the NLP model is much simpler to solve than the MINLP model in simultaneous approaches. One reason is that it is much easier to develop tight convex estimators to counteract the problems with local optima caused by nonconvexities in the model. The model for selecting HLDs in figure ! is termed vertical MILP due to our previous research trying to minimize total area by maximizing vertical heat transfer between the composite curves. As will be described later, work to improve this model has lead to a switch from transshipment to transportation models, however, the idea of vertical heat transfer is still applied, although in a less explicit manner. Thus, we decided to continue referring to the model as a vertical MILP. The Vertical MILP

Final Network IfllAT

IAdjustUnits Fig. 1

the minimum number of units) and 3) adjust the value of HRAT slightly above and below the initial value.

Sequential framework for HENS.

The designer is now in a position to evaluate the first generated network with respect to Total Annual Cost, network complexity (number of units, stream splits, etc.), operability and controllability, before exploring other alternatives. The logical sequence of actions is indicated in figure 1 as the following nested loops: 1) derive networks for the second or third "best" HLDs, keeping the number of units and HRAT unchanged, 2) increase the number of units by one (the system always starts at

Model

The original vertical MILP model is described by Gundersen and Grossmann, 1990. Briefly, it is based on the idea that vertical heat transfer between the Composite Curves will improve the use of available driving forces and thus reduce total area. The MILP transshipment model for fewest number of units involves a Heat Cascade based on Temperature Intervals (TIs). When studying vertical heat transfer, the Composite Curves are also divided into Enthalpy Intervals (EIs) at every kink point of the curves. The model allows matches according to the TIs based on EMAT, while the quality of each match is determined from the EIs based on HRAT. A penalty term is added to the objective function that compares the actual heat exchange between two streams with the maximum amount that can be transferred vertically between the same streams. Since a match in the general case will involve heat exchange in more than one TI, the actual non-vertical heat transfer may be larger

PSE '97-ESCAPE-7 Joint Conference in the final network where each multi-interval match is transformed into a single exchanger. The penalty term (Si.j) is therefore only a lower bound on criss-crossing, but still a good indication of the ability of the final network to transfer heat vertically and save total area.

Example 1 Table 1 shows stream data for Example 1 which is used to test previous and new models for selecting HLDs. In case (a), two of the streams (HI and C1) have significantly lower film heat transfer coefficients than the other streams. In case (b), all process streams have equal film heat transfer coefficients. Table 2 shows the sequence of HLDs and the corresponding networks obtained by the original vertical MILP model for Example l a. The various designs have been named according to their capital cost (A is cheapest, G is most expensive). Heat exchanger cost is obtained from:

Note that only capital cost is considered in this paper, since the various HLDs and their corresponding networks are derived for the same value of HRAT. Operating cost is the same for all designs, and the optimization of HRAT (the outer loop in figure 1) is not the issue here. Even though the original vertical model worked nicely for problems with similar heat transfer coefficients, the results in table 2 are no____!satisfactory. t It is interesting to note that three designs (A, B and C) in table 2 have total area less than the approximate target, 726 m2, calculated from the Bath formulae (Townsend and Linnhoff, 1984). Stream Data for Example 1a and b

Str.

T&Y3

Td_°_) mCp(kW/°C)

H1 H2 H3 C1 C2 C3 ST CW

300 200 190 160 180 190 250 30

200 190 170 180 190 230 250 50

Table 2

10 100 50 50 100 25

h(kW/m 2°C) 0.1 1.0 1.0 0.1 1.0 1.0 4.0 2.0

1.0 1.0 1.0 1.0 1.0 1.0 4.0 4.0

Original Vertical MILP, Ex. 1a, U=4

Seq.

Des

1 2 3 4 5 6 7

B E D A C F G

Y~ZSij(kW) 0 2000 2000 2000 3000 3000 4000

area. Gundersen and Grossmann, 1990, explained the tw__._Qeffects that should be considered in such cases.

Shifting means that a hot stream with high film heat transfer coefficient should be used to heat up a cold stream at higher temperatures than strict vertical heat transfer indicates. Shifting is taken care of in the vertical model by individual stream contributions (ATi) to the approach temperature and use of modified temperatures. The second effect, pairing, is accounted for in an extended vertical model (Gundersen et al., 1996) by adding a new penalty term to the objective function, where designs involving matches between streams with similar heat transfer conditions have high score (low penalty term). This principle was put forward by Umeda et al., 1978. The objective function of the extended vertical model is: mint a. YEyij+ 6" EES~j+ ~'. EEy~j.[ 1-min(hflhj, hJh~)]

Cost = 8600 + 670. (A) 0.83 (m2 and US$)

Table 1

$61

Area(m2)

Cost(kS)

699.5 1352.5 832.8 519.7 696.3 1385.1 1485.9

214.1 331.9 241.0 177.1 219.4 343.0 358.7

For cases with significant differences in heat transfer coefficients (Example l a), strict vertical heat transfer does not always result in minimum total heat transfer

a is set to a small value (the model is used with a fixed number of units, but keeping this term in the objective function improves the search). ~ and y are set to values that make the second and third term (0-1) normalized. Table 3 shows results for the improved vertical MILP model. The sequence of HLDs found by this model is almost perfect. Only design F is identified "out of order". The explanation is that, for this network, two exchangers have to operate at AT=10°C throughout the unit, which is considerably less than HRAT=20°C.

Table 3

Improved vertical MILP, Ex. 1a, U=4

Seq.

Des

Y.YSij

f(hi,h~)

1 2 3 4 5 6 7

A B C F D E G

1092 1030 2366 2456 2440 2500 3910

1.250 3.050 2.375 2.600 3.050 3.050 3.050

Obj.f.

Cost(kS)

0.59 1.02 1.19 1.26 1.37 1.39 1.74

177.1 214.1 219.4 343.0 241.0 331.9 358.7

A New Transportation Model Table 3 and the incorrect position of design F in the ranking of HLDs shows the obvious fact that non-vertical heat transfer at low driving forces is worse than such heat transfer at large driving forces. One potential solution would be to add weights to each non-vertical term in the objective function, such as the inverse of the temperature difference for the actual streams. In order to find the driving forces for a certain match, however, it is necessary to keep track of where (in which interval) heat is delivered from a hot stream and where (in which interval) heat is consumed by a cold stream. Further, with modified temperatures to account for the effect of shifting, it is impossible to have corresponding Enthalpy and Temperature Intervals. In conclusion, the

$62

PSE '97-ESCAPE-7 Joint Conference

extended vertical MILP model is unable to account for where non-vertical heat transfer occurs.

additional EIs are introduced by adding (subtracting) EMAT to (from) cold (hot) stream temperatures.

The need to keep track of where (at what temperatures) heat is introduced and removed from the heat cascade in addition to the corresponding hot and cold streams, means that the number of subscripts for transferred heat becomes four (Qim,jn , where i,j are hot and cold streams and m,n are intervals). The advantage of transshipment models versus transportation models is thus lost.

Figure 2 shows heat flows in the transportation model. A hot stream (i) is divided into smaller amounts of heat that flow into the various EIs (m-I, m and m+l). Likewise, a cold stream (j) is divided into smaller amounts of heat required in the various EIs (n-I, n and n+l). The EIs are not used as intermediate storage (transshipment models), only to make the division of streams possible. Whether heat can flow from a hot sub-stream (i,m) to a cold substream (j,n-1) in figure 2, depends on the temperatures.

Another disadvantage with the extended vertical MILP transshipment model is the need to assign heuristic values to the weights 13 and ~/. There are also heuristic elements in the use of stream individual AT contributions as well as in the third term accounting for pairing. This paper describes a new transportation M I L P model for selecting HLDs, where the objective function is:

min Zi Y,j Zm Zn [ Qim,jn/ ( Uij. ATLM,mn) ] Uij is the heat transfer coefficient for a match between hot stream i and cold stream j, and ATLM.m,is the log mean temperature difference for heat transfer between intervals m and n. It is important to note that Uij and ATLM,mnare constants that can be calculated ahead of time, thus the model is still linear. It is, of course, possible to extend Uij to account for variations with temperature by including subscripts for intervals m and n (Uim,j n ). Also note that the objective function is now a very realistic estimate of area. To go one step further would involve fixed charge terms (involving binary variables) and an economy of scale type exponent to convert area into capital cost. The fixed charge terms could easily be implemented, however, in the sequential framework, the MILP model for selecting matches (HLDs) is applied with a fixed number of units. Using economy of scale type cost law rather than area in the objective function would transform the MILP model into an MINLP model, which is what we wanted to avoid in the first place. Shifting and pairing are no longer treated as two distinct phenomena, but taken care of as a combined effect (as it should be) of Uij and ATLM.mnin the objective function.

Fig. 2

New Transportation Model for HENS

In line with the basic idea of vertical heat transfer, the intervals used are Enthalpy Intervals (EIs) derived from the composite curves drawn according to HRAT. To allow for non-vertical heat transfer, where appropriate,

Choosing a value for EMAT is not as easy in this model as it was in the vertical transshipment model, where EMAT = 0 was used in the TIs. In the transportation model, the value of EMAT used to create additional EIs has to be balanced. If it is set too low, the corresponding values of ATLM,mnbecome very small, and such nonvertical heat transfer will face large penalties. If, on the other hand, EMAT is set too high, potential HLDs will be excluded from the feasible set of solutions. In this paper, we have used EMAT = HRAT / 4. Table 4 shows the results when the new transportation model is applied to Example la. As can be seen, the sequence of solutions produced by the model (using integer cuts to find subsequent HLDs) is perfect. In table 5, the same model is applied to Example lb (where all process streams have identical heat transfer coefficients). Table 4

Transportation Model, Ex. 1a, U=4

Seq.

De..._~s

Obj.func

Areafm2)

Cost(kS)

1 2 3 4 5 6 7

A B C D E F G

547.7 698.9 716.5 832.4 1352.6 1405.0 1486.0

519.7 699.5 696.3 832.8 1352.5 1385.1 1485.9

177.1 214.1 219.4 241.0 331.9 343.0 358.7

Transportation Model, Ex. lb, U=4

Table 5 Seq.

De..._ss

Obj.func

Area(m2)

Cost(kS)

1 2

B D

163.9 175.8

164.0 175.9

91.2 94.1

3 4

E A

282.0 293.0

281.9 272.5

119.5 116.7

5 6

G C

293.9 311.0

293.8 291.3

122.4 121.6

7

F

461.4

441.6

159.4

We have kept the names of the various networks (A,...,G) from case (a). It is interesting to note how dependent the relative costs of the networks are upon the differences in heat transfer coefficients. Also note that the area target from the Bath formulae in this case is 159 m2, while the best designs have areas of 164.0 and 175.9 m 2.

$63

PSE '97-ESCAPE-7 Joint Conference In order to see how the new model performs, it is necessary to look in the column for cost. Clearly, the model is still able to distinguish good HLDs from bad ones. The fact that the various designs are produced in a sequence which is slightly "imperfect", should not be paid too much attention. Designs that have switched places according to the perfect ranking, are very similar in both area and cost. We have grouped the designs in table 5 according to their costs to highlight this point.

Example 2

-

U~j.y~ ~ 0

can be reduced by assigning the smallest possible value to the upper bounds Uij. _Onefeasible value that has been used in the past is: Uij = min [Y-,mQim , 53nQjn] This value, however, only considers the amount of heat, not the temperature where heat is available or required. Much lower values can be obtained by considering mCpvalues and temperatures for each pair of streams:

mCp

15oo 8o. D

[20]

°

[so]

90 °

[ 10001075 kW

I

1550 kW

<

400kW

1

25°

@ 1800 kW

,30,

450 kW

Fig. 3 Optimal Network for Ex. 2 (HRAT=20°C, U=5). Figure 3 shows the best design (A) for Example 2, when HRAT is specified at 20°C, and the number of units is 5. Table 6

Stream Data for Example 2

Str.

~

!,(°)

mCp(kW/°C)

h(kW/m:°C)

H1 H2 C1 C2 ST CW

150 90 20 25 180 I0

60 60 125 100 180 15

20 80 25 30

0.1 0.1 0.1 0.1 0.1 0.1

Table 7

1 2 3 4

We have looked at two ways to make the formulation of the MILP model tighter. First, the gap introduced in the logical constraints that relate binary variables y~j and continuous variables Qim,jn: XmZnQim,jn

Table 6 shows stream data for Example 2, a problem that has been thoroughly discussed in the literature. Again, we have only looked at solutions with minimum number of units (in this case Umin= 5 for HRAT = 20°C, which is used here). Table 7 shows results for Example 2 using the new transportation MILP model. The sequence of HLDs found are again perfect, and demonstrates the ability of the model to select HLDs that result in Heat Exchanger Networks with low total area and capital cost.

Seq.

factor is the difference (gap) between the values of the objective functions for the actual MILP model and its corresponding relaxed LP model. An MILP model with zero gap can be solved as an LP problem.

Transportation Model, Ex.2, U=5 De..___ss O~.func A B C D

3245.9 3290.1 4178.4 4230.9

Area(m:)

Cost(kS)

3146.6 3203.6 3613.4 4054.9

717.3 730.8 801.1 853.2

Uij = min [ -~'mQtm , ~n Qjn , max t min (mCpi, mCpj) • (Ts, i - Tsj- EMAT), 0 ~ ] The third term (max/min) will detect cases where the hot supply temperature (Ts,i) is less than the cold target temperature and cases where the cold supply temperature (Ts,j) is higher than the hot target temperature (adjusted for EMAT). Since this term sometimes gives negative duties, we have to include zero and a max operator. The first two terms are included since the third term gives larger duties than available or required in cases where the hot target temperature is higher than the cold supply temperature and cases where the cold target temperature is less than the hot supply temperature (adjusted for EMAT). Of course, the smallest mCp value of the streams considered determines the critical values. By implementing these Uij values for a larger problem with 13 streams and 17 units (not presented here), the gap was reduced by only 10%, which is disappointing. Other problems show similar reductions, which is not enough to make it possible to solve 30 streams problems. Keeping the new values for Uij and adding similar local constraints (using Uim,jn) on the heat exchanged:

Qimjn - Uim,jn'Yij

--< 0

resulted in gap reductions on the order of 20%, however, some of the corresponding computer time savings are canceled by the fact that the model is growing in size and more iterations are required to solve the LPs at each node of the branch and bound search tree.

Tighter MILP Models The difficulty of solving MILP problems is not only related to the number of binary variables. A more critical

Another idea that we have briefly tested, so far with little success, is to put upper bounds on the number of units for each stream, based on the total number of units and the

$64

PSE '97-ESCAPE-7 Joint Conference

duty of the actual stream relative to the total amount of heat exchanged.

Conclusion A sequential framework with the potential to identify near-optimal Heat Exchanger Networks is presented. A crucial part of the approach is a new MILP transportation model for selecting Heat Load Distributions for a given number of units and level of heat recovery. The model is based on the concept of vertical heat transfer that was introduced in an earlier MILP transshipment model. Preliminary results show excellent ability of this model to select HLDs which lead to networks with low total area requirements and low total cost. Attempts have also been made to reduce the gap of these MILP models and thus increase the size of the problems that can be solved in reasonable computer times, so far with limited success.

Future Work More testing of the new transportation model is required, however, the results presented in this paper are indeed promising. To fulfill the ambitions of the sequential framework to be able to handle industrial problems with up to 30 streams, considerable improvements in the formulation and solution of the MILP model have to be achieved. Crucial elements are the gap of the MILP model and the branch and bound search strategy. Additional use of thermodynamics, logic and heuristics is required. Further, the robustness of the NLP model to generate and optimize the actual networks should be improved by the addition of tight convex under-estimators and increased use of information available from the MILP model in setting up and initializing the NLP model. An example of the latter is the use of information about where (in which interval) heat is exchanged between the various streams when setting up and pruning the superstructure.

References

Floudas C.A., Ciric A.R. and Grossmann I.E., 1986, Automatic Synthesis of Optimum Heat Exchanger Network Configurations. AIChE Jl., 32, pp 276-290. Gundersen T., Duvold S. and Hashemi-Ahmady A., 1996, An Extended Vertical MILP Model for Heat Exchanger Network Synthesis. Comput. Chem. Engng., 20, (Proc. from ESCAPE6, Rhodes), pp $97-S 102. Gundersen T. and Grossmann I.E., 1990, Improved Optimization Strategies for Automated Heat Exchanger Network Synthesis through Physical Insights. Comput. Chem. Engng., 14, pp 925-944. Linnhoff B. and Ahmad S., 1990, Cost Optimum Heat Exchanger Networks - 1. Minimum Energy and Capital using Simple Models for Capital Cost. Comput. Chem. Engng., 14, pp 729-750. Quesada I. and Grossmann I.E., 1993, Global Optimization Algorithm for Heat Exchanger Networks. IEC. Res., 32, pp 487-499. Quesada I. and Grossmann I.E., 1995, A Global Optimization Algorithm for Linear Fractional and Bilinear Programs. Jl. Global Optim., 6, pp 39-76. Papoulias S.A. and Grossmann I.E., 1983, A Structural Optimization Approach to Process Synthesis - II. Heat Recovery Networks. Comput. Chem. Engng., 7, pp 707-721. Raman R. and Grossmann I.E., 1991, Logical Inference in Branch and Bound Search for Process Synthesis. Proc. from PSE'91, I, Montebello, Canada, pp 10.I-10.15. Ryoo H.S. and Sahinidis N.V., 1995, Global Optimization of non-convex NLPs and MINLPs with Applications in Process Design. Comput. Chem. Engng., 19, pp 551-566. Sagli B., Gundersen T. and Yee T.F., 1990, Topology Traps in Evolutionary Strategies for Heat Exchanger Networks. In Comp. Appl. in Chem. Engng., ed. by Bussemaker H.T. and Iedema P.D. (Proc. from CACE'90, The Hague, May 1990), Elsevier, Amsterdam, pp 51-58. Townsend D.W. and Linnhoff B., 1984, Surface Area Targets for Heat Exchanger Networks. IChemEAnnl. Res. Mtg., Bath, UK. Trivedi K.K., O'Neill B.K., Roach J.R. and Wood R.M., 1990, Systematic Energy Relaxation in MER Heat Exchanger Networks. Comput. Chem. Engng., 14, pp 601-611. Umeda T., Itoh J. and Shiroko K., 1978, Heat Exchange System Synthesis. Chem. Eng. Progr., 74, pp 70-76.

Ahmad S., Linnhoff B. and Smith R., 1990, Cost Optimum Heat Exchanger Networks - 2. Target and Design for Detailed Capital Cost Models. Comput. Chem. Engng., 14, pp 751-767.

Yee T.F. and Grossmann I.E., 1990, Simultaneous Optimization Models for Heat Integration - II. Heat Exchanger Network Synthesis. Comput. Chem. Engng., 14, pp 1165-1184.

Cerda J., Westerberg A.W., Mason D. and Linnhoff B., 1983, Minimum Utility Usage in Heat Exchanger Network Synthesis - A Transportation Problem. Chem. Eng. Sci., 38, pp 373-387.

Zhu X.X., O'Neill B.K., Roach J.R. and Wood R.M., 1995, A New Method for Heat Exchanger Network Synthesis using Area Targeting Procedures. Comput. Chem. Engng., 19, pp 197-222.

Cerda J. and Westerberg A.W., 1983, Synthesizing Heat Exchanger Networks having Restricted Stream/Stream Match using Transportation Problem Formulations. Chem. Eng. Sci., 38, pp 1723-1740. Ciric A.R. and Floudas C.A., 1991, Heat Exchanger Network Synthesis without Decomposition. Comput. Chem. Engng., 15, pp 385-396.