Synthesis of heat exchanger networks using genetic algorithms

Synthesis of heat exchanger networks using genetic algorithms

Available online at Applied Thermal Engineering 28 (2008) 1763–1773 Synthesis of heat excha...

248KB Sizes 4 Downloads 50 Views

Available online at

Applied Thermal Engineering 28 (2008) 1763–1773

Synthesis of heat exchanger networks using genetic algorithms Jean Dipama a, Alberto Teyssedou a,*, Mikhaı¨l Sorin b a

Nuclear Engineering Institute, Engineering Physics Department, E´cole Polytechnique de Montre´al, P.O. Box 6079, succ. Centre-ville, Montreal, Que´bec, Canada H3C 3A7 b CANMET Energy Technology Centre, 1615, boul. Lionel-Boulet, Varennes, Que´bec, Canada J3X 1S6 Received 18 July 2006; accepted 12 November 2007 Available online 23 November 2007

Abstract This paper presents a methodology for carrying out both the synthesis and the optimization of heat exchanger networks (HEN) based on a genetic algorithm (GA) technique. The proposed methodology allows the topology of the HEN and the heat load distribution that satisfy an energy target to be obtained from a single pass calculation. The technique is characterized by a high reliability in determining the best HEN structure. In addition, it also allows the generation of several HEN configurations to be obtained; thus, a best HEN selection in accordance to a given application can be easily achieved. The proposed approach permits the HEN topology to be obtained only if flow stream divisions are avoided. The simplicity of the proposed procedure makes this technique useful when a fast solution of a problem is looked for. The validation of the algorithm is carried out by solving benchmark tests taken from the open literature. In all the cases studied, the results compare quite well with those given in the open literature. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Heat exchanger networks; Optimization; Energy recovery; Genetic algorithms

1. Introduction The synthesis of a HEN is usually a very complex task that implies a combinatory problem for matching hot and cold flow streams in order to permit a maximum energy recovery to be achieved. In general, the HEN synthesis process can be summarized as follows: A set of hot flow streams must be cooled to specific temperature values, while another set of cold flow streams must be heated up to determined values. Each flow stream is characterized by its own specific heat capacity and mass flow velocity. Thus, the problem to be solved consists in finding the optimal topology (i.e., heat exchanger structure) of a HEN having the most appropriate heat load distribution in such a way that the maximum thermal power can be transferred between the flow streams. It is obvious that the optimization process also must reduce the number of external utilities (i.e., heat sources and sinks). Therefore, *

Corresponding author. Tel.: +1 514 3404711; fax: +1 514 3404192. E-mail address: [email protected] (A. Teyssedou).

1359-4311/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2007.11.014

the synthesis of the optimal HEN requires working in two different ‘‘solution spaces”: (a) a topological space where, according to the nature of the interaction between the flow streams different structures are possible and (a) a thermal load space where different thermal power distributions between heat exchangers are possible. Thus, the optimization process consists in simultaneously finding the HEN topology and the heat load distribution that permits a maximum energy recovery to be satisfied without violating thermodynamics principles. In most of the works carried out on this subject, the topology and the thermal power distributions are treated separately. The technique proposed in this paper allows the aforementioned two processes to be simultaneously handled by using a genetic algorithm (GA). The research work in this area have been essentially developed according to three different concepts: the pinch method that is based on thermodynamic principles [1];


J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

the use of mathematical methods [2–4] and the methods based on the use of artificial intelligence first used among others by Grimes et al. [5]. One of the most important concepts used in the first method cited above consists of constructing heat load composite curves. These curves provide detailed information on the energy recovery, the minimum heating and cooling power required by the whole process as a function of a minimum temperature admissible across the heat exchangers. Consequently, these curves permit the overall thermal power in conjunction with the heat transfer surface required by the process, to be determined. The pinch technology, successfully used in a great number of industrial applications, requires two steps [1]; first it targets the maximum energy recovery, and then it synthesizes the heat exchangers network. This method supposes that the flow streams, the supply temperatures and the flow heat capacities are known. Even though, this technique is based on thermodynamic principles, it also uses several heuristic rules. These heuristic rules constitute the major disadvantage of this method, because they are based on approximations which are sometimes difficult to validate. Moreover, in some cases, the application of the pinch method may be quite lengthy and cumbersome. Even though, the rules used for matching the hot and cool flow streams are simple, finding an initial network design is usually a very difficult task. The decomposition in two pinch sub-systems, limits the evolution of the design because it becomes conditioned by the initial halved design [6]. These limitations have motivated the search of more efficient mathematical techniques for handling HEN optimization problems. In general, a mathematical modeling technique consists in an ensemble of equality and inequality equations that must satisfy the minimization of an objective function. The resulting system of equations can be written under the following general form: Objective function to be minimized : f ðX ; Y Þ; Constrained equations : gðX ; Y Þ P 0; hðX ; Y Þ ¼ 0; ð1Þ where X 2 Rn and Y 2 [0, 1]m. Note that X and Y are vectors of variables while g and h are vectors of functions. For the case where there are no binary variables (i.e., dim(Y) = 0), and the f, g and h are linear functions, the mathematical problem is reduced to the simplest case, that is a ‘‘Linear Programming” (LP) system. If the model includes binary variables (i.e., dim(Y) 6¼ 0) and f, g, and h are linear functions, the problem is of the ‘‘Mixed Integer Linear Programming” (MILP) type. For the synthesis and optimization (maximum energy recovery) of a HEN without flow stream division then, the optimization problem is of the MILP type. However, the applicability of this technique requires two questions to be addressed simultaneously. The first one concerns the solution of an integer problem that results from the discretization of the temperatures and enthalpies. It is important to remark that moderate size networks may

