A Simple Method for Finding Optimal Paths of Hot and Cold Streams inside Shell and Tube Heat Exchangers to Reduce Pumping Cost in Heat Exchanger Network Problems

In this paper, a simple method is presented for the synthesis and retrofit of heat exchanger networks (HENs) by considering pressure drop as well as finding proper path of streams inside heat exchangers (HEs) to reduce the pumping cost of network. Generally, HEN problems lead to MINLP models which have convergence difficulties due to the existence of both continuous and integer variables. In this study, instead of solving these variables simultaneously, a combination of Genetic Algorithm (GA) with Quasi Linear Programming (QLP) and Integer Linear Programming (ILP) was used for solving the problem. GA was used to find optimal HENs structure and streams paths, whereas continuous variables were solved by QLP. For the retrofit of HENs, a modified ILP model was used. Results show that the proposed method has the ability to reduce the cost of annual pumping due to considering optimal paths for streams in the HEs compared to the literature.


Introduction
According to the onion diagram provided by Smith 1 to obtain optimal operating conditions of industrial units and assigning two stages of this diagram for heat recovery systems as well as hot and cold utility management, the importance of synthesis as well as retrofit of heat exchanger networks (HENs) is significantly clear. The study of synthesis and retrofit procedures for HEN problems has been the subject of numerous research, due to the large savings that can be achieved in utility and also electrical consumption costs 2,3 . These different approaches can be broadly classified into three groups: pinch analysis, mathematical programming and stochastic algorithms alone or combined with the two previous methods.
The first group consists of pinch analysis method based on thermodynamic and graphical techniques, which was first seen in the work of Linnhoff and Hindmarsh 4 , and Tjoe and Linnhoff 5 . Linnhoff and Hindmarsh 4 used a simple method and manual calculations for the synthesis of HENs in order to find the maximum amount of heat recovery. Despite the simplicity of the method, due to the need for heuristic choices in the heat exchanger configuration stage, it was insufficiently reliable to achieve the final optimal structure. Also, this method was unsuitable for designing large HENs due to it being time consuming. Tjoe and Linnhoff 5 introduced a two-step method for modifying the HEN. In the first step, the retrofit targeting was performed, while in the second step, modifications were made on the results obtained from the first step. One of the major drawbacks of this work was the failure to consider any general rules for the area distribution in the network at the design stage. In general, the most important disadvantage of these methods is that they do not accurately capture the three way tradeoff between exchangers, heat exchange area, and energy use. Furthermore, using these methods for HEN synthesis and retrofit of large industrial cases is particularly tedious because they entail a great amount of manual development.
The second category is based on mathematical programming approaches found in the works of Grossmann and Sargent 6 , Floudas et al. 7 , Ciric and Floudas 8 and Yee and Grossmann 9,10 . The models derived from these methods were non-linear programming (NLP) and/or mixed integer non-linear programming (MINLP), which were divided into two parts: sequential and simultaneous approaches.
Floudas et al. 7 divided the problem into two sub-networks with different temperature intervals, which were solved sequentially using a sequential method for HEN synthesis. Because of this segmentation and the solution method used, it failed to reach the overall optimal solution. To overcome this shortcoming, Yee and Grossmann 9 used a simultaneous method without decomposition for the synthesis of HENs. Their approach was based on the MINLP model. One of the features of the model presented in their work was that it did not depend on the prediction of the pinch point for the subnetworks, nor was it dependent on constant approximation of the fixed temperature approaches. The features of this method include the independence of the prediction of the pinch point for the subnetworks as well as the constant approximation of the fixed temperature approaches. The only simplification used in this model was the isothermal mixing assumption. A similar model was proposed by Ciric and Floudas 11 for the synthesis of HENs using a hyper-structure. The generalized Bender's Decomposition (GBD) method was used to solve the MINLP model obtained for hyper-structure. Yee and Grossmann 10 proposed a two-stage model to retrofit HENs. In the first step, they used a mixed integer linear programming (MILP) model to estimate the economic feasibility of the retrofit, and in the second step, they used an MINLP model to find the best HEN. To consider additional details, they needed to apply binary variables to the MINLP model, which made it more complex.
It has been shown that the MINLP models can yield better structures than sequential models in synthesis as well as retrofit of HENs. But due to complexity of these models, they grow to an unmanageable size even for medium-sized problems, see Floudas 2 . To overcome these obstacles, Asante and Zhu 12,13 , and Zhu and Asante 14 took advantage of both the pinch-based and mathematical optimization methods in their proposed model. In their method for achieving proper thermal recovery, several MILP models were solved. After finding the best HEN structure, an NLP model was employed to minimize the network›s annual cost. Although this model was a success compared to the previous MILP models, it still had relatively poor performance for large-scale problems.
As reported by Furman and Sahinidis 15 , solving HENs based on MINLP model, either sequentially or with a stage wise-simultaneous formulation, is NP-hard, which causes increased computation time exponentially with problem size when using deterministic approaches.
To overcome the problems of the aforementioned approaches, the use of stochastic approaches in both synthesis and retrofit of HENs have been widely considered by researchers due to their flexibility and efficiency 16 . Therefore, the third group includes methods based on stochastic and heuristic algorithms alone or alongside previous methods. The use of simulated annealing (SA) for retrofit of HENs 17 and genetic algorithm (GA) 18,19 , tabu sreach (TS) 20 , combination of GA and SA 21 , differential evolution algorithm (DEA) 22 and particle swarm optimization (PSO) 23 for the synthesis of HENs, are examples of these efforts.
Although heuristic algorithms are the best for dealing with discrete variables, they are very slow in finding optimal values for continuous variables. Therefore, synthesis or retrofit of HENs may take a long time even for smaller cases. It seems that the heuristic algorithms are suitable for structural optimization (i.e., for handling binary variables) because of their discrete nature and strong ability for complete search of the solution space without being trapped at local optima 18,24 .
Therefore, the combination of stochastic algorithms with pinch-based methods or mathematical optimization in recent years has attracted much attention by researchers to solve MINLP problems in the design or modification of chemical processes, like 25,26 .
Rezaei and Shafiei 24 presented a new method for HENs retrofitting by coupling GA, NLP and integer linear programming (ILP). In this work, discrete variables (i.e., configuration of HEN) were carried out by the GA, and continuous variables were handled using a modified NLP formulation for maximum energy recovery. At the end, the minimum investment cost of modifications determined by an ILP model. Sreepathi and Rangaiah 27 used a new method for retrofit of HENs using IDE. In this study, several exchanger reassignment strategies are proposed. All of these strategies were then examined by considering a single objective function and using IDE. Finally, a reassignment strategy, which includes practical considerations was developed and used in multi-objective optimization on the retrofit of HENs. Leandro et al. 28 presented a meta-heuristic two-level method based on SA and Rocket Fireworks Optimization for the synthesis of multi period HENs. To improve solutions quality, they used a new post-optimization scheme, which was coupled with the methodology. In their work, they did not consider the isothermal mixing assumption. To calculate the total annual cost of the HENs, however, maximum required areas among all periods were considered. Aguitoni et al. 29 considered a method with two levels for HENs synthesis. GA used to optimize discrete variables at upper level and lower level heat loads, along with stream split fractions of HENs handled by DEA.
However, in most of the methods presented, the pressure drop issue was not considered in the synthesis and retrofit of HENs. For the first time, Polley et al. 30 , and Polley and Panjeh Shahi 31 showed that failure to consider the pressure drop in the synthesis and retrofit of HENs leads to serious problems in the final network. The excess of pressure drop in the final HEN over the maximum pressure drop allowed for the streams, will be one of the problems that results in the rejection of the final network. On the other hand, if the flow pressure drop is less than the minimum permissible level, it will result in excessive use of the heat transfer area of the required value. Due to the importance of this topic, the researchers were eager to consider the pressure drop of the streams in the synthesis as well as retrofit of HENs, mentioned in 32-34 based on pinch technology, 35-37 based on mathematical programming topics, and 38 based on a combination of stochastic algorithms and mathematical optimization methods.
Nie and Zhu 32 proposed a decomposition strategy, which first identified HEs that required additional area using a unit-based model. Next, for the structure of the previous step area distribution, shell arrangements and use of heat-transfer enhancement were optimized. At the same time, HEs that required no additional area were also optimized using simple models. Finally, the HEs were regrouped after optimization. The calculations continued until no change in members of each stage were observed. Panjeshahi and Tahouni 34 study the association of pumps and compressors costs together with the required additional area and operating cost using the pinch technology to solve the debottlenecking of the HEN of crude oil pre-heat train. Silva and Zemp 35 considered an optimization model for the retrofit of HENs, which minimized the additional area, and satisfied the available pressure drops of streams. In their method, the traditional HEN retrofit procedure was modified so that the problem was described as a non-linear model. The additional area required for the new HEN as well as available pressure drop were estimated based on economical optimization. Ravagnani and Caballero 37 proposed an optimization model based on area, energy, and pumping costs. Their algorithm was made of two distinct MINLP models. Soltani and Shafiei 38 proposed a new method for retrofitting the HENs by considering the pressure drop of streams inside the HEs, taking into account the constant heat transfer coefficients, and using the combination of linear programming (LP), ILP and GA. Behrouzrand and Soltani 39 used the method presented in reference 38 along with the use of ASPEN HYSYS simulator, succeeded in considering the temperature variations of variables, such as heat capacity and phase change of streams inside the HEs in the HENs synthesis.
Some of the important work done on the synthesis and retrofit of HENs in recent years are the works of Akpomiemie and Smith 40 , Tarighaleslami et al. 41 , Bao et al. 42 and Pavão et al. 43 In the work of reference 40 , a seven-step method based on NLP model was used to retrofit HENs by considering suitable pressure drop along with increasing heat transfer coefficients for streams. Tarighaleslami et al. 41 presented a new utility HEN design procedure for total site heat integration. The unified method they used applied more firm restrictions on the utility HEN than previous conventional methods. Bao et al. 42 developed an optimum-protection strategy called random walk algorithm with compulsive evolution for the synthesis of HENs. Their method had strong global and local search capabilities. To reduce computational time in that work, a leader-follower optimization method was considered. Pavão et al. 43 presented a new framework based on an extended superstructure model and a meta-heuristic approach for HENs retrofit. The mentioned structure contained features that are not available in most HEN studies based on mathematical programming, which can be pointed out for streams sub-splitting, partial mixing, serial HEs at a stream branch, and the inclusion of heaters/coolers at intermediate positions in all streams.
In all the work done, the paths of hot and cold streams within the shell and tube HEs were pre-determined. Although choice of flow paths, in addition to pressure drop considerations, depends on some mandatory considerations such as acidity and corrosion of the fluid type, optimization of these paths in each heat exchanger of the HEN reduces the cost of pumping flows, and thus reduces the annual costs of the final HEN.
Therefore, this study presents a simple method for the synthesis as well as retrofit of HENs by considering the pressure drop of streams and determining the appropriate path (pipe or shell path) for the hot and cold streams in each exchanger in order to reduce the costs of pumping streams in the network.
In the presented work, a combination of GA as one of the efficient stochastic algorithms and quasi-linear programming (QLP) method were used for the synthesis of HENs. The network structure as well as the optimal paths of the hot and cold streams inside the HEs are obtained by GA. The QLP model consisted of an outer search loop and an LP model. The outer search loop was used to find the optimal values of nonlinear variables (such as split ratio of streams and minimum approach temperature difference), and the LP model was used to find the optimal values of ordinary variables, such as the temperature of streams, heat load, and heat transfer area of each exchanger, as well as pressure drop of streams. Using the results of the QLP model and applying a modified ILP model of reference 38 for deciding whether to eliminate or reuse old HEN heat exchangers and pumps and/or introduce new ones into the new HEN, retrofitting of the HENs was carried out.

