Synthesis of Cost-Optimal Heat Exchanger Networks Using a Novel Stochastic Algorithm and a Modiﬁed Stage-Wised Superstructure

: Facing the current energy structure urgently needs to be transformed, heat exchanger network (HEN) can implement heat recovery and cost reduction by the arrangement for heat exchanges between cold and hot streams. The plenty of integer and continuous variables involved in HEN synthesis cause the results to be easily trapped in local optima. To avoid this situation, the mechanism of accepting imperfect solutions is added in a novel algorithm called Random Walk Algorithm with Compulsive Evolution. However, several potential solutions maybe abandoned by accepting imperfect solutions. To maintain the global searching ability, and at the same time, protecting the potential solutions during the optimization process, the limitations of accepting imperfect solutions are investigated in this work, then a back substitution strategy and elite optimization strategy based on algorithm are proposed. The former is to identify and adjust the inferior individuals in long-term stagnation while the latter is to keep and perform a ﬁne search for the better solutions. Furthermore, a modiﬁed stage-wised superstructure is also developed to implement the ﬂexible placement of utilities, which efﬁciently enlarges the solution domain. The validation of strategies and model is implemented by three cases, the results are lower, with 2219 $/year, 1280 $/year, and 2M $/year than the best published result, revealing the strong abilities of the proposed method in designing more economical HENs.


Introduction
Heat exchanger network (HEN) synthesis can effectively realize heat recovery, reduce external utility consumptions, and improve economic benefits, so it has attracted extensive attention from many academic researchers and engineering workers in the past four decades and is widely applied in various fields including chemical process and crude refining industries [1,2]. HEN synthesis (HENS) is the heat integration between hot and cold process streams [3] with the aim of minimizing total annual cost (TAC) consisting of three sections i.e., utility consumption cost, fixed charge, and area cost of heat units.
Many methodologies for HENS have been developed and they can be broadly split into two categories, thermodynamic approaches and mathematical programming methods [4]. Pinch method [5] is a classical method used to evaluate optimal HEN using thermodynamics laws. In this method, the pinch point is first identified based on the given heat recovery approach temperature to predict maximal heat recovery, then the original HEN is decomposed according to the pinch point to predict minimal heat units. Subsequently, several methods deriving from pinch technology are proposed, including dual-temperature approach [6] and pseudo-pinch design method [7]. The pinch-based methods are well developed and demonstrated excellent efficiency [8][9][10]. promote structure mutation, which was implemented by accepting imperfect solutions (AIS) with a certain probability. The Metropolis acceptance rule i.e., the AIS behavior was an effective strategy applied to HENS problems [28,33], which was conducive to the local minimum avoidance. However, in RWCE, the AIS may cause individuals to get trapped in long-term evolutionary stagnation with the generation of numerous ineffective random walks, which is a limitation found in substantial case studies. We named this limitation L1. Additionally, AIS may also cause the original solutions before AIS to be replaced by imperfect solutions, so the original solutions cannot be exploited completely, which is another limitation called L2. Due to the two limitations, the RWCE can hardly evolve further and consequently miss the optimum.
The AIS plays a significant role in the applications of stochastic methods to HENS problems, the limitations of which should be investigated to help design better HENs. Therefore, the motivation of this work is to explore the specific limitations of AIS in HENS problems, based on which, two novel strategies including back substitution strategy and elite optimization strategy are proposed respectively to deal with L1 and L2. Then a novel algorithm combining RWCE and the two novel strategies is formed. Moreover, a modified SWS model is also designed with a flexible utility placement strategy, to embed more possible network structures. Finally, the novel algorithm and the modified SWS model are applied for the stream-splitting HENS problems.
In the following parts of the paper, the representation of the modified SWS model and the corresponding mathematical formulation are presented in Section 2, the novel algorithm is described in detail in Section 3, the investigative results of three illustrative cases are presented in Section 4 to verify the effectiveness of the proposed method, finally Section 5 provides the conclusions.

Mathematical Formulation and RWCE
The HENS process can be stated as the heat integration between a set of N H hot streams (to be cooled down to their required temperatures) and a set of N C cold streams (to be heated up to their required temperatures). The given conditions include heat capacity flow rate of each stream, overall heat transfer coefficients of each potential match, and cost data of heat exchangers and utilities.
A modified SWS (MSWS) model is proposed here, where utilities can be placed at the end of streams and each stream branch, so inner utilities can be obtained. Some similar superstructures can also be seen in the literatures [38,39]. Figure 1 depicts the comparison of SWS and MSWS based on a four-stream problem. Obviously, the MSWS includes the SWS. Moreover, the assumption of non-isothermal mixing [40][41][42] is used here. The simplification of isothermal mixing assumption requires that the outlet temperatures of all branches of a single stream are equal at each stage, which is a limitation to the generation of better candidate HEN structures.