require a large number of intervals; thus, resulting in very large linear system of equations to be solved. The second question concerns the parametric optimization itself that usually is handled by another mathematical solver, such as the Simplex technique. Duran and Grossmann [7] have proposed a method that optimizes simultaneously the process as well as the heat exchanger network. In this method, the flow streams, as well as the supply and target temperatures of the flow are not imposed, but they result from the optimization process itself. This method also permits the minimum approach temperature to be determined. A MILP model is then used to synthesize the optimal heat exchanger network for the maximal heat recovery, by using the minimum number of utility units. In order to reduce the size of the system of equations, different approaches have been proposed such as: the use of the stage concept [8] and the use of the block concept [9]. Even thought, the use of the aforementioned techniques can help in reducing the size of the problem to be handled, some industrial applications can still produce quite large optimization schemes that cannot be easily solved. Recently, new stochastic methods based on GA’s have been successfully used in several engineering applications. The strategy of these algorithms is founded on the principle of natural evolution as postulated for the biological spices, where a large number of possible solutions can be handled. This is in contrast to classic mathematical methods that converge through a solution by evaluating the local gradient of the functions. According to the authors’ knowledge, only few works have been published where the GA have been used for synthesizing a process including HEN, starting from pre-established known topologies [10]. In some cases the GA were used to find the optimal heat exchange topologies while other methods were used for obtaining the parametric optimization [11]. Similarly, Lewin et al. [12] have proposed a procedure based on the use of GA’s for the optimization of the heat exchanger topology from a compact and constant representation of the HEN. This representation is used for both, satisfying the requirements of the GA itself and for carrying out the parametric optimization. Recently, Ravagnani et al. [13] have used GAs for synthesizing HENs based on a previous optimization of the network pinch temperature that permitted the authors to determine the optimal power of the utilities. A design methodology of HEN including flow stream splitting based on the use of a randomized algorithm was developed by Pariyani et al. [14]). In these two works, the authors have used the decomposition technique across the HEN pinch point for carrying out the synthesis of HENs. Thereafter, each sub-network is optimized separately by using evolutionary algorithms. It is important to remark, however, that in both cases the best HEN solution cannot be achieved without some external intervention of the user. The method presented in this paper, however, assumes that the flow streams are completely defined including supply and target temperatures, minimum approach temperatures and flow heat capacities. We use the table algorithm

J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

(TA) [1] to target the Maximum Energy Recovery (MER), while the genetic algorithm allows both the topology and the parametric optimization of a HEN to be simultaneously determined. Even though, the proposed method does not support stream splitting, it can be applied independently of the complexity of the problem because it requires neither the network decomposition across the pinch nor the intervention of the user. The methodology is very robust, easy to implement and it can be used for synthesizing complex networks by using a personal computer.

2. General overview of the genetic algorithm In general, a GA represents a methodology that permits a solution to be searched in a ‘‘solution space” by following a procedure similar to that proposed for the natural selection and evolution of biological species [15,16]. Considered as an optimization technique a GA examines and manipulates simultaneously a large population of possible solutions. Each solution is identified using a symbolic mathematical representation. The search procedure starts from a set of initial possible solutions that are represented by ‘‘chromosomes”; this set is called the ‘‘initial population”. The solutions presented in a population are then chosen according to a fitness criterion and are then used in order to generate a new population of solutions. This multiplication and selection procedures permit the ‘‘quality” of the new population to be improved with respect to the initial one. The generation and the selection procedures of new populations are repeated until a given convergence condition is satisfied. The general form of a GA can be summarized as follows: 1. Start with a random generation of an initial population of N chromosomes. 2. Carry out a fitness evaluation, f(x), for each x-chromosome forming the population. 3. Apply a cross-over operation to the population in order to generate a new one according to the following steps: (a) Select two parent chromosomes according to their best fitness. (b) Use a cross-over probability in order to reproduce the two parents into two new chromosomes (offspring’s). Note that if crossing the parents is not carried out, then the offspring’s become an exact replica of their parents. (c) Use a mutation probability to modify the new chromosomes. (d) Relocate the new chromosomes in the population space. 4. Use this new population for continuing searching the best solution, i.e., continue the execution of the algorithm. 5. Carry out a test for satisfying a convenient convergent criterion, if this condition is achieved stop the procedure