Problem description
Because HEN synthesis as well as retrofitting is considered as an optimization problem, the following information was assumed to be given in the problem definition: (i) inlet and outlet temperatures of process and utility streams, (ii) specific heat capacity (C P ), heat transfer coefficients (h), density (ρ), viscosity (µ), thermal conductivity (k), and mass or volumetric flow rate of the streams, (iii) the layout of initial HENs, heat transfer area of exchangers, and pressure drop of streams (for retrofit cases), and (iv) some economic data. In addition, a number of simplifying assumptions were considered in solving the HEN problems. These assumptions include: (i) using only shell and tube exchangers in the HENs, (ii) considering constant thermo-physical properties for streams, and (iii) ignoring piping costs in order to avoid complexity of the problem and ensure convergence. Objective function of such systems is to achieve minimum total annual costs and the minimum investment costs for the synthesis and retrofit cases, respectively. The total annual costs can include the cost of purchasing HEs and pumps, the cost of electricity for pumping, and the cost of hot and cold utilities. Also, the investment costs can be considered for purchasing new HEs and pumps, relocation of the existing HEs and pumps, and the costs of addition of area to the current HEs. The results of the problem solving included: (i) heat load of process and utility HEs, (ii) minimum approach temperature (ΔT min ), (iii) optimal flow rate and temperature of streams, and finally, (iv) determining the optimal path of hot and cold streams inside the HEs.

