*2.1. Fundamentals*

The model for HEN optimization in this work is based on counterflow heat exchangers for the process-to-process heat exchangers as well as the utility heat exchangers. The heat load . *Q* in kW is calculated as the product of overall heat transfer coefficient *U* in kW/(m2K), the heat exchanger area *A* in m<sup>2</sup> and the logarithmic mean temperature difference (LMTD) Δ *Tm* in K:

$$
\dot{Q} = \mathcal{U} \cdot A \cdot \Delta T\_m \tag{1}
$$

The overall heat transfer coefficient is calculated via the individual heat transfer coefficients *hh* and *hc* in kW/(m2K) of the connected hot and cold process streams. The thermal resistance of the wall is neglected:

$$
\Delta U = \left(\frac{1}{h\_{\text{fl}}} + \frac{1}{h\_{\text{c}}}\right)^{-1} \tag{2}
$$

Due to the consideration of counterflow heat exchangers, the LMTD is used as temperature driving force for heat transfer:

$$
\Delta T\_m = \frac{\left(T\_h' - T\_c''\right) - \left(T\_h'' - T\_c'\right)}{\ln\left(\frac{T\_h' - T\_c''}{T\_h'' - T\_c'}\right)}\tag{3}
$$

The correlations for hot and cold outlet temperatures for the heat exchangers are as follows:

.

.

$$T\_h'' = T\_h' - \frac{\mathcal{Q}}{\dot{\mathcal{W}}\_h} \tag{4}$$

and:

$$T''\_{\mathcal{c}} = T'\_{\mathcal{c}} + \frac{\mathcal{Q}}{\dot{W}\_{\mathcal{c}}} \tag{5}$$

In order to represent different HEN solutions, a stage-wise superstructure proposed by Yee et al. [21] is used. As an example, a superstructure for two hot process streams *j* and three cold process streams *k* with stages *i* is shown in Figure 1 including the nomenclature used in the following equations.

**Figure 1.** Stage-wise superstructure of a heat exchanger network (HEN) with two hot and three cold process streams.

In each stage of the superstructure each hot stream is possibly connected with each cold stream. At the end of the streams a utility heat exchanger may be placed to ensure the desired target temperatures to ge<sup>t</sup> reached. Within the superstructure every heat exchanger has a specific index *ijk* depending on the positioning:

$$ijk = (i-1) \cdot N\_{\hbar} \cdot N\_{\hbar} + (j-1) \cdot N\_{\hbar} + k \tag{6}$$

Rearranging Equations (1)–(5) according to the outlet temperatures leads to a formulation which allows for the numerical calculation of the HEN temperatures:

$$
\begin{bmatrix} T''\_{h,ijk} \\ T''\_{c,ijk} \end{bmatrix} = \begin{bmatrix} \frac{(1-R\_{h,ijk})\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}}{1-R\_{h,ijk}\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}} & \frac{1-\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}}{1-R\_{h,ijk}\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}}\\ \frac{R\_{h,ijk}\left(1-\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}\right)}{1-R\_{h,ijk}\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}} & \frac{1-R\_{h,ijk}}{1-R\_{h,ijk}\varepsilon^{-NTl\_{ijk}^{'}(1-R\_{h,ijk})}} \end{bmatrix} \begin{bmatrix} T'\_{h,ijk} \\ T'\_{c,ijk} \end{bmatrix} \tag{7}
$$

with *Rh*,*ijk* being:

$$R\_{h,ijk} = \frac{\mathcal{W}\_{h,ijk}}{\dot{\mathcal{W}}\_{c,ijk}}\tag{8}$$

and:

$$NTU\_{ijk} = \frac{(UIA)\_{ijk}}{\dot{W}\_{h,ijk}}\tag{9}$$

The corresponding non-isothermal mixing temperature for each process stream after each stage is calculated based on an energy balance [22]. The maximum number of stages to represent a HEN during optimization is set to *Ns* = *max*{*Nh*, *Nc*} [22].

.

#### *2.2. Objective Function*

The objective function is structured in such way that a flexible cost functions framework can ge<sup>t</sup> incorporated in order to consider different implementation related factors like piping and pumping. A different cost function for every coupling of possible hot and cold process stream matches can be defined. In addition, different cost functions can ge<sup>t</sup> added for each utility heat exchanger. If necessary they can be further refined for every stage in the superstructure representation of the HEN. The objective function is based on the proposed structure by Rathjens and Fieg [1] and was further developed for the present work:

$$\begin{aligned} \min \left\{ \sum\_{n=1}^{N\_h + N\_c} \left[ \mathbf{C}\_{\text{CLI}} \cdot \max \left\{ \dot{W}\_n \left( T\_n^{\text{v}} - T\_{\text{out},n}^+ \right), 0 \right\} + \mathbf{C}\_{\text{HLI}} \cdot \max \left\{ \dot{W}\_n \left( T\_{\text{out},n}^- - T\_n^{\text{v}} \right), 0 \right\} \right] \\ \quad + \sum\_{i=1}^{N\_e} \sum\_{j=1}^{N\_h} \sum\_{k=1}^{N\_c} z\_{ijk} X\_{ijk} + \sum\_{n=1}^{N\_h + N\_c} \left[ z\_{\text{CLI},n} X\_{\text{CLI},n} + z\_{\text{HLI},n} X\_{\text{HLI},n} \right] \right\} \end{aligned} \tag{10}$$

Hot utility and cold utility are used to reach the desired target outlet temperatures *Tout*,*<sup>n</sup>* and cause costs based on their individual supply costs *CHU* and *CCU*. During optimization hot and cold utility can potentially be applied to each process stream depending on the maximum positive or negative deviation from the specified bounds of target temperatures. The respective utility which remains unused will make no contribution to the total annual costs (TAC). The binary variables *z* state the existence of a heat exchanger within the superstructure. The match-based costs *X* are calculated based on different cost function formulations with a wide range of possible correlations that can ge<sup>t</sup> implemented. The common cases of using a power function with the coefficients *a*0, *a*1 and *a*2 as well as a general polynomial representation up to 4th degree (coefficients *b*0, *b*1, *b*2, *b*3 and *b*4) are included. More specialized dependencies are covered using exponential (coefficients *c*0 and *c*1) or logarithmic expressions (coefficients *d*0 and *d*1):

$$\begin{array}{c} a\_{0,ijk} + a\_{1,ijk}\mathbf{z}^{\text{R2},ijk} + \\ X\_{ijk} = \begin{array}{c} b\_{0,ijk} + b\_{1,ijk}\mathbf{z}^{\text{R}} + b\_{2,ijk}\mathbf{z}^{\text{R}} + b\_{3,ijk}\mathbf{z}^{\text{S}} + b\_{4,ijk}\mathbf{z}^{\text{S}} \\ c\_{0,ijk} + e^{c\_{1,ijk}}\mathbf{z} + \\ d\_{0,ijk}\ln(\mathbf{z}^{\text{S}}) + d\_{1,ijk} \end{array} \end{array} \tag{11}$$

The different cost function formulations given in Equation (11) are chosen because they were able to represent the cost correlations for practical implementation that were faced when working together with our industrial partners. In order to be used in the objective function, the dependencies given in Equation (11) have to ge<sup>t</sup> expressed through the individual heat exchanger areas *A*, the temperature levels *T* and heat capacity flow rates . *W* of the involved streams and the heat load . *Q*:

$$\,\_5\exists\left\{\mathcal{A}\_{\mathrm{i}\bar{\mathrm{j}}\mathrm{k}\nu}\dot{\mathcal{W}}\_{\mathrm{h},\bar{\mathrm{j}}\mathrm{j}\nu}\dot{\mathcal{W}}\_{\mathrm{c},\bar{\mathrm{j}}\mathrm{j}\nu}\min\left\{\dot{\mathcal{W}}\_{\mathrm{h},\bar{\mathrm{j}}\mathrm{j}\nu}\dot{\mathcal{W}}\_{\mathrm{c},\bar{\mathrm{j}}\mathrm{j}\bar{\mathrm{k}}}\right\},\max\left\{\dot{\mathcal{W}}\_{\mathrm{h},\bar{\mathrm{j}}\mathrm{j}\nu}\dot{\mathcal{W}}\_{\mathrm{c},\bar{\mathrm{j}}\mathrm{j}\bar{\mathrm{k}}}\right\},\boldsymbol{T}\_{\mathrm{h},\bar{\mathrm{j}}\mathrm{j}\nu'}\_{\mathrm{c},\bar{\mathrm{j}}\mathrm{j}\bar{\mathrm{k}}'}\boldsymbol{T}\_{\mathrm{c},\bar{\mathrm{j}}\mathrm{j}\bar{\mathrm{k}}'}^{\prime}\dot{\mathcal{Q}}\_{\mathrm{c},\bar{\mathrm{j}}\mathrm{j}\bar{\mathrm{k}}'}\right\}\tag{12}$$

The expression of cost dependencies through distinct variables was chosen due to ensuring the universal applicability of the shown approach without extending the programming work for each different problem. Furthermore these variables can be used with minimal computational overhead. Because *ξ* can be chosen out of an amount of alternatives given in Equation (12) the combinatorial possibilities of representing different cost dependencies are immense. The costs of utility heat exchangers *XHU*,*<sup>n</sup>* and *XCU*,*<sup>n</sup>* are structured equivalently to *Xijk*.

#### *2.3. Genetic Algorithm*

The genetic algorithm uses the common genetic operations like selection, crossover and mutation. As optimization variables the parameters (*UA*)*ijk*, . *Wh*,*ijk* and . *Wc*,*ijk* are used like proposed by Fieg et al. [22]. The selection is based on the fitness value *F*. The fitness value represents the quality of an individual with respect to the objective function. It is based on the relative relation towards the average costs among all individuals and the minimum TAC of all individuals. The calculation of the fitness value is given in the following equation:

$$F = \frac{\mathbb{C}\_{TAC}^{-1} + \mathbb{C}\_{TAC, \text{min}}^{-1} - 2\mathbb{C}\_{TAC, \text{avg}}^{-1}}{\mathbb{C}\_{TAC, \text{min}}^{-1} - \mathbb{C}\_{TAC, \text{avg}}^{-1}} \tag{13}$$

The selection is carried out according to roulette wheel selection. Crossover is split in parameter crossover and structure crossover. Mutation is carried out on the variable parameters stated above. For details about the explicit formulations please refer to Fieg et al. [22]. The probability for the application of a crossover operation is 89% in this work and the probability for parameter crossover was chosen as 23% according to Brandt [23]. The probability for the general mutation is 1%, for parameter mutation it is 50% with a respective gene mutation probability of 1% [23].

Because of the optimization of the parameter (*UA*)*ijk* using the explicit temperature solution depicted in Equation (7) instead of the heat loads of the heat exchangers, every solution is thermodynamically feasible. Otherwise occurring temperature or heat load constraints can be omitted. The binary variable handling and corresponding constraints considering the continuous variables are adopted from Luo et al. [20]. During optimization the strategy of excessive use of utilities is used [24]. If the specified hot utility temperature is higher than any target temperature of all process streams and cold utility temperature is lower than any cold target temperature of all process streams, each network generated throughout optimization is feasible. As a result no outlet temperature constraints are considered here. Because . *Wh*,*ijk* and . *Wc*,*ijk* are also used as optimization variables, these values have to ge<sup>t</sup> constrained to be within a feasible region to hold the constraints given in Equations (14) and (15):

$$\sum\_{k=1}^{N\_{\mathcal{E}}} \dot{\mathcal{W}}\_{h,ijk} = \dot{\mathcal{W}}\_{h,j} \tag{14}$$

$$\sum\_{j=1}^{N\_b} \dot{\mathcal{W}}\_{\varepsilon,ijk} = \dot{\mathcal{W}}\_{\varepsilon,k} \tag{15}$$

During optimization . *W*<sup>∗</sup> *h*,*ijk* and . *W*<sup>∗</sup> *<sup>c</sup>*,*ijk* can occur that violate the above mentioned constraints. In order to keep these constraints, the normalization strategy from Fieg et al. is applied [22], which is shown by the following equations:

$$
\dot{\mathcal{W}}\_{h,ijk} = \frac{\dot{\mathcal{W}}\_{h,\dot{j}}}{\sum\_{k=1}^{N\_c} \dot{\mathcal{W}}\_{h,ijk}^\*} \dot{\mathcal{W}}\_{h,ijk}^\* \tag{16}
$$

$$
\dot{\mathcal{W}}\_{c,ijk} = \frac{\mathcal{W}\_{c,k}}{\sum\_{j=1}^{N\_h} \dot{\mathcal{W}}\_{c,ijk}^\*} \dot{\mathcal{W}}\_{c,ijk}^\* \tag{17}
$$

Furthermore, the structural control strategy proposed by Luo et al. [20] was adapted to ensure heterogeneity in the population and avoid local optimal solutions.

.

Structurally forbidden matches are removed automatically from a possible occurrence in the superstructure by setting the respective (*UA*)*ijk* to zero. A simple assignment of high costs for the specific match is not efficient against the background of fast convergence. This is particularly true for the second example shown below.

#### *2.4. Local Optimization*

Due to the strong nonlinearity of the stated objective function, an algorithm relying only on genetic operations would take too long to find an acceptable solution. Therefore, deterministic methods are used for local optimization. In order to generate promising HEN structures in the initialization step the approach of enhanced vertical heat transfer proposed by Stegner et al. [25] is used. The method proposed by Stegner et al. [25] adopts ideas known from the conventional vertical heat transfer concepts from traditional pinch approach but a new form of graphic depiction was implemented. Loop breaking as one of the heuristic approaches known from pinch technology is considered as well. This approach was successfully combined with stochastic optimization by Brandt et al. [26].

Furthermore, Newton's method is utilized for local parameter optimization [22]. Equation (18) shows the nomenclature in Newton's method for the variables (*UA*) ∗ *ijk*, . *W*<sup>∗</sup> *h*,*ijk* and . *W*<sup>∗</sup> *<sup>c</sup>*,*ijk* to ge<sup>t</sup> the optimized expression:

$$\left\{ (\mathcal{U}A)\_{ijk\prime} \dot{\mathcal{W}}\_{h,ijk\prime} \dot{\mathcal{W}}\_{c,jik} \right\} = \left\{ (\mathcal{U}A)^{\*}\_{ijk\prime} \dot{\mathcal{W}}^{\*}\_{h,ijk\prime} \dot{\mathcal{W}}^{\*}\_{c,jik} \right\} - \frac{\frac{\partial \mathcal{C}\_{\mathcal{I}AC}}{\partial \left\{ (\mathcal{U}A)^{\*}\_{ijk\prime} \dot{\mathcal{W}}^{\*}\_{h,ijk\prime} \dot{\mathcal{W}}^{\*}\_{c,jik} \right\}}{\frac{\partial^{2} \mathcal{C}\_{\mathcal{I}AC}}{\partial \left\{ (\mathcal{U}A)^{\*}\_{ijk\prime} \dot{\mathcal{W}}^{\*}\_{h,ijk\prime} \dot{\mathcal{W}}^{\*}\_{c,jik} \right\}^{2}}} \tag{18}$$

The partial derivatives needed in Equation (18) are calculated numerically using a step size of 10−4.

#### **3. Examples and Results**

In order to show the application of the presented approach two examples have been chosen. Example 1 is a mid-sized optimization problem for direct heat integration. Example 2 is an optimization problem for indirect heat integration using a heat recovery loop (HRL).

#### *3.1. Example 1*

Example 1 was first stated by Pho and Lapidus [27] and was named 10SP1. It consists of five hot and five cold streams. In the first part of the analysis in our work we focus on the cost functions already used in other publications to achieve comparable results (see Section 3.1.1). In the second part we focus on additional costs caused by the consideration of piping (see Section 3.1.2).

#### 3.1.1. Example 1a

The 10SP1 problem is well studied and has been used as an optimization case several times so far with decreasing costs throughout the years by many researchers (e.g., Nishida et al. [28] or Flower and Linnhoff [29].

During the last two decades the 10SP1 problem has still attracted a considerable amount of attention. The associated problem data is given in Table 1. Two different formulations for the cost calculation of the annual costs of heat exchangers for the 10SP1 problem are given in the literature. The first cost formulation is given in Equation (19):

$$X\_{i\bar{j}k/H\mathcal{U}\_n/\mathcal{C}\mathcal{U}\_n} = 140 \cdot A\_{i\bar{j}k/H\mathcal{U}\_n/\mathcal{C}\mathcal{U}\_n}^{0.6} \tag{19}$$

The second cost formulation is given in Equation (20):

$$X\_{\rm ijk/HLl\_n/CLl\_n} = 145.63 \cdot A\_{\rm ijk/HLl\_n/CLl\_n}^{0.6} \tag{20}$$

A summary of the achieved TAC is given in Table 2. Several different approaches were used to tackle the 10SP1 problem. Lewin et al. [24] used a two-level approach consisting of a genetic algorithm and a lower level transformation into a linear parametric optimization problem. Lewin [30] also used an alternative nonlinear programming approach for the lower level optimization. Lin and Miller [31] used a tabu search algorithm to tackle the 10SP1 problem. Pariyani et al. [32] used a randomized algorithm for problems with stream splitting and a modified version of a previous work from Chakraborty and Ghosh [33] for designing networks without stream splitting. Yerramsetty and Murty [34] used a differential evolution algorithm and Peng and Cui [35] utilized a two-level simulated annealing algorithm. The latest work of Aguitoni et al. [36] found the best solution for the cost formulation 1 with a combination of a genetic algorithm and differential evolution.


**Table 1.** Problem data for example 1a.


**Table 2.** Results comparison for example 1a.

1Solution without stream splits. 2 Solution with stream splits. 3 Revised by Pariyani et al. [32].

The optimization procedure with the chosen approach incorporates 375 optimization variables considering five hot streams, five cold streams and five stages for the superstructure. Each heat exchanger has the same area-related cost function.

In our work cost formulation 1 was used for the optimization. The minimum TAC found was 42,963 \$/yr. The solution is depicted in Figure 2 and shows a configuration of two stream splits as well as the use of two utility heat exchangers.

**Figure 2.** Grid-diagram representation of the optimal HEN configuration of example 1a (TAC: 42,963 \$/yr).

The cost formulation for the optimization is favoring maximum energy recovery (MER) networks. Therefore, the obtained network recovers the maximum amount of energy and only cold utility is used. The TAC obtained for example 1a is significantly lower than in the publications cited in relation to the progress in TAC made over the years. The HEN shown in Figure 2 was also evaluated using cost formulation 2 and the resulting costs are 43,321 \$/a, which is the lowest TAC found for formulation 2 so far. The results show that the chosen approach is capable of generating competitive results compared to other approaches presented in literature. Therefore, the consideration of even more complex problem formulations should be manageable, which is shown in the following parts.

#### 3.1.2. Example 1b

Example 1b is built on the problem definition from example 1a. In addition, piping costs are considered. The required coordinates for the corresponding streams are taken from Pouransari and Maréchal [18]. They added arbitrary Cartesian coordinates to the 10SP1 problem stated above. For the hot utility coordinates the values from Rathjens and Fieg [1] are used. The additional data is given in Table 3.


**Table 3.** Additional Cartesian coordinates for the 10SP1 problem.

The required inner pipe diameter *Di* in m for each possible connecting pipe between a heat source and heat sink is calculated with the correlation for optimal diameter given by Peters et al. [37] for turbulent flow and *Di* ≥ 0.0254 m:

$$D\_i = 0.363 \cdot \dot{V}^{0.45} \text{ g}^{0.13} \tag{21}$$

with the volumetric flow rate *V* in m3/s, the fluid density of 983 kg/m<sup>3</sup> and a specific heat capacity of 4.18 kJ/(kgK).

.

Peters et al. [37] give a diagram for the correlation between pipe diameter and the related costs per meter of pipe. Due to manually selecting and reading out values from this diagram there might be a slight deviation from the book values which cannot be further specified. Considering stainless-steel welded pipe of type 304 [37] results in the following cost correlations with respect to the stream heat capacity flow rates and heat exchanger areas:

$$\begin{aligned} X\_{\bar{i}\bar{j}k} &= -140 \cdot A\_{\bar{i}\bar{j}\bar{k}}^{0,6} + L\_{\bar{j}\bar{k}} \cdot \\ &\quad \left( -0.174 \cdot \left( \min \left\{ \dot{W}\_{h,\bar{j}j\bar{k}\prime} \dot{W}\_{c,\bar{j}\bar{k}} \right\} \right)^2 + 17.773 \cdot \left( \min \left\{ \dot{W}\_{h,\bar{j}j\bar{k}\prime} \dot{W}\_{c,\bar{j}\bar{k}\prime} \right\} \right) + 27.426 \right) \end{aligned} \tag{22}$$

and:

$$X\_{\rm CLI\_n} = 140 \cdot A\_{\rm CLI\_n}^{0.6} + L\_{\rm CLI\_n,j} \left( -0.174 \cdot \dot{W}\_{\rm CLI\_n}^2 + 17.773 \cdot W\_{\rm CLI\_n} + 27.426 \right) \tag{23}$$

with the distance *L* being:

$$L\_{jk\text{/CL}\_{n\text{/}}} = 2 \cdot \left( \left| \mathbf{x}\_{j} - \mathbf{x}\_{k/\text{CL}\_{n}} \right| + \left| y\_{j} - y\_{k/\text{CL}\_{n}} \right| + \left| z\_{j} - z\_{k/\text{CL}\_{n}} \right| \right) \tag{24}$$

The number of optimization variables is 375 as in example 1a. However, the number of potentially different cost functions is 25 (each match between hot and cold stream). The number of potentially different cost functions for the hot and cold utility usage is 10 for each utility. The real number of different cost functions is lower due to the fact that some distances between multiple hot and cold streams as well as utilities are equal.

The relatively short distances between the heat sources and heat sinks are compensated by the annualization of the piping costs over a one-year period. As a result the portion of piping costs becomes around 35% of TAC.

The solution with the lowest TAC found in this work is shown in Figure 3a. It has a TAC of 68,476 \$/yr. A structurally different solution which was found during optimization with a less complex structure but the same utility configuration is shown in Figure 3b.

**Figure 3.** Grid-diagram representations of: (**a**) the optimal HEN configuration of example 1b (TAC: 68,476 \$/yr); and (**b**) an alternative HEN configuration of example 1b (TAC: 68,843 \$/yr).

As opposed to the solution of example 1a only two stages within the superstructure are taken by heat exchangers. Furthermore, the number of utility heat exchangers increased to three compared to the case not considering piping. The results still represent MER networks. Figures 4 and 5 show the relative orientation of the heat sources and heat sinks towards each other.

**Figure 4.** (**a**) 3D layout representation of the optimal HEN configuration of example 1b, the chosen line thickness is proportional to the corresponding heat exchanger heat load; and (**b**) 3D layout representation of the alternative HEN configuration of example 1b, the chosen line thickness is proportional to the corresponding heat exchanger heat load.

The configuration containing three utilities was dominant throughout the optimizations. The consideration of piping costs and thus the changed objective function had a strong influence on the heat load distribution within the HENs. The comparison between the optimal results of example 1b and example 1a shows about 14.2% (4016 \$/yr) less piping costs, which is counterbalanced by an increase in installed heat exchanger area of about 16.2% (39.4 m2) and thus an increase in capital costs for the heat exchangers of 14.7% (1307 \$/yr).

**Figure 5.** (**a**) 2D layout representation of the optimal HEN configuration of example 1b; (**b**) 2D layout representation of the alternative HEN configuration of example 1b.

The favorable design of using cold utility for the streams H3 and H5 is still preferred when considering piping costs during optimization due to the coordinates chosen by Pouransari and Maréchal [18]. Furthermore the optimization is dominated by the close arrangemen<sup>t</sup> of the streams H3 and C4 as well as the streams H4, H5 and C3. Despite that, the solutions found considering piping costs already during the optimization are structurally different. Rathjens and Fieg [1] used a different cost functions formulation and were able to show significantly more local clustering of heat exchange similar to Figure 5b.

#### *3.2. Example 2*

The second example is a heat integration case with indirect heat integration utilizing a HRL with an intermediate fluid. It was taken from Chang et al. [13]. It is a case study of a heat integration project in the southern part of China. The example comprises two different plants with a distance *L* of 1000 m. One plant is the heat source plant and the other one the heat sink plant with seven process streams each. The problem data is given in Table 4.


**Table 4.** Problem data for example 2.

*CHU* = 20 \$/(kWyr), *CCU* = 8 \$/(kWyr) [38]; *CHX* = 11, 000 + 150·*AHX* [39].

In addition to the utility costs and the capital costs for the heat exchangers, the costs for pumping and piping are considered. The detailed equations are given in the following part.

For the calculation of the inner pipe diameter *Di* for the HRL, Equation (21) is used. The fluid density of 960 kg/m<sup>3</sup> and a specific heat capacity flow rate of 4.2 kJ/(kgK) is assumed.

The outer diameter *Dout* in m, the specific pipe weight *wtpipe* in kg/m and the resulting pipe capital costs *Pcul* in \$/m for schedule 80 steel pipes are calculated according to Stijepovic and Linke [39]:

$$D\_{\rm out} = 1.101 \cdot D\_{\rm in} + 0.006349 \tag{25}$$

$$wt\_{pipe} = 1330 \cdot D\_{in}^2 + 75.18 \cdot D\_{in} + 0.9268 \tag{26}$$

$$Pcul = 0.82 \cdot wt\_{pipc} + 185 D\_{out}^{0.48} + 6.8 + 295 \cdot D\_{out} \tag{27}$$

The resulting capital costs for piping in \$ are:

$$C\_{pipe} = 2 \cdot L \cdot P\_{cal} \tag{28}$$

For the calculation of the pumping capital costs the fanning friction factor *f* , the Reynolds number *Re* and the fluid velocity *u* in m/s are used to estimate the pressure drop Δ*p* in Pa and the resulting costs [13,39]. The fluid velocity is calculated according to Equation (29) using the fluids mass flow rate . *m* in kg/s of the intermediate fluid:

$$
u = \frac{4 \cdot m}{\varrho \cdot \pi \cdot D\_{in}^2} \tag{29}$$

The Reynolds number is calculated assuming a dynamic viscosity *μ* of 0.0002834 Pa s:

$$\mathcal{Re} = \frac{\varrho \cdot \mu \cdot D\_{in}}{\mu} \tag{30}$$

The resulting pressure drop Δ*p* is given by:

$$
\Delta p = 4 \cdot f \cdot \frac{L \cdot \varrho \cdot u^2}{2 \cdot D\_{in}} \tag{31}
$$

with the fanning friction factor of:

$$f = \frac{0.046}{Re^{0.2}}\tag{32}$$

The cumulated capital costs for the two pumps *Cpump* in \$ for the HRL are calculated according to Equation (33), utilizing the corresponding parameters for centrifugal pumps given in Jabbari et al. [40]:

$$\mathcal{C}\_{pump} = 2 \cdot \left( 8600 + 7310 \cdot \left( \frac{\dot{m} \cdot \Delta p}{\ell} \right)^{0.2} \right) \tag{33}$$

The operating costs for pumping *Cpump*,*op* in \$/yr are composed of the pump efficiency *ηpump* of 0.7, the price for electricity *Cel* of 0.1 \$/(kWh) [13], an operation duration *t* of 8000 h/yr and the parameters calculated beforehand [13]:

$$\mathbb{C}\_{pump,op} = 2 \cdot \mathbb{C}\_{cl} \cdot t \cdot \frac{\dot{m} \cdot \Delta p}{\varrho \cdot \eta\_{pump}} \tag{34}$$

The annualization of the capital costs for the heat exchangers, piping and the required pumps is done using an annualization factor *AF* in 1/yr. The annualization factor is calculated according to Chang et al. [13] with an operation time *n* of five years and a fractional interest rate *I* of 10% per year:

$$AF = \frac{I \cdot \left(1 + I\right)^n}{\left(1 + I\right)^n - 1} \tag{35}$$

In order to use the superstructure of the utilized optimization model in this work, two pseudo-streams are added to the problem definition to emulate the HRL. The respective heat transfer coefficients for the two streams are 1 kW/(m2K) each [13,14]. The resulting superstructure has eight stages as well as eight hot streams and eight cold streams. In order to represent the cost correlations described beforehand, they are transferred to a dependency in heat capacity flow rate to ge<sup>t</sup> integrated in the cost functions framework proposed in this work. Therefore, the costs for piping and pumping are added up and calculated for different mass flow rates. The resulting costs are plotted against heat capacity flow rates calculated based on the respective mass flow rates. A curve fitting is carried out afterwards. The deviation from calculating the exact costs based on the equations given above is below 0.0012% considering rounded values for the piping and pumping costs. The resulting correlations are shown in Equation (36):

$$X\_{ijk} = \begin{cases} 2910.77 + 39.57 \cdot A\_{ijk}^{1} & 2 \le i \le 7, 1 \le j \le 7, k = 1\\ (109, 107.62 + 300.52 \cdot \dot{W}\_{h,j} - 0.1019 \cdot \dot{W}\_{h,j} & \text{ijk} = \{57, 505\} \\ +3.7486 \cdot 10^{-5} \cdot \dot{W}\_{h,j}^{3} - 6.1293 \cdot 10^{-9} \cdot \dot{W}\_{h,j}^{4})/2, & \text{(36)} \\ 2910.77 + 39.57 \cdot A\_{ijk}^{1} & 2 \le i \le 7, j = 8, 2 \le k \le 8 \\ \text{forbidden}, & \text{otherwise} \end{cases} \tag{36}$$

The utility heat exchangers for the two plants are considered as already installed like it was done by Chang et al. [13]. Therefore, *XHU*,*<sup>n</sup>* and *XCU*,*<sup>n</sup>* are zero for this example. The first and last stages are used to include the costs for pumping and piping. The areas of the heat exchangers to emulate the HRL (*ijk* = {57, 505}) are chosen to be near infinite by the algorithm to avoid utility usage, thus create a self-adjusting temperature level and balance the heat transferred from the heat source plant to the heat sink plant. In order to prohibit direct heat integration between the two plants, all direct matches are forbidden. Altogether the applied superstructure yields 1536 optimization variables (compared to 1029 for direct heat integration without the pseudo streams). Due to the large number of forbidden matches (426) in relation to the number of possible matches (512), the number of optimization variables considered is reduced equivalently. The whole problem is described with only three different cost functions (see Equation (36)).

Iteration over the mass flow rate yields the optimal configuration for the HRL for example 2 which is shown in Figure 6. The superstructure representation of the solution in Figure 6 is shown in Figure A1 in Appendix A.

The TAC achieved for the solution shown in Figure 6 is 1.54 M\$/yr. The associated annualized capital costs for the heat exchangers are 0.574 M\$/yr (37.2% of TAC) and the utility consumption causes annual costs of 0.603 M\$/yr (39.1% of TAC). The overall piping costs are 0.326 M\$/yr (21.1% of TAC) and pumping costs are 0.038 M\$/yr (2.5% of TAC). The solution obtained in this work is therefore 4.4% cheaper than the solution reported by Chang et al. [13]. These differences are explained by the non-isothermal mixing model used and the relaxed temperature constraints in this work.

**Figure 6.** Plant representation of the optimal HEN configuration found of example 2 including a HRL (TAC: 1.54 M\$/yr).