and select the chromosome that has the best fitness as the solution of the problem. 6. If step 5 is not satisfied then go back to step 2.

3. The proposed approach In order to search the most appropriate solution within a ‘‘solution space”, it is obvious that a consistent representation of this space is required. It is important to remark that the transformation of a real world problem into a solution space necessitates the use of an appropriate symbolic representation. Thus, the search of the most ‘‘adapted” solution is carried out within a chain of symbols used for representing each chromosome. Note that each of these symbols represents a potential solution of the problem. This representation is defined in such a way that the transformation of the topology of the HEN can be easily associated to the real physical problem [12]. The approach presented here is based on the following assumptions: (a) All the heat exchangers are of counter-current type. (b) The heat capacities and the heat transfer coefficients are independent of flow temperature variations that occur along the heat exchangers. (c) There is no heat transfer across the HEN pinch point. (d) Flow stream splitting are not allowed.

3.1. Codification of the HEN structure Taking into account that independent variables used for the optimization are allocated to each heat exchanger, they will determine the topology of the HEN itself. Therefore, each unit (heat exchanger) can be identified by an appropriate integer value. This number represents a connection between a cold and a hot flow stream. Moreover, it is important to remark that the order used to locate the heat exchanger in the network is important. In order to take into account this order a ‘‘level” concept is introduced as shown schematically in Fig. 1 [12]. For a given level, each cold flow stream is restricted to receive heat only from one hot flow stream. Even though all possible interconnections are established in an arbitrary way, however, the aforementioned rule must be respected all along the application of the procedure. In this way the structure shown in Fig. 1 can be represented in a nL  nk ‘‘incidence” matrix form. The nL rows of the matrix correspond to the level number while the column number, nk, corresponds to the number of a cold flow stream. Thus, a different from zero value k of the (i, j) couple of the incidence matrix indicates that the flow stream j receives heat from the hot flow stream k across the HEN level i. This matrix representation of the HEN is shown schematically in Fig. 2. Note that the codification of the incidence matrix corresponds to a direct representation of a HEN having a real physical sense. Afterwards, this representation is trans-


J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773 Level 1

Level 2

Level 3 1 Hot Streams 2

1 Cold Streams 2

Fig. 1. Level concept applied to a simple HEN containing five heat exchanges.

Level Number

3.3. Genetic operators

Cold Stream Number 1











There is no heat exchanger on hot stream 2 at level 1 Matching between cold stream 2 and hot stream 1 at level 3

Fig. 2. Heat exchanger interconnection incidence matrix.

formed into a linear series of nL  nk integers that will correspond to each gene of a chromosome. This transformation provides a more convenient form for representing the information that must be handled by the GA. For instance, according to the matrix given in Fig. 2 the HEN shown in Fig. 1 can be represented by the chromosome (1, 0, 2, 2, 1, 1).

The heat load distribution is coded using a chromosome written in a binary form, i.e., formed only by 0 and 1 (bits) having a finite length L. For instance, a HEN having six heat exchanger units can be represented by the following binary chain: (0010, 1111, 1000, 0101, 0010, 1101), where each sub-group represents the heat load of a heat exchanger at a given position expressed in binary form. Depending on the precision required, each sub-group can have different lengths. The decoding of each sub-group (heat load) is carried out by using the following expression: bup  blow  C 10 ; 2L1

A ¼ a1 a2 a3 a4 a5 a6


B ¼ b1 b2 b3 b4 b5 b6 :


Thus, for a given arbitrary position kp = 4 the cross-over operation will result in the following two new chromosomes (i.e., new solutions):

3.2. Codification of the heat load distribution

X i ¼ blow þ

There are two major genetic operators: cross-over and mutation [16]. The cross-over operator is the most important one of any GA. This operator randomly chooses two chromosomes of a given population and mix their genetic information across a given location, kp. This location is randomly selected within the length of a chromosome, according to a uniform probability distribution function defined in the closed interval [1, L  1]. In this way, two new chromosomes are created by exchanging the initial genetic information encountered between the position 1 and kp of the original ones. Consider for instance two chromosomes A and B having the same length (=6), represented by


where Xi is the value of the thermal load for the heat exchanger i, blow and bup are the minimum and maximum values that Xi can reach, respectively. The last term in Eq. (2), C10, is the decimal conversion of the binary number Xi and L is the total length (number of bits) of the binary chain. Even though in this work a binary representation of the load distribution has been used, the proposed technique can be easily extended to a more convenient form that consists of using a real number representation.