Methodology
As previously mentioned, HENs synthesis and retrofit is inherently an MINLP model, in which the values of discrete variables (i.e., layout of HEs and the connections of hot and cold streams to them) along with the continuous variables (such as: flow rate, temperature and pressure drop of streams, heat load and heat transfer area of HEs, load of pumping and …) should be calculated simultaneously. The nature of the MINLP models is such that no one reached the overall optimum solution of objective function in them 1,2 .
Thus, all the methods suggested by the researchers are only intended to reduce the complexity of the solution, as well as increasing the probability of reaching the final and better solutions. The method presented in this paper also follows this trend. Furthermore, it presents a very simple and practical way to find the optimal paths for hot and cold streams in both shell and tube sides of exchangers. In this paper, instead of solving simultaneously the continuous and discrete variables in synthesis and retrofit of HENs, a combination of GA (to find optimal structures of HENs) and QLP for the synthesis cases, and QLP + ILP model for the retrofit cases (to find the values of continuous variables and decision making for elimination or reuse of current exchangers and pumps and/or introducing new ones to the network as well as objective function) were used. The following sections describes these models (i.e., GA and QLP+ILP).

GA performance
To understand the performance of the GA and how to find the best HENs and optimal paths of hot and cold streams inside the shell and tube exchangers, the first step was to explain how the structure of HENs was generated by this algorithm. The most appropriate way to show the position of the exchangers within the network was node representation, seen in the works of Levin et al. 18 , Rezaei and Shafiei 24 , and Soltani and Shafiei 38 . For the use of the GA operators, the structure of each HEN (or chromosome) was divided into several genes. In this paper, the location of each exchanger in a structure of HEN was considered by numerical codes that constituted a victor which is known as the heat exchanger address victor (HAV). Fig. 1 shows a typical HEN structure with 3 genes and 6 exchangers ( Fig. 1a) and the corresponding HAV ( Fig. 1b) which was randomly generated by the GA.
Each HAV has the following properties: i) It is considered as a chromosome, and therefore consists of several genes. Since, there is no general rule to determine the number of genes, it should be selected heuristically, which depends on the size of the network under study. Also, each gene consists of 10 integer numbers.
ii) The first element of each gene indicates the number of HEs in that gene. A maximum of 3 and a minimum of 1 heat exchanger were considered in each gene to avoid complicating the networks created by the genes. If this number is one, like the first gene in Fig. 1, there is no need to split streams. If this number is 2, one of the hot or cold streams must be split into two branches, like the second gene in Fig. 1. However, if this number is 3, one of the hot streams and one of the cold streams must be split into two branches, like the third gene in Fig. 1.
iii) After this number, the next three numbers in each gene (i.e., numbers 2 to 4) indicate the path (or paths) of the hot stream (or streams) inside the exchanger (or exchangers) in that gene. The number 1 indicates the path of hot stream through the tube side, and number 2 indicates the path of hot stream through the shell side of each heat exchanger. For example, in the third gene, where there are three exchangers, these three numbers are 1, 2, and 1 respectively. These numbers indicate that hot streams in this gene flow inside the tube, shell, and tube of the exchangers (from left to right), respectively.
iv) The next three numbers in each gene (i.e., numbers 5 to 7) indicate the number of hot streams. For example, in the second gene, which had two exchangers, these three numbers are 2, 1, and 0, respectively. These numbers indicate that the second hot stream (i.e., H 2 ) enters the first exchanger (i.e., E2 in Fig. 1) and the first hot stream (i.e., H 1 ) enters the second exchanger (i.e., E3 in Fig. 1) in this gene. Because this gene contained only two exchangers, the final number must be 0. The same process is repeated for the last three numbers in each gene (i.e., numbers 8 to 10), but this time for cold streams.
Since the number of elements in each gene is considered to be 10, if the number of HEs defined in each gene is less than 3, then some of the elements in this gene will be zero, which is clearly shown in Fig. 1b.
It should be noted that in all HEN structures, a utility exchanger is considered at the end of each stream (cold utility for hot streams and hot utility for cold streams). For each hot stream, if the outlet temperature of the last HE (for example, the second hot stream outlet temperature (H 2 ) from the E6 in Fig. 1a) is greater than the target temperature for that stream (i.e. 2 H out T ), a cold utility exchanger will be used at the end of that hot stream. Also, for each cold stream, if the outlet temperature of the last HE is smaller than the target temperature for that cold stream, a hot utility exchanger will be used at the end of that cold stream. If the outlet temperature of the last HE for each stream is equal to the target temperature of that stream, the heat load of the related utility exchanger will be set to zero.
According to the description provided, various HEN structures can be easily produced by the proposed method. This high structural variation will increase the probability of finding an optimal overall objective function. As mentioned, in this study, GA was used to generate and search for optimal HEN structures. Fig. 2 shows the flowchart of this algorithm. Details of the various stages of its operation are given further herein.
After providing the required information and parameters (i.e., number of members in each generation, number of generations, the rate of crossover and mutation, number of hot and cold streams, and number of genes in each structure), by producing random numbers, GA creates initial population members in the same way as the HAV shown in Fig. 1b. In the following algorithm, the objective function of these members (i.e., HEN structures) are calculated according to the method described in the next section. Next generation members are then produced according to the objective function of each structure, as well as the help of reproduction, crossover and mutation operators. In the production of new generation members from previous generation structures, there is a high likelihood of employing members that have a better objective function. Selection of these members is carried out by the roulette wheel. This enables the production of better structures in future generations than in older ones by the GA. Next, the fitness of new generation structures are recalculated and this process continues until the predetermined number of generation is reached. In this paper, the number of initial population was set at 100 members, and reaching 100 generations was considered a termination condition. The performance of reproduction, crossover and mutation operators were as follows: Reproduction: based of this operator, 50 percent of the next generation of population members (i.e., chromosomes) are produced based on the two processes of elitism (5 %) and random selection (45 %), which are applied to current generation members.
Crossover: First, in this operator, some members or chromosomes (in pairs) from the current generation members are randomly selected as the parents. These chromosomes are then randomly divided into one or two parts. Finally, offspring is produced using the single-point and two-point crossovers. In this study, the rate of crossover was considered to be 50 %. Note that in this operator, the probability of selecting parents who have a better objective function value is higher. Fig. 3a shows how this operator works for single point crossover mode. As this figure shows, the first gene from the first parent was transferred to the second gene from the second parent, and vice versa to produce new structures (or offspring).
Mutation: After the new generation members were produced from the previous generation members using the reproduction and crossover operators, the mutation operator was applied to all them. A random number is generated for each new generation member. If this number is less than the rate of mutation defined for that generation, then mutation is applied to that member. Otherwise, the member will remain unchanged. The mutation rate for each generation is calculated from equation 1 25 : where Mu is the mutation rate, which is calculated in each generation. Mu min and Mu max are the minimum and maximum mutation rates allowed and their values are 10 % and 100 %, respectively, and C w and C b are the value of the cost of the worst and the best chromosome, respectively. The mutation operator was applied to the selected member or chromosome as follows: i) Firstly, a gene from that chromosome was randomly selected.
ii) The whole structure of that gene was removed and a new structure produced instead.
The performance of the mutation operator on a selected member or chromosome is shown in Fig.  3b. As this figure shows, only the entire structure of the selected gene (i.e., 3 rd gene) was removed, and a new structure produced instead. The structure of the other genes in the selected chromosome had not changed.