Objective Function
The objective function of HENS is TAC: where A i,j,k is heat exchanger area, A HU,j and A HU,i,j,k are respectively the areas of hot utility at stream end and stream branch, Q HU,j and Q HU,i,j,k are the loads of hot utility respectively at stream end and stream branch. C F , C A , C CU , and C HU are respectively fixed charge of heat units, area cost coefficient, cost coefficient of cold and hot utility. E is an integer variable representing whether the heat unit exists (E = 1) or not (E = 0).

Objective Function
The objective function of HENS is TAC: , , +  , , , +   , , ,   +  , , , +  , , , +   , , , where , , is heat exchanger area, , and , , , are respectively the areas of hot utility at stream end and stream branch, , and , , , are the loads of hot utility respectively at stream end and stream branch. , , , and are respectively fixed Areas of heat exchanger can be calculated by Equation (2). Where LMTD represents logarithmic mean temperature difference and is calculated by Equation (5) when θ 1 is not equal to θ 2 , otherwise it is replaced by arithmetic mean temperature difference (AMTD) denned as Equation (6). Where T in i,j,k,h and T out i,j,k,h are inlet and outlet temperatures of a single exchanger at hot stream; T in i,j,k,c and T out i,j,k,c are inlet and outlet temperatures of a single exchanger at cold stream. Equation (7) represents the calculation of overall heat exchange coefficient K, where h is individual heat exchange coefficient of streams. The calculation for utility area is similar to Equation (2).

Constraints
The feasibility constraints are given as follows: (1) Overall heat balance for each stream where T in i and T out i are inlet and outlet temperatures of hot streams; W is heat capacity flow rate.
(2) Heat balance for each heat exchanger where sp i,j,k,h , and sp i,j,k,c are split ratios of hot and cold streams, representing the proportion of the heat capacity flow rate of a stream branch to the corresponding total stream.
(3) Heat balance for each stage where T out i,k and T out j,k are outlet temperatures of hot and cold stream at a certain stage. (4) Heat balance for the utility at stream branch T out where T in CU,i,j,k and T out CU,i,j,k respectively denote the inlet and outlet temperatures that a hot stream branch flows through a cold utility, T in HU,i,j,k and T out HU,i,j,k respectively denote the inlet and outlet temperatures that a cold stream branch flows through a hot utility.
(5) Heat balance for the utility at stream end Loads of heat exchangers and inner utilities are treated as optimization variables. If the stream cannot reach the target outlet temperature after each iteration, the utilities at stream end are used as additional energy to meet the requirements.
(6) Heat balance for each mixer (7) Mass balance for each mixer (8) Temperature constraints in stream branches (9) Minimum temperature approach constraints where ∆T min is the minimum approach temperature of heat units. T in CU and T out CU respectively denote the inlet and outlet temperatures of cold utility. The corresponding symbols for hot utility are T in HU and T out HU .

Methodology
A novel algorithm combining RWCE and two strategies to solve the limitations of AIS is designed to solve the complex MSWS model with multiple optimization variables. Several points presented as follows should also be paid attention to. In RWCE, the historical optimum of an individual means the best solution obtained by the individual in the past iterations. The current solution of an individual represents the solution obtained by the individual at the current iteration. The optimal solution indicates the best one in the historical optimums of all the individuals. The optimal TAC after AIS means the best TAC obtained since the individual accepts an imperfect solution at a certain iteration.