A0 ¼ b 1 b 2 b 3 b 4 a 5 a 6


B0 ¼ a 1 a 2 a 3 a 4 b 5 b 6 :


Each new chromosome generated in this manner results from a combination of the genetic information contained by the two initial chromosomes identified as ‘‘parents”. The main purpose of a GA consists in exploiting the genetic heritage by creating new genetic resources. It is important to note that the cross-over operation does not introduce a new variety of genetic information; in order to enrich the genetic information a mutation operation is necessary. The mutation operator corresponds to a random and local modification of one gene in a given chromosome. Therefore, this operation permits to enlarge the solution space by changing the genetic patrimony contained in a given population. This operation constitutes an efficient tool that allows local maxima or minima encountered during the solution searching process to be avoided. 4. Modeling HEN’s In order to reduce the consumption of thermal energy, i.e., to satisfy an overall energy target, a mathematical

QC 1

Δ Th (k )

Tho (k )

T 1



Cold Streams

QC 2

ΔT (k )


T (k ) i c


Thi (k )


X m (k )

ΔT i (k )


S 3


QH 1 T3T

T (k ) o c

ΔTc (k )

Hot Streams

J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

T4T QH 2


= starting temperature of flow stream m


= target temperature of flow stream m

Fig. 3. HEN modeling approach.

model of the HEN must first be developed. Before writing the equations that will allow the optimization problem to be treated, consider the simple HEN shown in Fig. 3. It is apparent that the optimization of this network will contain several constrained equations. The equations equal to zero represent the conservation of energy along the flow streams. These equations satisfy the inlet and outlet temperature boundary conditions. According to Fig. 3, the energy conservation equation can be expressed as X X m ðkÞ þ QðC;H Þm  Qm ¼ 0; ð5Þ where Qm represents the heat transported by the flow stream m, Q(C,H)m are the utilities (i.e., heat sources or sinks) required to satisfy the overall heat target and Xm(k) is the load of heat exchanger k located in the flow stream m. Moreover, in order to avoid a violation of the second law of thermodynamics and to avoid all possible unrealistic solutions, i.e., infinite heat transfer surfaces, a pair of constrained equations different from zero is necessary. These constrained equations are written as DT i ðkÞ ¼ T oh ðkÞ  T ic ðkÞ P DT min ; DT o ðkÞ ¼ T ih ðkÞ  T oc ðkÞ P DT min :


where DTmin is the minimum approach temperature of the heat exchangers. The supra indices i and o indicate the inlet and outlet temperatures respectively, h identifies the hot flow streams, c identifies the cold flow streams and k identifies a particular heat exchanger in the network. Since the main objective of the optimization problem consists in reducing the total energyP consumption, then the objective function is given by f ¼ QH , where QH are the external utilities required to satisfy the temperature boundaries imposed by the process itself. Therefore, the linear model program can be reduced to the following form:



Constrained to T ohðkÞ

T icðkÞ


X k ¼ Qm ;

6 DT min ;


T ihðkÞ  T ocðkÞ 6 DT min : For the HEN shown in Fig. 3, the required system of constrained equations is formed by ten inequalities and four equalities. It is obvious that the number of constrained equations increases for more complex problems involving several heat exchangers. For problems with many constraints, a genetic algorithm will lack finding just only one feasible solution during a first run. Therefore, many runs and an enormous computational time should be necessary before few feasible solutions will be found. In order to avoid this major drawback and to be able to treat an optimization problem submitted to several constraints as described by Eqs. (5)–(7), the system must be transformed in a non-constrained one. In order to satisfy this requirement, a penalty function is introduced in the objective function. In such a case, a model that is able to treat the proposed HEN can be expressed by X min F ¼ QH þ P ; ð8Þ where P represents the penalty function that has to be solved for each chromosome of a given population. This penalty function is evaluated as follows: P ¼ Ck 





8 > <0 DT ðkÞ UðkÞ ¼ 1  DT min > :

if DT ðkÞ P DT min ;


if DT ðkÞ < DT min :

The number of inequalities required by the driving forces at the inlet and outlet of each heat exchanger is given by N; thus, Ck is a penalty coefficient used to control the relative


J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

importance of each constraint. Note that the penalty function given by Eq. (8) does not completely eliminate the non-feasible solutions (bad solutions) from the initial solution space. This function, however, reduces their importance by conserving their genetic information that can be useful in improving the evolution process of the GA itself. This particular approach differs from others where relaxation functions are used for controlling the so called bad solutions that are completely rejected from the solution space [13]). Note that a rejection of non-feasible solutions will not contribute in improving the genetic heritage during the evolutionary process.