QLP and ILP performance
As obvious from the GA flowchart shown in Fig. 2, the QLP and ILP models were used to calculate the optimum values of variables (continuous variables in QLP and discrete variables in ILP), as well as objective function of synthesis and retrofit of HENs. In the following sections, the performance of the QLP and ILP models are explained. Also, their flowchart is shown in Fig. 4.

QLP procedure
To better understand the QLP procedure, it was necessary to write the energy and mass balance equations around the mixing and split points as well 2 Figure 2 Determine the initial parameters for GA as exchangers of each HENs. These equations included (i) energy balance for each exchanger; (ii) energy balance for mixing points; (iii) mass balance for splitters; (iv) monotonic decrease and increase of temperatures on hot and cold streams, respectively; and (v) minimum driving force temperature (ΔT min ) for hot and cold side supply, and exit temperatures in each exchanger. Some of these equations (i.e., energy balances for process HEs and mixing points) are nonlinear. Therefore, the resulting model will be an NLP model that is difficult to solve.
To reduce the complexity of the NLP, and make it a simpler QLP model, the following steps were taken.
1. The number of continuous variables in each structure was calculated based on the number of genes, hot and cold streams and HEs.
2. The variables were divided into two groups: 2-1 The first group variables included the temperature of the inlet and/or output of the HEs and their heat load.
2-2 The second group of variables, however, comprised the split ratio of the streams (y i ), and the minimum possible temperature difference between the hot and cold streams (ΔT min ) within the network.
Since the second group of variables caused the model governing equations to be nonlinear, an outer search loop combined with an LP model was used instead of solving all variables simultaneously.
In the outer loop, the values for the second group variables were considered from the corresponding intervals randomly. In this paper, the 3 In equation 2, Q k is the heat load of the HEs, ΔT is the approach temperature in hot and cold end of the exchangers, and S. F. is a scaling factor that must be chosen large enough so as not to affect the objective function, which is maximum energy recovery. This factor is included in the objective function to prevent some HEs from being pinched in ΔT min . This prevents using extra surface in HEs of the network 38 .
In the next step, using the guessed values of the second group variables and the calculated optimal values of the first group variables from the LP model, the heat transfer area and pressure drop of each heat exchanger in both its shell and tube side were calculated.
As mentioned in the structure of the HEN produced by the GA, for any heat exchanger, the hot stream may pass through the shell or tube. This situation is also true for the cold stream. So the equation of pressure drop for each heat exchanger are as below: For hot streams: 3.5 , , , These equations are valid for shell-tube HEs, streams without phase change, and turbulent flow conditions. These equations were first introduced by Poly et al. 30 to calculate the pressure drop of streams passing through the HEs, and later developed by other researchers such as 3,36,38 . In these equations, t and s stand for tube and shell, respectively, Δp t,k and Δp s,k are the pressure drop in the tube and shell side of k th exchanger, respectively, h i and h j are the film heat transfer coefficients of i th hot and j th cold stream, respectively. The parameters K t and K s are constants, which depend on the heat exchanger geometry and the physical properties of streams. Moreover, A k,i,j is the heat transfer area of k th exchanger between hot stream i and cold stream j which is given by equation 5: , , , , , In this equation, LMTD is the logarithmic mean temperature difference, and Q k,i,j is the heat load of k th exchanger between hot stream i and cold stream j. Furthermore, the parameters K t and K s are defined as equations 6 and 7 3 where μ, ρ, k and C P are the viscosity, density, thermal conductivity, and heat capacity of hot and/or cold stream, respectively, m is the mass flow rate, L tp is the tube pitch, D and D t are the internal and external diameters of the tubes, respectively, and D e is the equivalent diameter whose various relationships are listed in reference 3 . After calculating the pressure drop of each exchanger in the HEN, the pressure drop of the hot and cold streams is calculated according to the following method: 1-If the stream does not contain a splitter (or splitters), the total pressure drop is obtained by the sum of the pressure drop of all heat exchangers on that stream.
2-For cases where the stream contains a splitter, the pressure drop of the splitter section is considered to be the maximum pressure drop in one of the branches, like in Fig. 5.
Thus, until now, for the each HEN generated by the GA and taking into account the assumed values for the second group of variables, the pressure drop of the streams and the heat transfer area of each exchanger, as well as the temperature of the hot and cold streams inside the network, can be easily calculated by the QLP model and equations 2 to 7.
If the goal is synthesis of HENs, then the objective function is obtained from equation 8 39 . In these equations, C HE and C Pump are the annual costs for the exchanger and pump purchase, respec-tively, C utility is the total utility cost, and C elc is the electricity cost for the pump operation. Also, C HU and C CU are the hot and cold utilities costs, respectively, HU and CU are the load of hot and cold utilities, respectively, V · is the volumetric flow rate, η is the efficiency of the pump, ce is the unit cost of electricity, t is the operating time per year, I is the inflation rate, and the letters a to g as well as d* to f* are constant values that will vary by device type. But if the target is to modify or retrofit the HENs, then procedure enters ILP's model which is explained in the next section.