RWCE Method for HENS with Stream Splits
In consideration of MSWS model, an individual in RWCE represents a HEN consisting of heat exchanger loads, inner utility loads, and stream split ratios. The main evolution to implement is the "random walking" of each individual as presented in Equations (30)- (34). Q n,p,it+1 = Q n,p,it + (1 − 2 × r 1 ) × r 2 × S Q , n = 1, . . . , N, p = 1, . . . , NE QICU n,l,it+1 = QICU n,l,it + (1 − 2 × r 3 ) × r 4 × S ICU , n = 1, . . . , N, l = 1, . . . , NE (31) QI HU n,m,it+1 = QI HU n,m,it + (1 − 2 × r 5 ) × r 6 × S I HU , n = 1, . . . , N, m = 1, . . . , NE (32) SRH n,p,it+1 = SRH n,m,it + (1 − 2 × r 7 ) × r 8 × S SRH , n = 1, . . . , N, p = 1, . . . , NE (33) SRC n,p,it+1 = SRC n,p,it + (1 − 2 × r 9 ) × r 10 × S SRC , n = 1, . . . , N, p = 1, . . . , NE (34) where Q n,p,it is the initial heat exchanger load, Q n,p,it+1 is the heat exchanger load after walking. QICU n,l,it is the initial load of inner cold utility, QICU n,l,it+1 is the load of inner cold utility after walking. SRH n,m,it represents the initial split ratio of hot stream, SRH n,p,it+1 is the split ratio of hot stream after walking. S q , S ICU , and S I HU are respectively maximal walking step sizes of heat exchanger load, inner cold and hot utilities. S SRH and S SRC are respectively maximal walking step sizes of split ratios of hot and cold streams. r 1~r10 . are uniform random numbers distributed in the interval (0, 1). Lower bound Q min is set to handle integer variables to determine the generation and elimination of heat exchangers and inner utilities, which is described as Equation (35) (taking heat exchanger as an example). Q min is a parameter to control the existing of heat exchanger and its value is 10 kW in this paper. If Q n,p,it+1 is smaller than Q min , the corresponding heat exchanger will be eliminated, making the effective heat exchanger load Q n,p,it+1 equal to 0.
At the itth iteration, the initial solution for individual n can be expressed as X n,it . After walking, the new effective solution X n,it+1 is generated. Then the TAC(X n,it+1 ) is compared with TAC(X n,it ), and the lower one is selected. If TAC(X n,it+1 ) is not smaller than TAC(X n,it ), X n,it+1 can also be accepted with a certain probability δ. The specific operation is expressed as Equation (36), which is the mutation and selection criterion of RWCE.