4.1. Structure of the proposed algorithm In the proposed algorithm the solution search process is carried out in two steps. The first step consists in finding the network topology that has the potential of recovering the maximum thermal power. Herewith, two different populations are initially generated: a population that contains a potential solution for the HEN topology itself and a population that contains a potential solution for the thermal load distribution. Starting from these two populations a ‘‘funnel model”, schematically shown in Fig. 4 is used. The funneling technique is applied as follows: each HEN topology is used as a funnel that forces the random load population distribution to flow through in such a way that an evaluation of the maximum recoverable energy, for each possible network is carried out. The maximum energy recovery obtained for each structure is gathered in memory and only the topology that has the best fitness is retained. Then, by applying the cross-over and mutation operators those chromosomes characterized by a best fitness are used to generate a new population. This new population is after-

wards submitted to the same funneling process and the whole procedure is repeated until the evolution satisfies a convenient stop criterion. After the best HEN topology is obtained, the second step consists of refining the heat load distribution in order to optimize the maximum energy recovery. Each new generation (population) of heat load distribution is submitted to cross-over and mutation in order to evolve toward a solution that corresponds to a maximum of energy recovery. Note that all the constraints are rigorously satisfied during the aforementioned evolutionary process. The proposed algorithm is summarized as: Step 1: Generate random initial populations of HEN topologies and heat load distributions Step 2: For each HEN topology in the population For each heat exchanger Randomly assign a heat load distribution End loop Select the best heat load distribution that satisfies the maximum energy recovery Save the best heat load distribution of the HEN End loop Step 3: Select the HEN topologies that satisfy the best heat load distributions Step 4: Crossover the selected HENs by using a probability pc Step 5: Mutate the offspring’s by using a probability pm Step 6: New population of HEN offspring’s Step 7: If the stop criterion is not satisfied then Go to step 2 Else Stop Step 8: Pick up the best HEN topology having the best heat load distribution

Heat Load Crhomosomes Heat Load Crhomosomes






010111...11101 2







Linear Form of the Incident Matrix

100101...10001 . . . 011101...10101 Topology Crhomosomes 211120

Funnel Model




Apply: Selection Cross-over & Mutation

Fig. 4. Schematics of the funnel modeling technique.


J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

For this case, the minimum approach temperature (pinch) is DTmin = 10 °C. The topology of the HEN predicted with the proposed algorithm is exactly the same as the given in the literature. Table 2 shows a comparison between the heat load distributions obtained using the present method and the solution given in Ref. [17]. With the exception of the external utility QH2, the solution obtained with the proposed method is quite close to that given in the Ref. [17]. Even though, small differences in the load distribution are observed, it is important to note that the overall heat recovery (QREC) is almost the same as that of the benchmark. Fig. 5 shows the final synthesis of the HEN which also corresponds to that proposed in the literature. We present a second benchmark problem and compare our results with those proposed by Lewin et al. [12]. The temperature boundary conditions as well as the heat capacity of each flow stream are given in Table 3. The minimum approach temperature for this benchmark problem is set to 10 °C. As it is shown in this table, the problem consists of five hot and five cold flow streams with different heat capacities and inlet and outlet temperatures that must be satisfied for the aforementioned minimum pinch temperature. The results are summarized in Table 4 and a comparison of the heat load distributions is given in Table 5. The topology of the network obtained with the present algorithm is shown in Fig. 6. Although we do not have the same structure as that given in Lewin et al. [12], we

4.2. Application of the proposed method In order to demonstrate the capability of the proposed algorithm, three HENs problems taken from the open literature have been selected as benchmark tests. Each case has been synthesized by using the method described above. During the validation process an equal number of chromosomes (N = 100) were randomly generated for the initial population. This number was conserved after each crossover operation. Since no particular convergence criterion was applied for stopping the algorithm, a maximum number of 500 generations was used in all the cases for quitting the program. The software was fully implemented in Visual Fortran running in a 2 GHz PC. The temperature boundary conditions as well as the heat capacity of each flow stream, for the first benchmark test are given in Table 1 [17]. Table 1 HEN boundary and heat capacity conditions – benchmark 1 Flow stream

Inlet temperature (°C)

Outlet temperature (°C)

Heat capacity (MW/°C)

Hot 1 Hot 2 Cold 1 Cold 2

250 200 20 140

40 80 180 230

0.15 0.25 0.20 0.30

Table 2 Comparison between the predictions of the proposed method with the data of benchmark 1 Load/Source

GA solution (MW)

Ref. [17] (MW)

X1 X2 X3 X4 X5 QC1 QC2 QH1 QH2 QC min QH min QREC

6.56 17.44 12.49 7.97 7.04 9.94 0.00 0.97 6.54 10.004 7.504 51.495

6.5 17.5 12.5 8.0 7.0 10.0 0.0 0.0 7.5 10.0 7.5 51.5



Table 3 HEN boundary and heat capacity conditions – benchmark 2 Flow stream

Inlet temperature (°C)

Outlet temperature (°C)

Heat capacity (kW/°C)