ILP procedure
In ILP model, the goal is to find minimum investment costs for converting a current HEN into a network proposed by the GA (i.e., new HEN). Therefore, some of the current network HEs may be discarded or reused in their original location or elsewhere by adding or reducing their heat transfer area. Also, due to the increase or decrease in pressure drop of streams in the new HEN compared to the current one, the location of some of the pumps may also change.
Note that, due to the presence of HEs in the new HEN, there will be a pressure drop in the streams. Therefore, no pump is excluded from the current HEN. In addition, new pumps may be purchased to overcome the additional pressure drop in the new HEN, if required. These new pumps will be placed in series next to the current HEN pumps. Therefore, the objective function of ILP model and its associated constraints can be as follows:   indicates that if the exchanger or pump from the current HEN is used for the same streams in the new HEN, they will not reassign costs. Otherwise, these costs will be applied to the objective function according to the economic data.
Note that the calculated objective function for either the HEN synthesis (i.e., equation 8) or the retrofit (i.e., equation 13) is the best possible amount with respect to the guessed values of the second group variables. Therefore, to calculate the general best possible objective function for the HEN structure generated by the GA, it is necessary to re-consider the new values for the second group variables according to the defined intervals for them. In the following, all steps are repeated from the beginning to obtain a new value for the objective function. By comparing the new value of the objective function with the previous one, new values for the second group variables are considered again. This goes on until the best values of the second group variables, as well as the general best value of the objective function for that HEN are obtained. After calculating the best optimal value of the objective function for the HEN, this value is sent to the GA. As illustrated in the flowchart of GA, Fig. 2, this algorithm will find the best HEN structure using this value as well as its operators (i.e., reproduction, crossover, and mutation). The performance of the proposed method is examined by studying two examples in the next section.