Limitations of Accepting Imperfect Solutions
Two limitations of AIS are discussed in this section. L1: The individuals may not avoid local optimums or even walk to imperfect regions by AIS in the middle or late period, which will lead to the evolution stagnation and efficiency decrease. A nine stream problem listed in Table 1 was used to explain the L1. Figure 2 depicts the iteration of optimal TACs obtained by RWCE (N = 1) with the δ respectively set to 0 and 0.01. It is pronounced that the TAC corresponding to δ = 0 declines more slowly compared to that corresponding to δ = 0.01 in the early period. The final TACs corresponding to δ = 0 and δ = 0.01 were respectively 3,250,825 $/year and 2,934,609 $/year, so the AIS could accelerate the evolution process and improve the convergence precision. However, for δ = 0.01, the individual did not evolve (the TAC didn't decline) any more after about 9,500,000 iterations. By tracking the optimization process, the individual accepted an imperfect solution at the 9,601,186th iteration and the current TAC rose from 2,934,609 $/year to 2,939,384 $/year. Figure 3 displays the iteration of the historical optimal TAC and the optimal TAC after AIS from the 9,601,186th to the 9,999,186th iteration (after the 9,999,186th iteration, the two TACs keep constant). After AIS, the individual could not find better solutions for a long time and stopped at a TAC of 2,934,699 $/year which was inferior to the historical optimal TAC 2,934,609 $/year.     When the individual stopped evolving, its corresponding HEN structure was recorded to compare with its historical optimal structure before AIS. Figure 4 depicts the diagrammatic sketch of the comparison result. Figure 4a shows the best HEN structure before AIS, Figure 4b,c gives possible HEN structures when evolution stops after AIS. It can be found that Figure 4b has the same structure as Figure 4a, but their continuous variables are different. For instance, the heat loads of heat unit 1 in the two structures are respectively 100 kW and 90 kW. While Figure 4c has a structure different from Figure 4a, because the heat unit 2 is deleted in Figure 4c. Therefore, both the continuous variables distribution and the structure might be destroyed after AIS, making the individuals walk ineffectively for a long time and have difficulty in returning to their own best position.  When the individual stopped evolving, its corresponding HEN structure was recorded to compare with its historical optimal structure before AIS. Figure 4 depicts the diagrammatic sketch of the comparison result. Figure 4a shows the best HEN structure before AIS, Figure 4b,c gives possible HEN structures when evolution stops after AIS. It can be found that Figure 4b has the same structure as Figure 4a, but their continuous variables are different. For instance, the heat loads of heat unit 1 in the two structures are respectively 100 kW and 90 kW. While Figure 4c has a structure different from Figure 4a, because the heat unit 2 is deleted in Figure 4c. Therefore, both the continuous variables distribution and the structure might be destroyed after AIS, making the individuals walk ineffectively for a long time and have difficulty in returning to their own best position.  L2: The L2 has a negative effect on the convergence precision, which should also be paid attention to. After accepting an imperfect solution, the original solution before AIS is abandoned that cannot get further optimization, so the original solution is not exploited fully and then the precision declines.

Strategies
Removing the AIS behavior is obviously an ineffective approach to deal with the two limitations, because it can only keep the local optimums. Therefore, to enhance the search ability of RWCE, two novel strategies aiming at the limitations of AIS are designed in this section. The strategies are then combined with RWCE to form a novel stochastic algorithm named IRWCE.

Back Substitution of Optimums
To deal with the L1, a new strategy called back substitution of optimums is developed to compulsively make the individuals in long-term stagnation return to their own historical optimums to adjust their optimization direction in time.
The historical optimum of each individual is first recorded; then an individual condition monitoring mechanism is proposed, a counter is designed to count the evolution stagnation iterations of each individual after AIS. When the individual evolves, the counter resets; when the counter reaches a certain value (BSC), the corresponding individual is considered to be in a long-time stagnation; subsequently, the recorded historical optimum will be delivered to the corresponding individual to implement the back substitution. After back substitution, the count does not continue and becomes zero only if the corresponding individual evolves. If the counter reaches the value 2 × BSC, the back substitution of the recorded historical optimum to the corresponding individual will be performed again. The maximal iteration BSC between two back substitutions is called back substitution cycle. The strategy will be performed repeatedly based on the back substitu- L2: The L2 has a negative effect on the convergence precision, which should also be paid attention to. After accepting an imperfect solution, the original solution before AIS is abandoned that cannot get further optimization, so the original solution is not exploited fully and then the precision declines.

Strategies
Removing the AIS behavior is obviously an ineffective approach to deal with the two limitations, because it can only keep the local optimums. Therefore, to enhance the search ability of RWCE, two novel strategies aiming at the limitations of AIS are designed in this section. The strategies are then combined with RWCE to form a novel stochastic algorithm named IRWCE.

Back Substitution of Optimums
To deal with the L1, a new strategy called back substitution of optimums is developed to compulsively make the individuals in long-term stagnation return to their own historical optimums to adjust their optimization direction in time.
The historical optimum of each individual is first recorded; then an individual condition monitoring mechanism is proposed, a counter is designed to count the evolution stagnation iterations of each individual after AIS. When the individual evolves, the counter resets; when the counter reaches a certain value (BSC), the corresponding individual is considered to be in a long-time stagnation; subsequently, the recorded historical optimum will be delivered to the corresponding individual to implement the back substitution. After back substitution, the count does not continue and becomes zero only if the corresponding individual evolves. If the counter reaches the value 2 × BSC, the back substitution of the recorded historical optimum to the corresponding individual will be performed again. The maximal iteration BSC between two back substitutions is called back substitution cycle. The strategy will be performed repeatedly based on the back substitution cycle if the count continues.
Limited by the relatively small walking step size i.e., (1 − 2 × r 1 ) × r 2 × S Q , individuals can hardly walk with a long distance in the limited iterations. Hence, when the individual cannot evolve the structure after several operations of back substitution, a random perturbation is executed to achieve a long jump. The perturbation here is described as follows. Some heat loads Q n,best are selected randomly from the individual historical optimal HEN, which are then given perturbation as presented in Equation (37), finally the small heat units are also eliminated as Equation (35). The r 12 is a uniform random number distributed in the interval (0, 1), k is the perturbation factor. In this study, κ is directly set to 12. The perturbation may help realize long jumps by giving relatively big changes to the heat exchanger loads. However, this perturbation may also destroy the efficient matches of individual historical optimums, so it is executed only when several back substitutions of individual historical optimums cannot evolve the structure.
The working procedure of the back substitution strategy is depicted in Figure 5. The "M" represents the evolution stagnation iterations, "G" is the number of back substitution. The "G cr " is set to determine whether the perturbation can be implemented or not. The working procedure of the back substitution strategy is depicted in Figure 5. The "M" represents the evolution stagnation iterations, "G" is the number of back substitution. The " " is set to determine whether the perturbation can be implemented or not.

Elite Optimization
The L2 shows that original solutions before AIS may be replaced by imperfect solutions, so it cannot get complete optimization. Hence, an elite optimization strategy is proposed here to solve L2, where the population is divided into two sections, one section includes NP (NPV0.5N) elite individuals, the other section includes N-NP basic individuals. The elite individuals protect and optimize the solutions with better potentials i.e., elite solutions from being replaced by imperfect solutions, where elite solutions are selected in the optimization process of basic individuals.
The basic individuals evolve as the original RWCE, their historical optimums are sequenced according to TAC, where the top NP historical optimums with relatively lower TACs are selected as the elite solutions and delivered to elite individuals. An elite solution corresponds to an elite individual. If a basic individual updates its historical optimum which is better than the worst elite solution, the new generated optimum will replace the worst elite solution. Thus, the elite solutions are kept without being replaced by imperfect solutions and also have the opportunity to get further optimization.
For the elite solutions, a fine-search (FS) strategy [43] is introduced. The FS strategy can improve the precision in both the continuous and integer variables, which is described in the following. In Equation (38)

Elite Optimization
The L2 shows that original solutions before AIS may be replaced by imperfect solutions, so it cannot get complete optimization. Hence, an elite optimization strategy is proposed here to solve L2, where the population is divided into two sections, one section includes NP (NPV0.5N) elite individuals, the other section includes N-NP basic individuals. The elite individuals protect and optimize the solutions with better potentials i.e., elite solutions from being replaced by imperfect solutions, where elite solutions are selected in the optimization process of basic individuals.
The basic individuals evolve as the original RWCE, their historical optimums are sequenced according to TAC, where the top NP historical optimums with relatively lower TACs are selected as the elite solutions and delivered to elite individuals. An elite solution corresponds to an elite individual. If a basic individual updates its historical optimum which is better than the worst elite solution, the new generated optimum will replace the worst elite solution. Thus, the elite solutions are kept without being replaced by imperfect solutions and also have the opportunity to get further optimization.

Calculation Steps of the Improved RWCE
An improved RWCE (IRWCE) based on the original RWCE and the two new strategies is developed, the process of which is exhibited in Figure 7. (

1) Initialization The population including N-NP basic individuals and NP elite individuals is set. Each individual is initialized as Equations (44)-(48). Where
, , , , , , , , , , , , and , , respectively represents the initial heat exchanger load, inner cold, and hot utility loads, and split ratios of hot and cold streams; , , , , and are respectively the initial maximal heat exchanger load, inner cold and hot utility loads, and split ratios of hot and cold streams, they are set as 200 kW, 0 kW, 0 kW, 10 and 10 in this work.
, , = × ， = 1, … , , = 1, … , (2) Back substitution The back substitution strategy is used for basic individuals. The individual condition monitoring mechanism works to count the evolution stagnation iterations after AIS. If the individuals do not satisfy the requirements of back substitution, this step will be skipped. The strategy is executed based on the back substitution cycle and ends only if the corre-

Calculation Steps of the Improved RWCE
An improved RWCE (IRWCE) based on the original RWCE and the two new strategies is developed, the process of which is exhibited in Figure 7.

Results and Discussion
In this section, three different-sized cases available in the literature were solved with the proposed method IRWCE based on the MSWS model. The evaluation for the performance of algorithms depends on the precision and efficiency. While for HENS problems, the precision of the solution obtained within acceptable computation effort is the primary concern [39]. Therefore, the HENs obtained by IRWCE were compared with the previ-  (44)- (48). Where Q n,p,0 , QICU n,l,0 , QI HU n,m,0 , SRH n,p,0 , and SRC n,p,0 respectively represents the initial heat exchanger load, inner cold, and hot utility loads, and split ratios of hot and cold streams; Q max , QICU max , QI HU max , SRH max , and SRC max are respectively the initial maximal heat exchanger load, inner cold and hot utility loads, and split ratios of hot and cold streams, they are set as 200 kW, 0 kW, 0 kW, 10 and 10 in this work. Q n,p,0 = Q max × r 24 , n = 1, . . . , N, p = 1, . . . , NE

(1) Initialization The population including N-NP basic individuals and NP elite individuals is set. Each individual is initialized as Equations
QICU n,l,0 = QICU max × r 25 (45) QI HU n,m,0 = QI HU max × r 26 (46) SRH n,p,0 = SRH max × r 27 , n = 1, . . . , N, p = 1, . . . , NE SRC n,p,0 = SRC max × r 28 , n = 1, . . . , N, p = 1, . . . , NE (2) Back substitution The back substitution strategy is used for basic individuals. The individual condition monitoring mechanism works to count the evolution stagnation iterations after AIS. If the individuals do not satisfy the requirements of back substitution, this step will be skipped. The strategy is executed based on the back substitution cycle and ends only if the corresponding individual evolves its structure after back substitutions.
(3) Individual random walk The basic individuals walk according to the Equations (30)- (34), while the elite individuals walk by the Equations (39)- (43). The Equation (35) is used to deal with integer variables, the minimal heat loads for the basic and elite individuals are respectively Q min,b and Q min,e .
(4) Selection and mutation Equation (36) presents the specific implementation of selection and mutation. After random walking, the generated candidate solutions are compared with the initial ones, where the better ones are selected. If the candidate solution is not better than the initial one, the candidate solution will also be accepted with a certain mutation probability. The mutation probabilities for the basic and elite individuals are respectively δ b and δ e .
(5) Termination If the maximal iteration is satisfied, the iteration ends; otherwise go to step (2).

Results and Discussion
In this section, three different-sized cases available in the literature were solved with the proposed method IRWCE based on the MSWS model. The evaluation for the performance of algorithms depends on the precision and efficiency. While for HENS problems, the precision of the solution obtained within acceptable computation effort is the primary concern [39]. Therefore, the HENs obtained by IRWCE were compared with the previously reported ones on TAC that can describe the optimization precision objectively. The code was implemented using the Compaq Visual Fortran (Version 6.6) and executed on a Windows Server system with an Intel(R) Xeon(R) CPU E31235 @ 3.20 GHz.

Case Study 1
The case with the problem data listed in Table 2 have six hot streams and four cold streams. Ravagnani et al. [44] developed a method combining GA and Pinch Method to solve this case, finally a TAC 5,672,821 $/year was presented. The application of IRWCE to the case based on the MSWS model gave a TAC 5,586,395 $/year, the corresponding HEN structure is given in Figure 8 (split ratios are in parentheses). Table 3 gives the comparison of results; our solution is better than the reported ones and saved 2219 $/year as compared to that by Chen et al. [26]. An inner utility with 6723.5 kW at C4 can be seen in our structure, which cannot be achieved by the original SWS. It indicates that the MSWS enlarges the solution region and can give more satisfactory HEN.   The fixed charge and area cost exponent for the case are respectively zero and one, due to which the obtained HENs tend to include many heat units. To make the case more realistic, Huang et al. [40] provided a new fixed charge 8000 $/year, then they gave a result of 5,737,274 $/year. The modified case was then solved in some literatures. Huang and Karimi [46] gave a HEN with the TAC of 5,733,679 $/year, who used an improved outerapproximation algorithm without generating feasible starting points. Pavao et al. [33] also presented a TAC 5,715,026 $/year by employing a two-level algorithm. In this work, IR-WCE and MSWS found a network structure given in Figure 9, its TAC was 5,713,746 $/year and the computational time is 28,810.59 s, which was superior to all the reported ones. Table 4 presents the results comparison, the significant reduction in area leads to the  The fixed charge and area cost exponent for the case are respectively zero and one, due to which the obtained HENs tend to include many heat units. To make the case more realistic, Huang et al. [40] provided a new fixed charge 8000 $/year, then they gave a result of 5,737,274 $/year. The modified case was then solved in some literatures. Huang and Karimi [46] gave a HEN with the TAC of 5,733,679 $/year, who used an improved outer-approximation algorithm without generating feasible starting points. Pavao et al. [33] also presented a TAC 5,715,026 $/year by employing a two-level algorithm. In this work, IRWCE and MSWS found a network structure given in Figure 9, its TAC was 5,713,746 $/year and the computational time is 28,810.59 s, which was superior to all the reported ones. Table 4 presents the results comparison, the significant reduction in area leads to the cost advantage of our solution.

Case Study 2
The case involving nine streams was further investigated here, which was originally proposed by Linnhoff and Ahmad [47]. The problem data are presented in Table 1. Many scholars tested their methods by solving this case. Huo et al. [45] presented a modified SWS in which the stages with or without stream splits can be determined flexibly, then a two-level algorithm combining GA and PSO was utilized, finally reporting two HENs whose TACs were respectively 2,922,585 $/year and 2,925,634 $/year. Pavao et al. [39] presented a new SWS model with different placements of heaters and coolers, then a modified method combining SA and Rocket Fireworks Optimization was designed and applied to the case, annually giving a TAC 2,909,906 $/year whose structure included an inner utility.
The IRWCE and MSWS were applied to the case, BSC=2,000,000 was set, then the HEN in Figure 10a was obtained with the TAC of 2,922,234 $/year. When the BSC was set to 3,000,000, a TAC 2,921,233 $/year was obtained, the corresponding structure is presented in Figure 10b. Figure 10c displays the optimal HEN obtained by setting BSC to 2,500,000, whose TAC is 2,907,007 $/year within 8150.20 s and is 2899 $/year less as compared to that by Pavao et al. [39]. The selection of back substitution cycle has significant impact on the optimization process. When BSC is set too big, the individuals in long-term stagnation will not be adjusted in time; when BSC is set too small, the positive effect of AIS will be weakened.
The comparison of results is put in Table 5. Obviously, our optimal HEN is best, an inner utility (2556.8 kW) can also be seen at a split branch of C1. Furthermore, by observing the temperatures of C1 and HU, the match between C1 and HU is infeasible when C1 is heated to over 250 °C, so the HU at the end of C1 can only be 0 kW (C1 is directly be heated to the target temperature by the matches with hot streams) and over 5000 kW (the outlet temperature of C1 is lower than 250 °C).

Case Study 2
The case involving nine streams was further investigated here, which was originally proposed by Linnhoff and Ahmad [47]. The problem data are presented in Table 1. Many scholars tested their methods by solving this case. Huo et al. [45] presented a modified SWS in which the stages with or without stream splits can be determined flexibly, then a two-level algorithm combining GA and PSO was utilized, finally reporting two HENs whose TACs were respectively 2,922,585 $/year and 2,925,634 $/year. Pavao et al. [39] presented a new SWS model with different placements of heaters and coolers, then a modified method combining SA and Rocket Fireworks Optimization was designed and applied to the case, annually giving a TAC 2,909,906 $/year whose structure included an inner utility.
The IRWCE and MSWS were applied to the case, BSC = 2,000,000 was set, then the HEN in Figure 10a was obtained with the TAC of 2,922,234 $/year. When the BSC was set to 3,000,000, a TAC 2,921,233 $/year was obtained, the corresponding structure is presented in Figure 10b. Figure 10c displays the optimal HEN obtained by setting BSC to 2,500,000, whose TAC is 2,907,007 $/year within 8150.20 s and is 2899 $/year less as compared to that by Pavao et al. [39]. The selection of back substitution cycle has significant impact on the optimization process. When BSC is set too big, the individuals in long-term stagnation will not be adjusted in time; when BSC is set too small, the positive effect of AIS will be weakened.    The comparison of results is put in Table 5. Obviously, our optimal HEN is best, an inner utility (2556.8 kW) can also be seen at a split branch of C1. Furthermore, by observing the temperatures of C1 and HU, the match between C1 and HU is infeasible when C1 is heated to over 250 • C, so the HU at the end of C1 can only be 0 kW (C1 is directly be heated to the target temperature by the matches with hot streams) and over 5000 kW (the outlet temperature of C1 is lower than 250 • C).

Case Study 3
This is a large scale case including thirteen hot streams and seven cold streams, its data are given in Table 6. Escobar and Trierweiler [48] proposed an initialization strategy to generate better starting points, then they reported a result 1,461,006 $/year that was revised to 1,537,086 $/year by Pavao et al. [25]. Then the later authors provided a TAC 1,516,482 $/y with an improved PSO. By employing the FS strategy, Xiao et al. [43] gave a TAC 1,447,482 $/year. Applying the proposed method IRWCE and MSWS to the case, a HEN depicted in Figure 11 with the TAC 1,434,297 $/year was obtained, saving 13,185 $/year as compared to the result by Xiao et al. [43]. The computational time of this large-scale case is 28,517.52 s. Table 7 presents the comparison of the results. Obviously, our solution has less utility consumptions and higher precision.  Moreover, an inner hot utility with 825.3 kW is located at C7. Figure 12 gives the comparison of two structures with different utility placement, Figure 12a is the inner utility from Figure 11, Figure 12b is the structure where the utility is moved to the stream end. The two structures have the same heat exchange loads, but total areas in Figure 12a is 286.1 m 2 less than that in Figure 12b, indicating that the placement of inner utility can help save the areas. In effect, the placement of utility at the stream end reduces the heat exchange temperature difference of the match between HU and C7 and then causes the area to increase from 24.0 m 2 to 814.9 m 2 .  Moreover, an inner hot utility with 825.3 kW is located at C7. Figure 12 gives the comparison of two structures with different utility placement, Figure 12a is the inner utility from Figure 11, Figure 12b is the structure where the utility is moved to the stream end. The two structures have the same heat exchange loads, but total areas in Figure 12a is 286.1 m 2 less than that in Figure 12b, indicating that the placement of inner utility can help save the areas. In effect, the placement of utility at the stream end reduces the heat exchange temperature difference of the match between HU and C7 and then causes the area to increase from 24.0 m 2 to 814.9 m 2 . The parameters of IRWCE used for the three cases are given in Table 8. Furthermore, for all the three cases, the and presented in Equation (48) are set based on whether the corresponding heat exchanger load is zero or not. The values of for the The parameters of IRWCE used for the three cases are given in Table 8. Furthermore, for all the three cases, the S SRH and S SRC presented in Equation (48) are set based on whether the corresponding heat exchanger load is zero or not. The values of φ for the three cases are respectively given in Equations (49)-(51).

Conclusions
A novel stochastic method IRWCE consisting of RWCE, back substitution strategy, and elite optimization strategy is proposed to solve the HENS problems based on a modified SWS model considering inner utilities. The two novel strategies can maintain the positive effect of AIS and also solve the limitations of AIS, making the IRWCE satisfy the needs of relatively strong global and local search ability. The back substitution strategy is an effective solver of L1, which can deliver the historical optimums to the corresponding individuals in long-term stagnation condition to adjust their evolution direction in time. The elite optimization strategy can deal with L2 by performing a fine search for the selected elite solutions. It was also found that AIS might destroy the efficient structure and continuous variables, significantly influencing the optimization process. The optimal TACs of the three cases obtained by IRWCE were respectively 5,586,395 $/year (5,713,746 $/year when the fixed charge was set to 8000 $/year), 2,907,007 $/year, and 1,434,297 $/y, which saved 2219 $/year (1280 $/year), 2899 $/year, and 13,185 $/year as compared to the reported optimal ones. Some HEN structures with inner utilities were also obtained and significantly reduced the TACs, demonstrating the effectiveness of the MSWS model that enlarged the solution space. The investigative results demonstrated IRWCE and MSWS could achieve better savings in energy consumptions and design cost for HENS problems.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
A heat unit area (m 2 ) AMTD arithmetic mean temperature difference ( • C) BSC back substitution cycle C A area cost coefficient ($/m 2 ) C CU cost coefficient of cold utility ($/(kW y)) C F fixed charge of heat units C HU cost coefficient of hot utility ($/(kW y)) E 0-1 binary variable representing whether the heat unit exists or not h individual heat exchange coefficient of streams (kW/(m 2 K)) K overall heat exchange coefficient (kW/(m 2 K)) LMTD logarithmic mean temperature difference ( • C) N population size N C number of cold streams N H number of hot streams N S number of stages NE maximal heat units (inner cold utilities or inner hot utilities) NP number of elite individuals Q heat load (kW) Q HU,j hot utility at stream end, kW Q HU,i,j,k hot utility at stream branch, kW Q max initial maximal heat load (kW) Q min,b lower bound of heat unit loads for the basic individuals (kW) Q min,e lower bound of heat unit loads for the elite individuals (kW) Q n,best heat loads of historical optimum for individual n (kW) Q n,best heat loads of individual historical optimum after perturbation (kW) Q n,p,it initial load of the pth heat exchanger for individual n at the ztth iteration (kW) Q n,p,it+1 heat load of the pth heat exchanger after random walking for individual n at the ztth iteration (kW) Q n,p,it+1 effective heat load of the pth heat exchanger after random walking for individual n at the itth iteration (kW) QICU max initial maximal load of inner cold utility (kW) QICU n,l,it initial inner cold utility at the it iteration, kW QICU n,l,it+1 inner cold utility after walking at the itth iteration, kW QI HU max initial maximal load of inner hot utility (kW) QI HU n,m,it initial inner hot utility at the ztth iteration, kW QI HU n,m,it+1 inner hot utility after walking at the itth iteration, kW r random numbers uniformly distributed in (0, 1) S ICU maximal walking step size of inner cold utility (kW) S I HU maximal walking step size of inner hot utility (kW) S Q maximal walking step size of heat exchanger load (kW) S SRC maximal walking step size of split ratios of cold streams S SRH maximal walking step size of split ratios of hot streams SP split ratio SRC max initial maximal split ratio of cold streams SRC n,p,it initial split ratio of cold stream at the ztth iteration SRC n,p,it+1 split ratio of cold stream after walking at the ztth iteration SRC n,p,it+1 effective split ratio of cold stream after walking at the ztth iteration SRH max initial maximal split ratio of hot streams SRH n,p,it initial split ratio of hot stream at the ztth iteration SRH n,p,it+1 split ratio of hot stream after walking at the ztth iteration SRH n,p,it+1 effective split ratio of hot stream after walking at the ztth iteration T in