Hot 1 Hot 2 Hot 3 Hot 4 Hot 5 Cold 1 Cold 2 Cold 3 Cold 4 Cold 5

433 490 490 478 472 333 389 311 355 366

366 411 339 422 339 433 495 460 450 478

8.79 10.55 14.77 12.56 17.73 7.62 6.08 8.44 17.28 13.9

149.9 ºC

196.8 ºC 250.0 ºC

40.0 ºC 9.94 MW

80.2 ºC

Hot Streams

150.0 ºC

200.0 ºC

80.0 ºC


140.0 ºC

175.2 ºC

0.97 MW

180.0 ºC

20.0 ºC Cold Streams

6.56 MW

7.97 MW 181.6 ºC

17.44 MW

140.0 ºC

12.49 MW

Fig. 5. Optimized HEN (benchmark 1).

208.2 ºC

230.0 ºC 7.04 MW

6.54 MW


J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

Table 4 Summary of the results obtained for the benchmark 2 Maximum energy recovery: 5771.24 kW Minimum cold utility requirement: 942.86 kW Minimum hot utility requirement: 91.2 kW DTmin = 10 °C Cold utility used (kW)

Hot utility used (kW)

Number of heat exchanger units

Number of cold utility units

Number of hot utility units

5771.24 5666.3

942.86 1047.8

91.2 196

8 7

4 5

1 3

Table 5 Comparison between the predictions of the proposed method with the data of benchmark 2 Load/source

GA solution (kW)

Ref. [12] (kW)

X1 X2 X3 X4 X5 X6 X7 X8 Qc1 Qc2 Qc3 Qc4 Qc5 QH1 QH2 QH3 QC min QH min QREC

535 471.2 1625.2 1541.2 227 480.4 786.3 – 54 47.16 217.8 223 506 164 16.4 15.57 1047.8 196 5666.3

588.93 553.28 559.44 1641.6 39.69 173.07 698.12 1517.11 – 280.17 153.72 5.24 503.73 91.2 – – 942.86 91.2 5771.24

can see that the recovered energy compare well with his results. It must be pointed out that for this benchmark, Lewin et al. have formulated the optimization as a MILP problem and they have used a classical method (i.e., Simplex) for solving it. If the problem is not linear (i.e., MIN-

Table 6 HEN boundary and heat capacity conditions – benchmark 3 Flow tream

Inlet temperature (°C)

Outlet temperature (°C)

Heat capacity (MW/°C)

Hot 1 Hot 2 Hot 3 Hot 4 Cold 1 Cold 2 Cold 3 Cold 4 Cold 5

327 220 220 160 100 35 85 60 140

40 160 60 45 300 164 138 170 300

0.10 0.16 0.06 0.40 0.10 0.07 0.35 0.06 0.20

LP), then the Simplex optimization method cannot be used. However, for the approach proposed in this paper this is not the case because the algorithm is able to handle any kind of problems, whatever they are linear or not. The third case tested consists of four hot and five cold flow streams (see Table 6). This case is given in Lewin [18]. It is important to remark that this author has optimized the HEN using a GA for the topology and the Simplex method for the heat load distribution. Note that for this case the minimum approach temperature (DTmin = 26 °C) is much higher than that used in the former cases. The comparison of the result obtained by using the proposed method with that given in Lewin [18] is shown in

372º C

433 ºC

54 kW

490 ºC

385.6 º C

490 ºC 459.2 º C

472 ºC

Cold Flow Streams

367.5 º C

422 ºC

339 ºC

535.0 kW

1625.2 kW

227.0 kW 480.4 kW

471.2 kW

1541.2 kW

403.2 ºC

433 ºC 164 kW

468 ºC

16.4 kW

15.6 kW

333 ºC

389 ºC 366.8º C

450 ºC 478 ºC

223.0 kW

411 ºC

339 ºC

506.0 kW

786.3 kW

460 ºC

218.8 kW

339.8 º C

478 ºC

495 ºC

47.2 kW

353.7 ºC

366 ºC

449.0 ºC

477.0 º C

Fig. 6. Optimized HEN (benchmark 2).

311 ºC

355 ºC 366 ºC

Hot Flow Streams

Ref. [12] Proposed algorithm

Recovered energy (kW)

J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773 Table 7 Comparison between the predictions of the proposed method with the data of benchmark 3 Load/source

GA solution (MW)

Ref. [18] (MW)

X1 X2 X3 X4 X5 X6 X7 X8 QC1 QC2 QC3 QC4 QH1 QH2 QH3 QH4 QH5 QC min QH min QREC

3.93 17.15 2.48 8.64 1.60 2.94 19.94 1.63 6.27 0.96 3.43 24.92 0.00 1.87 1.40 1.18 23.36 35.58 27.87 58.31

4.6 15.1 5.0 3.1 20.0 2.3 3.4 6.2 4.1 0.0 1.6 28.6 0.0 2.1 0.0 1.6 22.7 34.2 26.5 59.7