Illustrative examples
In this section, two case studies have been considered to evaluate the performance of the proposed method for the synthesis and retrofit of HENs with regard to pressure drop and variable flow path on the shell and tube sides in order to reduce the pumping costs.
Programming was done in MATLAB software, which contains some of the optimization functions and a number of codes for the process of the GA and decoding of the heat exchanger address victor to form the constraints of the optimization. But the run time in MATLAB environment is too long. Therefore, to solve this problem, the original codes were compiled to C language and final codes were loaded.
The number of generations as well as members of the initial population were considered 100 for both case studies. Furthermore, in all the HENs shown in the following figures, the underlined numbers indicate the heat load of HEs. Also, the time required to reach optimal results in each case study was approximately 2 and 3.5 hours.

Case study one: synthesis of HENs
This case study is a synthesis problem from Zhu and Nie 33 that includes four hot streams and five cold streams. In reference 33 , hot and cold streams are considered inside tube and shell, respectively. The thermodynamic information and physical properties of the streams as well as the cost data are given in Tables 1 to 3, respectively. Ta b l e 1 -Thermodynamic data for case study one The proposed Zhu and Nie 33 network had thirteen HEs with an annual cost of 3.081 M$ y -1 . For better comparison of the results between the proposed method and reference 33 , no split was considered in the design of HENs by GA. The best network obtained using the proposed method had 11 exchangers, as shown in Fig. 6. As may be seen from this figure, the path of hot streams in the four HEs shifted from the tube to the shell. The same thing happened in reverse for cold streams. Fig. 7 shows the performance of the GA in which the best and the average costs of each generation is plotted vs. the number of generations. The downward behavior of both curves in Fig. 7 confirms the correct performance of the GA. Non-repetition of the best answer after the 55 th generation also indicates the correct choice of 100 generation for this problem.  The results presented in these tables include the proposed method considering the optimum path of streams inside the HEs (case B), the proposed method when the hot streams flow inside the tube and the cold streams flow inside the shell, as in reference 33 (case A) and reference 33 . According to Table   4, the best solution obtained by the proposed method (i.e., case B) shows about 0.08 % and 9.5 % decrease in annual costs of the best HEN compared to case A and reference 33    to the last column of Table 5, show that finding the optimal paths of hot and cold streams within the tube and shell (case B) can be useful in reducing the pumping costs of HEN compared to the predetermined case of these paths (case A). According to the P/C system costs column in Table 4, this decrease is approximately equal to 0.002 %. The next important point is that the presented method in this paper (i.e., combining GA and QLP), yields better results despite its simplicity than the reference 33 which used NLP method. This may indicate high performance of the proposed method in the synthesis of HENs.