must manipulate several high dimensional tables that necessitate both high computational power and a large memory space. 5. Sensitivity analysis

Table 7. Fig. 7 shows the predicted topology for this HEN; even though the target is satisfied, the proposed topology differs from that recommended by Lewin. It is important to remark that due to the large number of constrained equations required for solving this network, the use of a standard technique such as Simplex becomes quite cumbersome. The degree of difficulty increases with increasing the number flow streams. In such a case, the Simplex method

The first benchmark case presented in Section 4 (see Fig. 5) was also used in order to study the sensitivity of the proposed algorithm with respect to the value of the penalty constant Ck introduced in Eq. (9). Moreover, the same HEN was synthesized by changing the cross-over (Pc) and mutation (Pm) probabilities. For a constant value of these probabilities (Pc = 10% and Pm = 1%), the constant Ck was changed from 10 to 10000; the overall thermal power target as a function of the number of generations is shown in Fig. 8. It is apparent that the convergence rate increases very quickly with increasing the value of Ck. This observation is obvious because a large value of this constant helps the algorithm to limit the effect of those chromosomes (i.e., solutions) that do not satisfy the constraints imposed to the HEN. Even though a large number of different simulations were carried out using the proposed technique, we are not able to recommend appropriate values for Ck. In addition some care should be taken when Ck is too big, because it may require starting the optimization procedure from a very large initial population space. Eq. (9) has also been used for evaluating the overall ‘‘quality” (fitness) of each population. Fig. 9 shows the penalty calculated with Eq. (9) as a function of the number 127.6 ºC

102.8 ºC

327.0 ºC

40.0 ºC 166.0 ºC

220.0 ºC

160.0 ºC 0.96 MW 117.1 ºC

192.8 ºC

143.8 ºC

60.0 ºC

220.0 ºC 3.43 MW 107.3 ºC

117.1 ºC

160.0 ºC

45.0 ºC 24.92 MW

300.0 ºC

100.0 ºC

35.0 ºC Cold Flow Streams

19.94 MW 113.8 ºC

91.0 ºC

137.0 ºC 164.0 ºC

3.92 MW

1.63 MW

1.60 MW

1.87 MW

134.0 ºC

138.0 ºC

85.0 ºC

17.15 MW

1.40 MW

101.3 ºC

60.0 ºC

150.3 ºC

2.94 MW

2.48 MW

170.0 ºC 1.18 MW

183.2 ºC

140.0 ºC

8.64 MW

Fig. 7. Optimized HEN (benchmark 3).

300.0 ºC 23.36 MW

Hot Flow Streams

6.27 MW


J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773


and increasing the cross-over probability seems to have a limited effect on the overall genetic quality of the population. Thus, it can be concluded that a small value of the mutation probability is necessary for generating chromosomes having new genetic characteristics, while the value of the crossover probability does not have a significant impact on the overall performance of the proposed technique.

Target 51.5 MW

51.6 51.4 51.2

10000 1000






6. Conclusions


A methodology based on the use of a genetic algorithm for synthesizing and optimizing large heat exchanger networks has been presented. The proposed technique allows both the topology and the heat load distribution to be solved simultaneously. The major drawback of solving a huge number of constrained equations is avoided by introducing a penalty function. Three benchmark cases taken from the open literature have been selected for validating the proposed methodology. For relatively small HENs the GA is able to reproduce quite well the results obtained by using conventional optimization techniques. For very large networks, however, some differences were observed. In all the cases studied, the proposed method is able to satisfy the overall heat power target. A sensitivity analysis has shown that increasing the importance of the penalty function considerably accelerates the optimization procedure. In addition, even though the performance of the technique seems to be strongly affected by mutation probability it is quite insensitive to the cross-over probability. Large HENs, which are almost impossible to be solved by using standard optimization techniques, can easily be treated by using the proposed GA method. In this work, the objective function is focused only on the Maximum Heat Recovery (MER) and do not tackle simultaneously the problem of reducing the number of units used to satisfy the MER, or the total cost involved. Since the knowledge of the MER may not be sufficient in the decision making process, then it is necessary to perform a multi-objective optimization in order to find solutions that simultaneously satisfy maximum heat recovery, minimum number of units and

50.2 50.0 1