Case study two: Retrofit of HENs
This case study is a HEN for a pre-heat train of a crude oil distillation unit which has 6 hot streams and one cold stream. These hot streams heat the crude as a cold stream before entering the distillation column. As mentioned in 34 , the objective of this case study is the reduction in hot utility load, and is a debottlenecking problem after a 20 % increase in network throughput, shown in Fig. 8. As may be seen from this figure, all the hot streams flow inside the shell and the cold stream flows inside the tube. Tables 6 to 8 show stream data, thermo-physical properties, and cost data for this case study, respectively 34 . Fig. 9 shows the best HEN using the method proposed in this paper. As may be seen from this figure, it has used 9 process HEs instead of 6 process HEs in the current HEN (Fig. 8)     HEs hot streams, and in 5 remaining HEs the cold stream flows inside the tube. Tables 9 and 10 show uses of current exchangers and pumps of initial HEN (i.e., Fig. 8) in the new HEN (i.e., Fig. 9), as well as added area and purchased HEs and pumps. As Table 9 shows, not only have all the current network HEs been used in the new one, but only five HEs needed to be added certain area to the heat transfer surface. In addition, except number 1 and 7 HEs, which are used precisely in their original location, the reassignment costs of the remaining HEs have been added to the investment costs.
According to Table 10, the same reassignment costs for all pumps, except pump number 1 (belonging to the first hot stream H 1 ), were also considered.
In addition, as shown in Table 10, numbers 1, 4, and 6 of the hot streams as well as the cold stream, require new pumps to overcome excess pressure drop due to the HEN retrofit operations. These purchased pumps should be used in series with the initial pumps of the mentioned streams. Table 11 gives a comparison between two methods. As may be seen from this table, the results obtained using the proposed method in this paper represent a 0.57 % increase in annual utility cost savings, and more than a 9.5 % decrease in annual investment costs compared to reference 34 .
In order to investigate the effect of finding the optimal path of hot and cold streams inside HEs in modified HEN, a comparison was made between the costs of purchasing new pumps and also pump-  ing streams in the retrofitted HENs obtained by the presented method in this article (i.e., Fig. 9) and reference 34 , shown in Table 12. According to this table, by optimizing the flow paths of hot and cold streams inside the HEs, the costs of purchasing new pumps and reassigning them, as well as the annual costs of pumping in the optimal network obtained in this paper compared to reference 34 , were reduced by 0.6 % and 1.5 %, respectively.

Conclusion
This paper presents a simple and useful method of combining GA, QLP, and ILP models to consider the optimal paths of hot and cold streams inside shell and tube HEs in retrofit as well as synthesis of HENs thereby considering the pressure drop of the streams in order to reduce pumping costs. How to address HEs in the network using the victor method enables the generation of various structures by GA along with the variation of flow path for tube and shell in each exchanger for streams. In this victor, numerical codes 1 and 2, which represent the tube and shell paths for the hot streams, respectively, were used to determine the path of the hot streams inside the HEs. This provides the possibility to predetermine the paths of the streams in the shell and/ or tube in any exchanger if there are constraints such as corrosion and such. Instead of solving all continuous variables in an NLP model, dividing their variables into two groups and obtaining their optimal values using the QLP model (consists of an outer search loop + an LP model), although it would take time for the operation, increases the probability of reaching the final solution.
Comparison of the results with the references showed that the proposed method could obtain better results -albeit at a relatively lower annual pumping cost -than the references. This method works better when there are more HEs and streams in the HEN. An important point to be made in analyzing the results is that, these results were obtained using a simple method (combining GA, QLP, and ILP models) compared to other methods. Furthermore, despite the simplicity of the method, it shows relatively better results than in the references. Therefore, it can be used with a few modifications in large industrial cases including a variety of HEs, phase change of streams inside the HEs, pressure drop of streams, prohibited connections (due to economic or environmental reasons), etc. This is intended as future work for researchers. The only drawback of this method is that it takes a relatively long time to obtain the answer. This is common because of the nature of stochastic algorithms such as GA. However, this problem can be largely resolved by using a parallel computer system.