NUMBER of GENERATIONS Fig. 8. Performance of the proposed method as a function of the penalty coefficient Ck (benchmark #1, Pc = Pm = 10%).

of generations for a constant value of Ck = 10,000 and different cross-over and mutation probabilities. In Fig. 9a a cross-over probability of 50% and a mutation probability of 10% were used for the first 500 generations; thereafter, the cross-over probability was kept constant while mutation probability was reduced from 10% to 1%, keeping the cross-over probability constant. In general, it is observed that reducing the mutation probability considerably reduces the value of the penalty function. Thus, a high mutation probability tends to destroy the genetic quality of the initial population. However, as was stated before, mutation is a very convenient operator required for generating new solutions when various local minima or maxima are encountered in an optimization problem. Fig. 9b shows the same case treated by using an initial mutation probability of 10% and a cross-over probability also of 10%. This probability was kept constant for the first 500 generations after which it was increased to 80% by maintaining the mutation probability constant. It is observed that keeping the mutation probability constant 1.4


Pm = 10%

Pm = 1%






0.8 0.6 0.4 0.2

Pc = 80%

Pc = 10%

1.0 0.8 0.6 0.4 0.2


0.0 1

100 200 300

400 500 600 700 800 900 1000



100 200 300

400 500 600 700 800 900 1000


Fig. 9. Penalty vs. number of generations (a) Pc = 50%; (b) Pm = 10%.

J. Dipama et al. / Applied Thermal Engineering 28 (2008) 1763–1773

minimum cost objectives. This corresponds to a search for the set of Pareto optimal solutions with respect to theses objectives. Work is under way based on the multi-objective optimization aspect of complex thermal systems [19]. Acknowledgements This work was funded by the discovery grant RGPIN 41929 of the Natural Sciences and Engineering Research Council of Canada (NSERC). References [1] B. Linnhoff, E. Hindmarsh, The pinch design method for heat exchanger networks, Chem. Eng. Sci. 38 (5) (1983) 745–763. [2] S.A. Papoulias, I.E. Grossman, A structural optimization approach in process synthesis, III. Total processing systems, Comput. Chem. Eng. 7 (1983) 723. [3] C.A. Floudas, A.R. Ciric, I.E. Grossmann, Automatic synthesis of optimum heat exchanger network configurations, AIChE J. 32 (1986) 276. [4] A.K. Saboo, M. Morari, R.D. Colberg, RESHEX: an interactive software package for the synthesis and analysis of resilient heat exchanger networks: I. Program description and application, Comput. Chem. Eng. 10 (1986) 577. [5] L.E. Grimes, M.D. Rychener, A.W. Westerberg, The synthesis and evolution of networks of heat exchangers that feature the minimum number of units, Chem. Eng. Commun. 14 (1982) 339–360. [6] B. Sagli, T. Gundersen, T. Yee, Topology traps in evolutionary strategies for heat exchanger network synthesis, in: H. Bussemaker, P. Iedema (Eds.), Process Technology Proceedings, vol. 9, 1990, pp. 51–58. [7] M.A. Duran, I.E. Grossmann, Simultaneous optimization and heat integration of chemical process, AIChE J. 32 (1) (1986).


[8] T.F. Yee, I.E. Grossmann, Simultaneous optimization models for heat integration, II. Heat exchanger networks synthesis, Comput. Chem. Eng. 14 (10) (1990) 1165–1184. [9] X.X. Zhu, B.K. O’Neill, J.R. Roach, R.M. Wood, A method for automated heat exchanger network synthesis using block decomposition and non-linear optimization, Chem. Eng. Res. Des. Part A 73 (11) (1995) 919–930. [10] C. Stair, E.S. Fraga, Optimization of unit operating conditions for heat integrated process using genetic algorithms, in: Proc. 1995 Int. Chem. Eng. Research Event, 1995, pp. 95–97. [11] I.P. Androulakis, V. Venkatasubramanian, A genetic algorithm framework for process design and optimization, Comput. Chem. Eng. 15 (4) (1991) 217–228. [12] D.R. Lewin, H. Wang, O. Shalev, A generalized method for HEN synthesis using stochastic optimization, I general framework and MER optimal synthesis, Comput. Chem. Eng. 22 (10) (1998) 1503– 1513. [13] M.A.S.S. Ravagnani, A.P. Silva, P.A. Arroyo, A.A. Constantino, Heat exchanger network synthesis and optimisation using genetic algorithm, Appl. Therm. Eng. 25 (2005) 1003–1017. [14] A. Pariyani, A. Gupta, P. Ghosh, Design of heat exchanger networks using randomized algorithm, Comp. Chem. Eng. 30 (2006) 1046– 1053. [15] J.H. Holland, Adaptation in Natural and Artificial Systems: An Introduction Analysis with Applications to Biology, Control and Artificial Intelligence, University of Michigan Press, 1975. [16] D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, 1989. [17] R. Smith, Chemical Process – Design and Integration, John Wiley, 2005. [18] D.R. Lewin, A generalized method for HEN synthesis using stochastic optimization – II: The synthesis of cost-optimal networks, Comput. Chem. Eng. 22 (10) (1986) 1387–1405. [19] M. Benali, A. Hammache, F. Aube´, J. Dipama, R. Cantave, Dynamic multiobjective optimization of large-scale industrial production systems: an emerging strategy, Int. Energy Res. 31 (12) (2007) 1202–1225.