Enhanced superstructure optimization for heat exchanger network synthesis using deterministic approach

Heat Exchanger Network (HEN) synthesis is primarily formulated as a mixed integer non-linear programming (MINLP) problem based on method of stage-wise superstructure (SWS). Approaches to obtain an optimal HEN configuration can adopt deterministic algorithms. But, as a large scale problem, it is difficult to solve due to the complexities arisen from nonlinearities of SWS. To overcome the model nonlinearities, stochastic algorithms and meta-heuristic approaches have been proposed in the literature to tackle the problem. However, it reaches a near-optimal HEN configuration from a series of stochastic solutions, which are obtained by the execution of many computational procedures with extensive time resource, and a global optimum will not be guaranteed from only randomly generated results. In this paper, an enhanced SWS approach is proposed, in which new temperature and heat duty constraints are updated to reduce the redundant combinations and avoid conflicted calculation of non-isothermal mixing energy balance. The present model is also extended to allow flexible stream splitting for practical applications. Then, a deterministic-based global solver (GAMS/BARON) is applied in solving three case studies. The results show that the proposed approach can provide cost-efficient HEN solutions with lower TAC than using existing stochastic and deterministic algorithms.

Introduction

Nowadays, the global energy usage is still largely based on fossil fuels, including coal, oil, and natural gas, and shows an increase by 4.6% in 2021 compared with that in 2020. In the industrial sector, the chemical industries are particularly energy intensive, accounting for about 30% of the total energy usage and associated CO2 emissions. Heat exchanger networks are essential for the chemical process industries, since they can reduce energy consumption and pollutant emissions by heat integration of process hot and cold streams. However, heat exchangers also require significant capital investment at the same time. Industrial processes demand not only an efficient energy recovery system, but also a better economic performance, which leads to widespread studies on HEN synthesis.

The approaches used in HEN synthesis can be divided into pinch methods and automated methods.

Pinch analysis is based on sequential thermodynamic analysis, and was proposed by Linnhoff and Townsend (1982), Linnhoff and Hindmarsh (1983). This method defines the minimum process approach temperatures as a pinch point to implement maximum heat recovery, and has been widely employed to provide a cost-efficient solution, following the sequences: (1) Set energy target; (2) Minimizing heat exchanger fixed cost; (3) Minimizing exchanger area cost. However, it requires experienced designers and such sequential procedures may miss promising solutions, since simultaneous optimization cannot be achieved considering trade-off between capital and energy costs.

Mathematical programming has been developed to synthesize HEN automatedly. It commonly formulates the design problem through a superstructure, in which exchanger placements with or without split, heat exchanger inlet and outlet temperatures, utilities and heat duties can be optimized simultaneously. One of the widely used superstructures was proposed by Yee and Grossmann (1990), defined as stage-wise superstructure. Based on the use of SWS, the HEN synthesis is formulated as an MINLP problem, targeting the minimum total cost. Due to non-linearity, non-convexity and raised problem sizes for large-scale HEN optimization problems using SWS, it is challenging to find an optimal HEN solution with an acceptable computational time. Many approaches have been reported to solve the HEN synthesis problem using different algorithms.

Stochastic algorithms treat the HEN design with randomization, not limited by model non-convexity and discontinuity. Some papers have applied these methods in solving HEN synthesis, such as Genetic Algorithms (GA) (Holland, 1992), Simulated Annealing (SA) (Kirkpatrick et al., 1983), Particle Swarm Optimization (PSO) (Kennedy and Eberhart, 1995), Differential Evolution (DE) (Storn and Price, 1997) and Ant Colony Optimization (ACO) (Dorigo and Caro, 1999). Subsequently, hybrid-approaches have been widely studied with developed stochastic algorithms. As the basis of the two-level approach (Lewin, 1998), study by Huo et al. (2012) applied GA at the upper level to find optimal network topologies and used PSO to optimize the continuous variables simultaneously. The PSO showed good performance when employed at the lower-level optimization. The approaches by Pavão et al. (2016), Pavão et al. (2017a) respectively applied GA-PSO and SA-PSO algorithms, and they introduced parallel processing techniques considering phase change. A two-level SA method was also proposed by Peng and Cui (2015) without considering stream splitting. To improve the performance of approaches for solving large-scale problems, Pavão et al. (2017b), Pavão et al. (2018b), Pavão et al. (2018a) employed the SA coupled with a Rocket Fireworks Optimization (RFO) meta-heuristic method to solve large-scale problems, including additional options of multiple utilities, cross flow and integration with pump. Rathjens and Fieg (2020) studied to embed GA in structure identification and change of reference system (SIR) as tracking large-scale problems. Such research was enhanced by using structure identification to reduce the optimization variables in their work. GA was further incorporated with the modified quasi-linear programming model in the study by Feyli et al. (2022), in which GA and modified quasi-linear programming handle HEN structure discrete and continues variables separately.

Nevertheless, for adopting stochastic algorithms, extensive computational efforts are commonly required to determine near-optimal solutions. Deterministic approaches have been widely employed to solve HEN synthesis problems in early stage, such as Branch and Bound (BB), Outer Approximation (OA) (Zamora and Grossmann, 1997), SMIN- αBB and GMIN-αBB algorithms (Adjiman et al., 1997). Björk and Westerlund (2002) solved the HEN synthesis problem with non-isothermal mixing assumption. Their work adopted the global optimization approach through introducing variable inverse transformation function. Bergamini et al. (2008) applied the OA algorithm in their method with piece-wise linear under-estimators for nonconvex constraints. Hasan et al. (2010) proposed a new SWS with substages to consider phase change by adding nonlinear property correlations in their model.

In recent years, the development of a computational platform has brought applications of using deterministic approach to solve industrial-scale HEN synthesis problems. Studies by Huang et al. (2012), Huang and Karimi (2013) reported a modified superstructure to consider non-isothermal mixing and solved by commercial solver BARON/GAMS. Faria et al. (2015) developed the Bound Contraction method to solve the HEN synthesis problem, in which variables are partitioned to construct a linear or convex lower bound. Kim et al. (2017) extended the work by Faria et al. (2015) to present a deterministic-based global optimization. In their paper, a lower bound model and an upper bound model were solved first by BARON/GAMS or ANTIGONE/GAMS, then lifting partitioning was added to find an optimal solution. Chang et al. (2020) proposed a new global optimization approach by establishing the concept of minimal structure (MSTR) to solve iso-thermal mixing SWS. By adopting CPLEX/GAMS solver, their system can generate a cost-efficient solution. Nevertheless, their solution qualities were limited by the assumption of isothermal mixing, and demanded extensive computational time.

In this work, based on the SWS developed by Huang et al. (2012), an enhanced non-isothermal mixing superstructure is proposed. Additional variables and constraints are defined to improve the model performance and achieve the flexible splitting requirement. A global solver, BARON/GAMS is employed as the deterministic approach to solve the problem toward minimizing the total annual cost of HEN. Parallel computation is introduced to intensify the ability for solving large-scale mixed integer trees. The proposed optimization route presents efficient procedures, in which appropriate heat transfer stages used in SWS are explored.

Mathematical model

The mathematical formulation is based on the well-known stage-wise superstructure (SWS) proposed by Yee and Grossmann (1990), and Figure 1 presents the basic superstructure that includes two hot and two cold streams in two stages. Stages are introduced, in which all possible matches between hot streams and cold streams are optimized simultaneously. Hot and cold utilities are given to satisfy stream target temperatures if heating and cooling duties provided by the heat recovery system in HEN are not enough, whose locations are optimized at the end of process streams.

FIGURE 1

Figure 1. Stage-wise superstructure in two stages: k = 1 and k = K.

Huang et al. (2012) indicated that the iso-thermal mixing for sub-streams assumed by Yee and Grossmann (1990) can overestimate the total area as stage outlet temperatures are equal to all exchanger outlet temperatures in all splits, thus a non-isothermal mixing (NIM) SWS was presented by updating non-isothermal energy balance. However, their SWS includes additional stages (stage k = 0 and k = K + 1), which leads to an increase of problem sizes. The work by Pavão et al. (2017b) offered a simple NIM-SWS with some modifications. But the feasibility constraints applied in Pavão et al. (2017b) are not rigorous, which causes an issue that, the sum of split fractions of a stream in a single stage is always equal to 1, even if this stream bypasses the stage without any existed match. This leads heat exchanger energy balance to be invalid for exceptional cases when deterministic approaches are employed.

In this work, the formulation removes the iso-thermal mixing assumed in the original SWS method. The actual outlet temperature of hot split streams and cold split streams are calculated through heat exchanger energy balance by introducing additional continues variables. Several logical constraints and continues variables proposed by studies from Huang et al. (2012), Huang and Karimi (2013) are introduced to allow a stream to fully bypass a single stage and eliminate redundant combinations in SWS. The new continues variables are further used to control the total branch number for an individual stream. Furthermore, additional convex and concave constraints are adopted to improve the solution quality of non-isothermal mixing model.

Basically, the set of hot process streams (i; i ∈ HP), cold process streams (j; j ∈ CP), and stages (k; k ∈ ST) are defined in the present NIM-SWS, as shown in Figure 1. The initial number of stages is assumed as ST = max[i, j]. The number of stages can be adjusted according to the problem size, particularly for large-scale problems. This is because large number of stages lead to a great number of binary variables and combinations between hot streams i and cold streams j. Although an increase of possible combinations means that a better optimal solution/global optimum may be found, the searching process demands a considerable time resource as arisen by a high number of binary variables and non-convex terms, and thus become unacceptable for industrial applications. If a near global optimum is acceptable, a small stage number can be explored to reduce computational resources. Hence, an appropriate stage number is determined through an analysis in case studies based on trade-offs between time resource and solution quality. The following assumptions are made:

a. Constant stream flow rate, heat transfer coefficients and thermodynamic properties.

b. Split stream which contains heat exchangers in series and by-pass case are not considered.

c. Utility inlet and outlet temperatures are given.

d. Utilities are used only to target the final stream.

Objective function

The objective of this optimization model is to minimize the total annual cost (TAC) of a heat exchanger network associated with heat exchangers and utility cost. Heat exchanger capital cost is calculated by a short-cut model containing the sum of unit fixed cost and area cost. The required heat transfer area and number of units are determined in the optimal HEN configuration and heat transfer duties, with the trade-off between capital and energy usage. The objective function is formulated as:

$\begin{array}{l}TAC={\displaystyle \sum _i}{\displaystyle \sum _j}{\displaystyle \sum _k}{z}_{i,j,k}\left(C_f+C_a{A_{i,j,k}}^{\beta }\right)\\ +{\displaystyle \sum _i}zc{u}_i\left(C_f+C_a{Ac{u}_i}^{\beta }\right)+{\displaystyle \sum _j}zh{u}_j\left(C_f+C_a{Ah{u}_j}^{\beta }\right)\\ +{\displaystyle \sum _i}qc{u}_i{C}_{cu}+{\displaystyle \sum _j}qh{u}_j{C}_{hu}& (1)\end{array}$

where, C_f is the unit fixed cost coefficient and C_a is the heat exchanger area cost coefficient, while β is the area exponent. C_cu and C_hu are the annual cost for cold utility and hot utility energy flows. qcu_i and qhu_j are the cold utility and hot utility. z_{i, j, k}, zcu_i and zhu_j are binary variables defined as the existence of that match between hot and cold streams in stage k, hot streams and cold utilities, cold streams and hot utilities respectively. A_{i, j, k}, Acu_i and Ahu_j are heat transfer areas provided by heat exchangers, heaters and coolers.

Stage energy balance

Stage energy balances are required to achieve feasible heat transfer in each stage, as formulated in Eqs 2, 3:

$\begin{array}{ll}\left(T{h}_{i,k}-T{h}_{i,k+1}\right)Cp{h}_i={\displaystyle \sum _{j\in CP}}{q}_{i,j,k}i\in HP,k\in ST& (2)\end{array}$

$\begin{array}{ll}\left(T{c}_{j,k}-T{c}_{j,k+1}\right)Cp{c}_j={\displaystyle \sum _{i\in HP}}{q}_{i,j,k}j\in CP,k\in ST& (3)\end{array}$

where Th_{i, k} and Tc_{j, k} are defined as the temperature of a hot stream as it enters stage k and temperature of a cold stream as it leaves stage k, respectively; Cph_i and Cpc_j are the total heat capacities for hot and cold streams. q_{i, j, k} is the heat duty for each possible match.

Overall energy balance

The heat exchanged from an individual hot (cold) stream should be balanced with total heat duties transferred in all stages and utilities, as formulated in Eqs 4, 5:

$\begin{array}{ll}\left(Thi{n}_i-TThou{t}_i\right)Cp{h}_i={\displaystyle \sum _{j\in CP}}{\displaystyle \sum _{k\in ST}}{q}_{i,j,k}+qc{u}_{i}i\in HP& (4)\end{array}$

$\begin{array}{ll}\left(TTcou{t}_j-Tci{n}_j\right)Cp{c}_j={\displaystyle \sum _{i\in HP}}{\displaystyle \sum _{k\in ST}}{q}_{i,j,k}+qh{u}_{j}j\in CP& (5)\end{array}$

where Thin_i and TThout_i are the supply and target temperatures for hot stream i; Tcin_j and TTcout_j are the initial inlet and target temperatures for cold stream j.

To calculate the loads of hot and cold utility, Eqs 6, 7 are presented according to the last-stage temperatures of hot streams and the first-stage temperatures of cold streams, as shown in:

$\begin{array}{ll}\left(T{h}_{i,lastk}-TThou{t}_i\right)Cp{h}_i=qc{u}_{i}i\in HP& (6)\end{array}$

$\begin{array}{ll}\left(TTcou{t}_j-T{c}_{j,firstk}\right)Cp{c}_j=qh{u}_{j}j\in CP& (7)\end{array}$

The first k and the last k are the temperature location at k = 1 and k = K respectively. Notably, the total number of temperature locations = stage number (ST) + 1.

Energy balance of individual heat exchangers

Different from the iso-thermal mixing (within one stream and single stage, outlet temperatures for split streams are equal to stage outlet mixing temperature), energy balance of individual exchangers is applied to express the non-isothermal mixing heat transfer. As illustrated in Figure 2, for a hot stream (cold stream), the inlet temperature of a split stream that enters a heat exchangers is the same as the stage inlet temperature Th_{i, k} (Tc_{j, k+1}), while the outlet temperature for each split stream can be different according to the varying heat capacities and heat duties. Continues variables fh_{i, j, k} and fc_{i, j, k}(0 ≤ fh_{i, j, k} ≤ 1;0 ≤ fc_{i, j, k} ≤ 1) are defined as split fractions of hot stream i and cold stream j to denote the changed heat capacities.

FIGURE 2

Figure 2. Non-isothermal mixing at stage k.

The true outlet temperatures of exchangers can be calculated by exchanger energy balance, as formulated by Eqs 8, 9:

$\begin{array}{ll}\left(T{h}_{i,k}-Thou{t}_{i,j,k+1}\right)Cp{h}_{i}f{h}_{i,j,k}\ge q_{i,j,k}i\in HP,j\in CP,k\in ST& (8)\end{array}$

$\begin{array}{ll}\left(Tcou{t}_{i,j,k}-T{c}_{j,k+1}\right)Cp{c}_{j}f{c}_{i,j,k}\ge q_{i,j,k}i\in HP,j\in CP,k\in ST& (9)\end{array}$

where, Thout_{i, j, k} and Tcout_{i, j, k} are the actual outlet temperatures of exchangers before mixing at the end of a stage. In this energy balance, if the match does not exist, the split stream flow rate must be zero, as shown in Eqs 10, 11:

$\begin{array}{ll}f{h}_{i,j,k}\le z_{i,j,k}i\in HP,j\in CP,k\in ST& (10)\end{array}$

$\begin{array}{ll}f{c}_{i,j,k}\le z_{i,j,k}i\in HP,j\in CP,k\in ST& (11)\end{array}$

The previous approach by Pavão et al. (2017b), provides a heat transfer stage, in which the sum of split fractions of a single stream is always equal to 1 even if a match does not exist, which violates the energy balance constraints presented in Eqs 10, 11 when a stream fully bypasses a heat exchanger during the optimization. Therefore, this work allows streams to fully bypass a stage if no exchangers exist, defined through two variables, as follows:

$\begin{array}{ll}y{h}_{i,k}\{\begin{array}{cc}1& \text{if\hspace{0.17em}hot\hspace{0.17em}stream\hspace{0.17em}i\hspace{0.17em}bypass\hspace{0.17em}the\hspace{0.17em}stage\hspace{0.17em}k}\\ 0& \text{otherwise}\end{array}i\in HP,k\in ST& (12)\end{array}$

$\begin{array}{ll}y{c}_{j,k}\{\begin{array}{cc}1& \text{If\hspace{0.17em}cold\hspace{0.17em}stream\hspace{0.17em}j\hspace{0.17em}bypass\hspace{0.17em}the\hspace{0.17em}stage\hspace{0.17em}k}\\ 0& \text{otherwise}\end{array}j\in CP,k\in ST& (13)\end{array}$

Based on Eqs 12, 13, yh_{i, k} (yc_{j, k}) is 1 if stream i fully bypass stage k and there is no existence of a match. Otherwise it is equal to 0. Then, the constraints of Eqs14–19 to define feasible stream split fraction were proposed by Huang et al. (2012). For a hot stream i (cold stream j), if any matches exist in stage k, yh_{i, k} (yc_{j, k}) is equal to 0, and the sum of hot stream split fraction fh_{i, j, k} (fc_{i, j, k}) is equal to 1. If hot stream i (cold stream j) completely bypasses this stage, it leads the sum of fh_{i, j, k} (fc_{i, j, k}) to be 0.

$\begin{array}{ll}{\displaystyle \sum _j}f{h}_{i,j,k}+y{h}_{i,k}=1i\in HP,k\in ST& (14)\end{array}$

$\begin{array}{ll}{\displaystyle \sum _i}f{c}_{i,j,k}+y{c}_{j,k}=1j\in CP,k\in ST& (15)\end{array}$

$\begin{array}{ll}y{h}_{i,k}+z_{i,j,k}\le 1i\in HP,j\in CP,k\in ST& (16)\end{array}$

$\begin{array}{ll}y{c}_{j,k}+z_{i,j,k}\le 1i\in HP,j\in CP,k\in ST& (17)\end{array}$

$\begin{array}{ll}{\displaystyle \sum _j}{z}_{i,j,k}+y{h}_{i,k}\ge 1i\in HP,k\in ST& (18)\end{array}$

$\begin{array}{ll}{\displaystyle \sum _i}{z}_{i,j,k}+y{c}_{j,k}\ge 1j\in CP,k\in ST& (19)\end{array}$

Constraints for flexible splits

To ensure feasible design, a HEN configuration should satisfy industrial requirements. Generally, if the number of splits are more than four, the HEN design is often considered as impractical because of increases operational and control issues. In the existing superstructure approaches, the split limit is generally based on the constraint of the number of total binary variables for a stream, as shown in Eqs 20, 22, which cannot restrict the number of total splits for a stream across all stages. For example, for hot stream i, consider using $\sum _{j\in CP}{z}_{i,j,k}$ and $\sum _{j\in CP}\sum _{k\in ST}{z}_{i,j,k}$ to represent splits in a single stage k and total stages, respectively. But for a condition, only one heat exchanger ($\sum _{j\in CP}{z}_{i,j,k}=1$) is placed at a single stage k, or some heat exchangers are placed at stages with one exchanger allocating in one stage ($\sum _{j\in CP}{z}_{i,j,k}=1and\sum _{j\in CP}\sum _{k\in ST}{z}_{i,j,k}>1$), which assigns a positive value (1 or more than 1) to total splits even if there is no split to be considered. If consider to use $\sum _{j\in CP}{z}_{i,j,k}-1$ and $\sum _k\left(\sum _{j\in CP}{z}_{i,j,k}-1\right)$ to correct the above two expressions, $\sum _{j\in CP}{z}_{i,j,k}-1$ is assigned to −1 when no match exists in a stage k, which leads to a negative number for actual total split.

To represent the number of total stream splits, this work proposes two linear constraints t by using continuous variables, as presented in Eqs 21,23:

$\begin{array}{ll}{\displaystyle \sum _{j\in CP}}{z}_{i,j,k}\le N_{si}i\in HP,k\in ST& (20)\end{array}$

$\begin{array}{l}{\displaystyle \sum _{j\in CP}}{\displaystyle \sum _{k\in ST}}{z}_{i,j,k}-\left({ST}_{max}-{\displaystyle \sum _{k\in ST}}y{h}_{i,k}\right)\\ +1\le N_{si}i\in HP& (21)\end{array}$

$\begin{array}{ll}{\displaystyle \sum _{i\in HP}}{z}_{i,j,k}\le N_{sj}j\in CP,k\in ST& (22)\end{array}$

$\begin{array}{l}{\displaystyle \sum _{i\in HP}}{\displaystyle \sum _{k\in ST}}{z}_{i,j,k}-\left({ST}_{max}-{\displaystyle \sum _{k\in ST}}y{c}_{j,k}\right)\\ +1\le N_{sj}j\in CP& (23)\end{array}$

where, N_si and N_sj are demanded total splitting number for each hot stream (cold stream).

Therefore, Eq. 21 uses $\sum _{j\in HP}{z}_{i,j,k}-\left(1-y{h}_{i,k}\right)$ to express the total split number in stage k. If there are multiple matches in stage k, bypass variable yh_{i, k} is equal to 0 and $\sum _{j\in CP}{z}_{i,j,k}-1$ become the actual split number in stage k. If $\sum _{j\in HP}{z}_{i,j,k}=0$, hot stream i fully bypasses this stage and yh_{i, k} = 1, and the left side of Eq. 21 gives a value of 1 that expresses the basic stream without branch. Now, the correction further derives real total stream split number by $\sum _{j\in CP}\sum _{k\in ST}{z}_{i,j,k}-\left({ST}_{max}-\sum _{k\in ST}y{h}_{i,k}\right)$, and value 1 is added at the left side of Eq. 21 as a basic stream number. Similar to hot streams, Eq. 23 corresponds to the constraint for cold streams to define total split numbers.

Assessment of temperatures

The inlet temperatures, Thin_i, Tcin_j, are assigned to the first stage temperature of hot streams and the last stage temperature of cold streams, as:

$\begin{array}{ll}Thi{n}_i=T{h}_{i,firstk}i\in HP& (24)\end{array}$

$\begin{array}{ll}Tci{n}_j=T{c}_{j,lastk}j\in CP& (25)\end{array}$

The same inlet temperatures also are given to split stream outlet temperatures of hot stream Thout_{i, j, k} and cold stream Tcout_{i, j, k} respectively, through Eqs 26, 27.

$\begin{array}{ll}Thi{n}_i=Thou{t}_{i,j,firstk}i\in HP,j\in CP,k\in ST& (26)\end{array}$

$\begin{array}{ll}Tci{n}_j=Tcou{t}_{i,j,lastk}i\in HP,j\in CP,k\in ST& (27)\end{array}$

Area calculation

To calculate heat exchanger area, logarithmic mean temperature difference (LMTD) and overall heat transfer coefficients (U_{i, j}) are required. Equation 28 is used to calculate the overall heat transfer coefficient:

$\begin{array}{ll}{U}_{i,j}=1/(\frac{1}{h{h}_i}+\frac{1}{h{c}_j})i\in HP,j\in CP& (28)\end{array}$

where, hh_i and hc_j are the film heat transfer coefficients for hot streams and cold streams, respectively. For heaters and coolers, Eqs 29, 30 are formulated:

$\begin{array}{ll}Uh{u}_j=1/(\frac{1}{hhu}+\frac{1}{h{c}_j})j\in CP& (29)\end{array}$

$\begin{array}{ll}Uc{u}_i=1/(\frac{1}{h{h}_i}+\frac{1}{hcu})i\in HP& (30)\end{array}$

where, hhu and hcu are the film heat transfer coefficients for hot utilities and cold utilities, respectively. Constant film heat transfer coefficients are assumed in this work. Fouling resistance and detailed exchanger geometries are not considered here.

Several papers proposed approximations of LMTD widely used in HEN synthesis to avoid the numerical difficulties caused by equal approach temperatures. Huang et al. (2012) presented and evaluated these available LMTD approximations based on different design cases, and indicated the LMTD approximation by Paterson (1984) can lead to cost-efficient solutions, as shown in Eq. 31. But it is known that Eq. 31 may overestimate the exchanger area, thus Eq. 32 is used to correct LMTD according to the approach temperatures that are obtained from the optimal solution.

$\begin{array}{l}LMT{D}_{i,j,k}={\left[d{t}_{i,j,k}·d{t}_{i,j,k+1}\frac{\left(d{t}_{i,j,k}+d{t}_{i,j,k+1}\right)}{2}\right]}^{1/3}\\ j\in CP& (31)\end{array}$

$\begin{array}{ll}LMT{D}_{i,j,k}=\frac{d{t}_{i,j,k}-d{t}_{i,j,k+1}}{ln\frac{d{t}_{i,j,k}}{d{t}_{i,j,k+1}}}i\in HP& (32)\end{array}$

where, dt_{i, j, k} and dt_{i, j, k+1} are approach temperatures at the hot end and cold end respectively. Based on the split stream temperatures, two constraints are given to formulate the approach temperatures:

$\begin{array}{l}d{t}_{i,j,k}\le T{h}_{i,k}-Tcou{t}_{i,j,k}+G_{i,j}\left(1-z_{i,j,k}\right)\\ i\in HP,\in CP,k\in ST& (33)\end{array}$

$\begin{array}{l}d{t}_{i,j,k+1}\le Thou{t}_{i,j,k+1}-T{c}_{j,k+1}+G_{i,j}\left(1-z_{i,j,k}\right)\\ i\in HP,j\in CP,k\in ST& (34)\end{array}$

where, G_{i, j} is the upper bound for these inequalities. When z_{i, j, k} is equal to 1, these constraints are activated, otherwise the upper bound G_{i, j} makes these be inactive. The value of G_{i, j} can be determined by given inlet and outlet streams temperatures:

$\begin{array}{l}{G}_{i,j}=max[0,\left(Tci{n}_j-Thi{n}_i\right),\left(Tci{n}_j-TThou{t}_i\right),\\ \left(TTcou{t}_j-Thi{n}_i\right),\left(TTcou{t}_j-TThou{t}_i\right)]\\ i\in HP,j\in CP& (35)\end{array}$

For heaters (coolers), cold end (hot end) approach temperature is replaced by dthu_j (dtcu_i), as shown in Eqs 36, 37, and LMTDhu_j (LMTDcu_i) of heaters (coolers) can be calculated in Eqs 38, 39.

$\begin{array}{ll}dth{u}_j\le Thuout-T{c}_{j,firstk}j\in CP& (36)\end{array}$

$\begin{array}{ll}dtc{u}_i\le T{h}_{i,lastk}-Tcuouti\in HP& (37)\end{array}$

$\begin{array}{l}LMTDh{u}_j={\left[\left(Thuin-TTcou{t}_j\right)dth{u}_j\frac{\left(Thuin-TTcou{t}_j+dth{u}_j\right)}{2}\right]}^{1/3}\\ j\in CP& (38)\end{array}$

$\begin{array}{l}LMTDc{u}_i={\left[\left(TThou{t}_i-Tcuin\right)dtc{u}_i\frac{\left(TThou{t}_i-Tcuin+dtc{u}_i\right)}{2}\right]}^{1/3}\\ i\in HP& (39)\end{array}$

where, Tcuin and Tcuout are inlet and outlet temperatures of cold utilities. Thuin and Thuout are inlet and outlet temperatures of hot utilities. Now, the heat exchanger area can be calculated as:

$\begin{array}{ll}{A}_{i,j,k}=\frac{q_{i,j,k}}{LMT{D}_{i,j,k}{U}_{i,j}}i\in HP,j\in CP,k\in ST& (40)\end{array}$

For heaters (coolers), q_{i, j, k}, LMTD_{i, j, k} and U_{i, j} are replaced by qhu_i(qcu_i), LMTDhu_j (LMTDcu_i) and Uhu_j (Ucu_i). The formulated exchanger area (Eq. 40) directly substitutes the process variables, A_{i, j, k}, Acu_i and Ahu_j in the objective function (Eq. 1).

Feasibility constraints

To ensure a feasible HEN configuration, some process constraints reported by Yee and Grossmann (1990), and Huang et al. (2012) [basic feasibilities of temperatures, and feasible minimum approach temperature difference (MATD)] are adopted in this model, are described in Eqs 41–47:

$\begin{array}{ll}T{h}_{i,lastk}\ge TThou{t}_{i}i\in HP& (41)\end{array}$

$\begin{array}{ll}TTcou{t}_j\ge T{c}_{j,firstk}j\in CP& (42)\end{array}$

$\begin{array}{ll}T{h}_{i,k}\ge T{h}_{i,k+1}i\in HP,k\in ST& (43)\end{array}$

$\begin{array}{ll}T{c}_{j,k}\ge T{c}_{j,k+1}j\in CP,k\in ST& (44)\end{array}$

$\begin{array}{ll}d{t}_{i,j,k}\ge MATDi\in HP,j\in CP,k\in ST& (45)\end{array}$

$\begin{array}{ll}dth{u}_j\ge MATDj\in CP& (46)\end{array}$

$\begin{array}{ll}dtc{u}_i\ge MATDi\in HP& (47)\end{array}$

For non-isothermal mixing, the supply temperature of hot stream i (cold stream j) should be associated to actual exchanger outlet temperature at the first k. In addition, the stage inlet temperature of hot stream i must be no less than the split stream temperature that enters the next stage k + 1. Similarly, the split stream temperatures of cold stream j as it leaves stage k must be equal to or larger than the stage inlet temperature. Hence, additional constraints of temperatures are presented as follows:

$\begin{array}{ll}Thi{n}_i=Thou{t}_{i,j,firstk}i\in HP,j\in CP,k\in ST& (48)\end{array}$

$\begin{array}{ll}Tci{n}_j=Tcou{t}_{i,j,lastk}i\in HP,j\in CP,k\in ST& (49)\end{array}$

$\begin{array}{ll}T{h}_{i,k}\ge Thou{t}_{i,j,k+1}i\in HP,j\in CP,k\in ST& (50)\end{array}$

$\begin{array}{ll}T{c}_{j,k+1}\le Tcou{t}_{i,j,k}i\in HP,j\in CP,k\in ST& (51)\end{array}$

For heat transfer in split streams, the temperature drop must be 0 if a match does not exist. An upper bound of temperature drop is required to formulate these constraints, as presented:

$\begin{array}{l}Tuppe{r}_{i,j,k}=max\left(0,Thi{n}_i-Tci{n}_j-MATD\right)i\in HP,j\in CP,\\ k\in ST& (52)\end{array}$

Then, the upper bounds participate in Eqs 53, 54:

$\begin{array}{l}T{h}_{i,k}-Thou{t}_{i,j,k+1}-z_{i,j,k}Tuppe{r}_{i,j,k}\le 0i\in HP,j\in CP,\\ k\in ST& (53)\end{array}$

$\begin{array}{l}Tcou{t}_{i,j,k}-T{c}_{j,k+1}-z_{i,j,k}Tuppe{r}_{i,j,k}\le 0i\in HP,j\in CP,\\ k\in ST& (54)\end{array}$

The constraints associated with heat duties significantly affect the quality of HEN solutions. Some heat duty constraints have been exploited by Björk and Westerlund (2002) are based on the stream temperature drop that is bounded by [(Thin_i − TThout_i), (TTcout_j − Tcin_j)]. However, Huang et al. (2012) indicated that Tupper_{i, j, k} (Eq. 52) is more appropriate for non-isothermal mixing, because using MPAT is better for estimating the actual temperature of split streams. Thus, the upper bound of heat duties Qupper_{i, j, k} can be formulated according to Tupper_{i, j, k}:

$\begin{array}{l}{\displaystyle Quppe{r}_{i,j,k}=min[\left(Thi{n}_i-TThou{t}_i\right)Cp{h}_i,\left(TTcou{t}_j-Tci{n}_j\right)}\\ {\displaystyle Cp{c}_j,Tuppe{r}_{i,j,k}·min\left(Cp{h}_i,Cp{c}_j\right)]i\in HP,j\in CP,k\in ST}& (55)\end{array}$

Next, the logical constraints for Q_{i, j, k} are defined as follows:

$\begin{array}{ll}{q}_{i,j,k}-{z_{i,j,k}Qupper}_{i,j,k}\le 0i\in HP,j\in CP,k\in ST& (56)\end{array}$

$\begin{array}{ll}qh{u}_j-{zhu}_{j}Cp{c}_j(TTcou{t}_j-Tci{n}_j)\le 0j\in CP& (57)\end{array}$

$\begin{array}{ll}qc{u}_i-zc{u}_{i}Cp{h}_i\left(Thi{n}_i-TThou{t}_i\right)\le 0i\in HP& (58)\end{array}$

To enhance the heat duty constraints, some convex and concave constraints are given by Huang et al. (2012) to reduce redundant combinations:

$\begin{array}{l}{\displaystyle q_{i,j,k}\ge Cp{h}_i\left[T{h}_{i,k}-Thou{t}_{i,j,k+1}+({fh}_{i,j,k}-1){Tupper}_{i,j,k}\right]}\\ i\in HP,j\in CP,k\in ST& (59)\end{array}$

$\begin{array}{l}{\displaystyle q_{i,j,k}\ge Cp{c}_i\left[T{cout}_{i,j,k}-T{c}_{j,k+1}+({fc}_{i,j,k}-1){Tupper}_{i,j,k}\right]}\\ i\in HP,j\in CP,k\in ST& (60)\end{array}$

$\begin{array}{ll}{q}_{i,j,k}\le Cp{h}_i\left(T{h}_{i,k}-Thou{t}_{i,j,k+1}\right)i\in HP,j\in CP,k\in ST& (61)\end{array}$

$\begin{array}{ll}{q}_{i,j,k}\le Cp{c}_j\left(Tcou{t}_{i,j,k}-T{c}_{j,k+1}\right)i\in HP,j\in CP,k\in ST& (62)\end{array}$

Optimization approach

In this study, the proposed model is formulated as a non-convex MINLP problem, by applying the global solver BARON/GAMS to target a near global optimal solution under a specified time resource. The up-to-date version of BARON (GAMS. 38.2) can provide a rapid multi-start local search by solving bounding LPs to find a feasible solution that is taken as the initial value of the upper bound solution in the branch and bound algorithm, which can avoid consuming extensive computational resources in searching feasible solutions for large-scale problems in the past versions. Also, the solver introduces an option to solve subsequent branch and bound optimization by adjusting subordinate NLP and LP solvers automatically dealing with different problem features, which improves the time-efficiency in solving process.

Figure 3 illustrates the optimization route. The key block is to reduce the heat transfer stage (k ∈ ST) used to a large-scale problem, since obtaining a global optimum is difficult for a large-scale MINLP problem. To find a high-quality solution within an acceptable time, the required stage number should be tested and analyzed through results comparison. The LMTD equation is employed at the last step to avoid overestimated exchanger area obtained by using LMTD approximation. Furthermore, parallel computation can be applied and serviced by CPLEX embedded into BARON to deal with mixed integer trees using multiple CPU processors, which dramatically reduces computational time.

FIGURE 3

Figure 3. Block diagram of optimization route.

Case studies

In this section, five case studies are presented to demonstrate the performance of the new stage-wise superstructure for HEN synthesis using deterministic methods. Process details and cost coefficients are provided by the literature, as attached in Supplementary material. The lower bound of minimum approach temperature difference MATD is assumed to be 1°C to release a “free” search during the optimization process, in which the optimal process minimum approach temperature corresponds to the optimal trade-off between capital and energy. This work sets 10 hours as the maximum time resource for large-scale cases and records corresponding iteration numbers. Algorithm performance is demonstrated by evaluating the total annual cost obtained by optimization solvers vs. required iteration numbers (Figure 4). Notably, to show distinct movements, the figures only show the maximum iteration number that matches the least change of TAC, rather than the iteration at 10 hours. The software version is based on GAMS 38.2 operated on a computational platform Intel ® Core™, I7-11800H, 2.30GHz with 16 GB RAM. The computer operating threads in BARON are 10.

FIGURE 4

Figure 4. Iteration process in BARON/GAMS for case studie. (A) Case study 1. (B) Case study 2. (C) Case study 3.1. (D) Case study 3.2. (E) Case study 4. (F) Case study 5.

Case study 1

This case is a well-known aromatics plant problem summarized from an actual industrial case. It was originally proposed by Linnhoff and Ahmad (1990), including four hot streams and five cold streams. Hot oil and cooling water are adopted as utilities. They used Pinch Analysis to find a HEN configuration, but showed a relatively high cost. Subsequently, a series of papers proposed to solve the problem using different algorithms. Peng and Cui (2015) applied Simulated Annealing (SA) to obtain an optimal HEN (TAC = 2,935,000 $/year), in which no-splitting constraints are considered in their superstructure. Pavão et al. (2017b) employed non-isothermal mixing SWS for HEN synthesis through SA-RFO meta-heuristic algorithm, and investigated the effect of stage numbers vs solution quality. They found the best HEN configuration of using 4 stages, with 80 binary variables and CPU time of 380 s, which indicated that the appropriate stage number is not always equal to max[i, j]. This problem was resolved by Pavão et al. (2018a) using an enhanced SWS allowing intermediate placements of heaters and coolers. They raised the stage number to 6 with an increased problem size of 183 binary variables and obtained a HEN solution (TAC = 2,909,906 $/year), in which an intermediate heater was installed on cold stream one rather than placed at the end of this stream. Xu et al. (2019) applied a relaxation strategy for fixed capital cost into a random walk algorithm with compulsive evolution to tackle this problem, then generated a HEN configuration with TAC = 2,917,682 $/year. Their method also placed the utility at an intermediate position rather than fixing it at the end of the stream. Rathjens and Fieg (2020) proposed a hybrid approach using GA coupled with structure identification and a change of the reference system (SIR), and found the best solution so far with TAC = 2,908,767 $/year. But their method needs initial guessed split fractions demanded by a structure identification procedure, which causes large time resource in solving large-scale problems.

In this work, 5 and 4 stages are used to handle the superstructure, which form a problem size with 109 discrete variables and 2,175 single equations, and 89 discrete variables and 1,759 single equations, respectively. The changes in network TAC vs. iteration number through the above two operating conditions are shown in Figure 4A. The better solution (TAC = 2,883,847 $/year) is found at iteration 23,940 by 5 stages, consuming CPU 2,205 s. Figure 5 and Table 1 present the details of the optimal HEN configuration and result comparison. Compared with the best reported answer so far, the present HEN requires total process heat transfer area of 16,960 m², which is smaller than 18,099 m² proposed in early work (Rathjens and Fieg, 2020). Although a high exchanger number is proposed, both the exchanger area cost and the fixed unit cost (exchanger number) contribute to the total capital cost. Therefore, the proposed solution can obtain a TAC saving of 0.86% with reduced area cost.

FIGURE 5

Figure 5. Optimal HEN solution for case study 1 (TAC: 2,883,847 $/year).

TABLE 1

Table 1. Results comparison for case study 1.

Case study 2

This case contains six hot streams, four cold streams, a hot utility (steam) and a cold utility (cooling water). Its synthesis was originally studied by Ravagnani et al. (2005), in which they applied the GA algorithm to synthesize the problem based on an optimal minimum approach temperature difference that is obtained through Problem Table from Pinch Analysis, then GA was used to design HEN at above the pinch and below the pinch separately. The best solution proposed in their work required the TAC of 5,672,821 $/year using 3,360 s. Subsequently, this case was dealt with by Yerramsetty and Murty (2008) and Khorasany and Fesanghary (2009), who employed a differential evolution algorithm and a combination of HS algorithm and SQP model, respectively, requiring the TAC of 5,666,756 $/year and 5,662,366 $/year. The solution was then improved by Pourfarhady Myankooh and Shafiei (2016) with 0.38% of TAC saving, in which they adopted an ant colony optimization algorithm in real number domain (ACOR) without splitting. Feyli et al. (2022) obtained the best reported solution (TAC = 5,624,661 $/year) so far by using GA to handle discrete variables as structure match variables, and a modified quasi-linear programming (MQLP) to optimize continuous variables.

In this work, this problem is solved through 4, 5 and 6 stages (max[i, j]), respectively, as shown in Figure 4B. As the solving process illustrated in this figure, for each operation the network TAC decreases in this work by BARON/GAMS from about 5,680,000 $/year to 5,630,000 $/year using 10,000 iterations, which shows that the local search function provided by BARON/GAMS is time-efficient to determine a feasible and good solution as a subsequent upper bound that differs from the best solution (TAC = 5,612,881 $/year) found through stages 4 with only 0.3%. The final HEN structure and result comparison are shown in Figure 6 and Table 2, respectively. As can be seen from the table, the proposed solution has lower TAC (0.21% TAC saving) than the best reported solution so far (Feyli et al., 2022).

FIGURE 6

Figure 6. Optimal HEN solution for case study 2 (TAC: 5,612,881 $/year).

TABLE 2

Table 2. Results comparison for case study 2.

Case study 3

Heat integration of crude oil distillation is always an attractive problem since much energy can be saved. This case is a crude oil fractionation process proposed by Kim et al. (2017), including 11 hot streams, two cold streams (crude oil), one hot utility and one cold utility. The heat capacity for each stream provided by Kim et al. (2017) is corrected by Pavão et al. (2018b), who found a deviation of heat duty between that was calculated by stream temperature and calculated by energy balance from the solution. They indicated that this deviation can cause a structure mistake of HEN configuration. According to the work (Pavão et al., 2018b), an annualizing factor is set to 0.1.

This problem was originally studied by Kim et al. (2017) by a bounded method. Their work calculated a lower bound model first, and then the HEN MINLP superstructure was optimized though adding lifting partitioning equations to find an upper bound. To satisfy the practicality with industrial crude oil units, the total split numbers for crude oil stream, C1 and C2 should not exceed four, respectively. Kim et al. (2017) reported their best solution (TAC = 3,456,649 $/year) using CPU time 47,783 s. Pavão et al. (2018b) applied SA-RFO method to solve this problem with a new superstructure, in which cross flow, sub-splits and serial heat exchangers are considered. An optimal solution with lower TAC 3,391,066 $/year after CPU time of 22,630 s is found in their work. They noted that the solution obtained by SA-RFO method is better than that in the paper (Kim et al., 2017). Nevertheless, Pavão et al. (2018b) ignored the split limitation, which causes that the total split streams of C2 is larger than four and becomes impractical for industrial crude oil systems even if large stream splits may obtain better heat recovery performance. Therefore, to evaluate the performance of the proposed method, two solutions are conducted and compared with the results of Kim et al. (2017) and (Pavão et al., 2018b), respectively. Case study 3.1 is based on limited stream splits as introducing new constraints. Case study 3.2 is obtained through the present method directly without the limit of splits.

Case study 3.1: Limited split number

This problem is solved based on 4, 6, 8 and 11 stages (max[i, j]), and the MATD is set to 10°C to be consistent with the literature (Kim et al., 2017). Equations 20, 21 are updated in the proposed SWS in this case, as:

$\begin{array}{r}{\displaystyle \sum _{j\in CP}}{z}_{i,j,k}\le 4\\ {\displaystyle \sum _{j\in CP}}{\displaystyle \sum _{k\in ST}}{z}_{i,j,k}-\left(ST-{\displaystyle \sum _{k\in ST}}y{h}_{i,k}\right)+1\le 4\end{array}$

Figure 4C shows the corresponding trend of TACs in the investigation. After a rapid drop of TAC within 1,000 iterations, the best solution is found with eight stages, in which the programming includes 277 discrete variables and 5317 single equations, using CPU time of 21,183 s. Figure 7 shows the details of the HEN configuration. From Table 3, a better HEN solution is obtained with lower TAC (TAC = 3,455,358 $/year) and the required time resource is reduced significantly (saving 55% CPU time) than the work by Kim et al. (2017).

FIGURE 7

Figure 7. Optimal HEN solution for case study 3.1 (TAC: 3,455,358 $/year).

TABLE 3

Table 3. Results comparison for case study 3.

Case study 3.2: Free split number

Here, this model is optimized without the limits of splitting number to exploit a better solution by adding more possibilities of matches. This superstructure is also defined with 4, 6, 8 and 11 stages (max[i, j]), and the results are shown in Figure 4D. It is noticed that the best solution is led by 11 stages with iterations of 9,987, using 4,564 s, and the other three routes are trapped after initial falls. Figure 8 shows the HEN configuration. The result comparisons are presented in Table 3. According to these results, the TAC of 3,365,441 $/year is obtained, which is lower than the best reported solution so far (TAC = 3,391,066 $/year) with TAC redution of 0.76% and time resource saving of 79% compared with the best reported solution (Pavão et al., 2018b).

FIGURE 8

Figure 8. Optimal HEN solution case study 3.2 (TAC: 3,365,441 $/year).

Case study 4

A sugarcane biorefinery model, 1G/2G anhydrous ethanol and electricity production process was proposed by Oliveira et al. (2018). In their work, a multi-period operation for HEN synthesis was exploited through three periods (Oliveira et al., 2018), including an order of detailed process simulation, data extraction (details of process streams, cost information) and HEN synthesis. Because each of the three periods shows a changed fraction of hot steam-bagasse that is sent to 2G ethanol production, and the existence or inexistence of cold streams (soaking water, cellulose + lignin), thus period 1 operation, consisting of eight hot streams, seven cold streams, one hot utility and one cold utility is taken as a benchmark in the previous literature (Aguitoni et al., 2019), which is also used in this work. The optimal HEN configuration reported by Oliveira et al. (2018) was initially synthesized using SA-RFO algorithm with 600 s. Their HEN solution demands to install 19 exchangers in total, with total heat exchanger area of 29,152 m² and TAC of 12,424,312 $/year. Aguitoni et al. (2019) applied SA-DE to find a solution (TAC = 12,389,890 $/year) with lower utility cost, even with increased exchangers and 7,200 s for algorithm execution.

This work deals with the problem through 4, 6 and 8 stages (max[i, j]), respectively, as shown in Figure 4E. From this figure, the best solution is found with eight stages, where the TAC reduces rapidly from about 20,284,800 $/year to final 11,910,112 $/year using 759 iterations and CPU time of 4,867 s, then remaining unchanged afterwards. The corresponding HEN configuration is shown in Figure 9 and the result comparison is presented in Table 4. Compared with the best reported answer so far (Aguitoni et al., 2019), this work leads to 3.9% TAC reduction with lower utility cost, and the same number of units of 19 as reported by Oliveira et al. (2018). Due to industrial relevance of heat integration in biorefinery process, the present work makes a comparison of the developed deterministic method with scholastic-based algorithms, which can be developed to a multi-period-based SWS for practical considerations.

FIGURE 9

Figure 9. Optimal HEN solution for case study 4 (TAC: 11,910,112 $/year).

TABLE 4

Table 4. Results comparison for case study 4.

Case study 5

This case focuses on an industrial-scale aromatics plant problem taken from the Bandar Imam aromatics plant, proposed by Khorasany and Fesanghary (2009). The process includes six hot streams and 10 cold streams, two types of hot utility and one cold utility. The existing plant demands 122.19 MW hot utility and 524.72 MW cold utility, with about 8,856,412 $/year of TAC. In total, 18 heat exchangers are installed in the process.

Khorasany and Fesanghary (2009) optimized the split flow rates and heat duties simultaneously through a hybrid algorithm (Harmony Search and Sequential Quadratic method) and a lower TAC 7,435,740 $/year was obtained in that work with CPU time of 2,564 s. Zhaoyi et al. (2013) applied GA and PSO to conduct a HEN solution with and without stream splitting, in which the best configuration requires TAC of 7,361,190 $/year and CPU time of 2,860 s. Pavão et al. (2017b) continued to reduce the energy consumption to 34.21 MW for hot utility and 438.78 MW for cold utility. A hybrid approach SA/PSO was used in their work and 4022 s was needed to achieve the optimal solution. Then, Aguitoni et al. (2018) investigated the performance of GA and differential evolution used in this case. They carried out a sensitivity analysis with changes of stage number and the best HEN solution (TAC = 7,102,786 $/year) was found by three stages, after 15 executions (CPU 6,871 s). Their solution reduced the process area with an increase in utility consumption. However, their superstructure contained three stages to reduce the total binary variables in optimization. According to the method proposed by Yee and Grossmann (1990), dividing the HEN problem into three stages means a lot of proper matches were missed in solutions. Since they used 6,871 s, based on their method, a great time resource may be required if the stage number is increased. Bao et al. (2018) resolved this problem using a random walk algorithm with compulsive evolution and reduced the TAC to 6,861,111 $/year with significant utility saving (70.3% of hot utility saving). Xu et al. (2019) used a novel relaxation strategy consisting of a fixed capital cost using a heuristic rule to avoid sudden change of TAC during the optimization. They obtained the best solution with 6,832,780 $/year.

In the present work, 4, 6, 8, 10 stages (max[i, j]) are investigated respectively, as illustrated in Figure 4F. The design with eight stages shows a dramatic improvement of solution quality than others, and leads to the best solution with 76,745 iterations, consuming CPU time of 3,797 s. Figure 10 and Table 5 present the optimal HEN configuration and a list of previous solutions. Compared with the best reported solution (Xu et al., 2019), the present HEN shows a lower TAC 6,816,475 $/year (0.24% of TAC saving) with an efficient exchanger arrangement, thus a smaller requirement of heat transfer area.

FIGURE 10

Figure 10. Optimal HEN solution for case study 5 (TAC: 6,816,475 $/year).

TABLE 5

Table 5. Results comparison for case study 5.

Conclusions

The industrial application of heat exchanger network synthesis demands to handle large-scale problems, with acceptable time and cost-efficient solutions. This paper proposes an enhanced stage-wise superstructure for HEN synthesis. New variables and constraints are introduced in a modified stage-wise superstructure to reduce its computational difficulties without adding nonlinearities. Compared with the existing models, this method not only eliminates the assumption of iso-thermal mixing by adding heat exchanger energy balance constraints, but also can also optimize splitting numbers efficiently. The global MINLP solver BARON/GAMS is incorporated, using parallel computation for dealing with mixed integer trees in large-scale problems.

Consequently in all case studies carried out in this work, the proposed method shows better results than the existing approaches. A cost-efficient HEN solution can be obtained by the enhanced deterministic approach with acceptable time. Different from the existing stochastic algorithms that may have inconsistent solutions in different executions, the proposed deterministic approach can give better solutions, with lower TACs, including 0.86 % TAC saving for case 1, 0.21% TAC saving for case 2, 0.76% TAC saving for case 3.2, 3.9% TAC saving for case 4 and 0.24 % TAC saving for case 5, respectively.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

NZ proposed the presented idea. ZY developed the mathematical model and performed the computations, with the support of NZ and RS. NZ and RS co-supervised the project. All authors discussed the results and contributed to the final manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsus.2022.976717/full#supplementary-material

References

Adjiman, C. S., Androulakis, I. P., and Floudas, C. A. (1997). Global optimization of MINLP problems in process synthesis and design. Comput. Chem. Eng. 21, S445–S450. doi: 10.1016/S0098-1354(97)87542-4

CrossRef Full Text | Google Scholar

Aguitoni, M. C., Pavão, L. V., and Antonio da Silva Sá Ravagnani, M. (2019). Heat exchanger network synthesis combining simulated annealing and differential evolution. Energy 181, 654–664. doi: 10.1016/j.energy.2019.05.211

CrossRef Full Text | Google Scholar

Aguitoni, M. C., Pavão, L. V., Siqueira, P. H., Jiménez, L., and Ravagnani, M. A. D. S. S. (2018). Heat exchanger network synthesis using genetic algorithm and differential evolution. Comput. Chem. Eng. 117, 82–96. doi: 10.1016/j.compchemeng.2018.06.005

CrossRef Full Text | Google Scholar

Bao, Z., Cui, G., Chen, J., Sun, T., and Xiao, Y. (2018). A novel random walk algorithm with compulsive evolution combined with an optimum-protection strategy for heat exchanger network synthesis. Energy 152, 694–708. doi: 10.1016/j.energy.2018.03.170

CrossRef Full Text | Google Scholar

Bergamini, M. L., Grossmann, I., Scenna, N., and Aguirre, P. (2008). An improved piecewise outer-approximation algorithm for the global optimization of MINLP models involving concave and bilinear terms. Comput. Chem. Eng. 32, 477–493. doi: 10.1016/j.compchemeng.2007.03.011

CrossRef Full Text | Google Scholar

Björk, K.-M., and Westerlund, T. (2002). Global optimization of heat exchanger network synthesis problems with and without the isothermal mixing assumption. Comput. Chem. Eng. 26, 1581–1593. doi: 10.1016/S0098-1354(02)00129-1

CrossRef Full Text | Google Scholar

Chang, C., Liao, Z., Costa, A. L. H., and Bagajewicz, M. J. (2020). Globally optimal synthesis of heat exchanger networks. Part II: non-minimal networks. AIChE J. 66, e16264. doi: 10.1002/aic.16264

CrossRef Full Text | Google Scholar

Dorigo, M., and Caro, G. D. (1999). “Ant colony optimization: a new meta-heuristic,” in: Proceedings of the 1999 Congress on Evolutionary Computation-CEC99 (Cat. No. 99TH8406), 6-9 July, Vol. 2 (Washington, DC), 1470–1477. doi: 10.1109/CEC.1999.782657

CrossRef Full Text | Google Scholar

Faria, D. C., Kim, S. Y., and Bagajewicz, M. J. (2015). Global Optimization of the Stage-wise Superstructure Model for Heat Exchanger Networks. Ind. Eng. Chem. Res. 54, 1595–1604. doi: 10.1021/ie5032315

CrossRef Full Text | Google Scholar

Feyli, B., Soltani, H., Hajimohammadi, R., Fallahi-Samberan, M., and Eyvazzadeh, A. (2022). A reliable approach for heat exchanger networks synthesis with stream splitting by coupling genetic algorithm with modified quasi-linear programming method. Chem. Eng. Sci. 248, 117140. doi: 10.1016/j.ces.2021.117140

CrossRef Full Text | Google Scholar

Fieg, G., Luo, X., and Je zowski, J. (2009). A monogenetic algorithm for optimal design of large-scale heat exchanger networks. Chem. Eng. Process. Process Intensif. 48, 1506–1516. doi: 10.1016/j.cep.2009.10.003

CrossRef Full Text | Google Scholar

Hasan, M. M. F., Jayaraman, G., Karimi, I. A., and Alfadala, H. E. (2010). Synthesis of heat exchanger networks with nonisothermal phase changes. AIChE J. 56, 930–945. doi: 10.1002/aic.12031

CrossRef Full Text | Google Scholar

Holland, J. H. (1992). Adaptation in Natural and Artificial Systems: An Introductory Analysis With Applications to Biology, Control, and Artificial Intelligence. Cambridge, MA: MIT press. doi: 10.7551/mitpress/1090.001.0001

CrossRef Full Text | Google Scholar

Huang, K. F., Al-Mutairi, E. M., and Karimi, I. A. (2012). Heat exchanger network synthesis using a stagewise superstructure with non-isothermal mixing. Chem. Eng. Sci. 73, 30–43. doi: 10.1016/j.ces.2012.01.032

CrossRef Full Text | Google Scholar

Huang, K. F., and Karimi, I. A. (2013). Simultaneous synthesis approaches for cost-effective heat exchanger networks. Chem. Eng. Sci. 98, 231–245. doi: 10.1016/j.ces.2013.05.023

CrossRef Full Text | Google Scholar

Huo, Z., Zhao, L., Yin, H., and Ye, J. (2012). A hybrid optimization strategy for simultaneous synthesis of heat exchanger network. Korean J. Chem. Eng. 29, 1298–1309. doi: 10.1007/s11814-012-0007-2

CrossRef Full Text | Google Scholar

Kennedy, J., and Eberhart, R. (1995). “Particle swarm optimization,” in: Proceedings of ICNN'95 - International Conference on Neural Networks, 27 Nov.-1 Dec, Vol. 4 (Perth, WA), 1942–1948. doi: 10.1109/ICNN.1995.488968

CrossRef Full Text | Google Scholar

Khorasany, R. M., and Fesanghary, M. (2009). A novel approach for synthesis of cost-optimal heat exchanger networks. Comput. Chem. Eng. 33, 1363–1370. doi: 10.1016/j.compchemeng.2008.12.004

CrossRef Full Text | Google Scholar

Kim, S. Y., Jongsuwat, P., Suriyapraphadilok, U., and Bagajewicz, M. (2017). Global optimization of heat exchanger networks. part 1: stages/substages superstructure. Ind. Eng. Chem. Res. 56, 5944–5957. doi: 10.1021/acs.iecr.6b04686

CrossRef Full Text | Google Scholar

Kirkpatrick, S., Gelatt, C. D. Jr., and Vecchi, M. P. (1983). Optimization by simmulated annealing. Science 220, 671–680. doi: 10.1126/science.220.4598.671

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewin, D. R. (1998). A generalized method for HEN synthesis using stochastic optimization — II.: the synthesis of cost-optimal networks. Comput. Chem. Eng. 22, 1387–1405. doi: 10.1016/S0098-1354(98)00221-X

CrossRef Full Text | Google Scholar

Linnhoff, B., and Ahmad, S. (1990). Cost optimum heat exchanger networks−1. Minimum energy and capital using simple models for capital cost. Comput. Chem. Eng. 14, 729–750. doi: 10.1016/0098-1354(90)87083-2

CrossRef Full Text | Google Scholar

Linnhoff, B., and Hindmarsh, E. (1983). The pinch design method for heat exchanger networks. Chem. Eng. Sci. 38, 745–763. doi: 10.1016/0009-2509(83)80185-7

CrossRef Full Text | Google Scholar

Linnhoff, B., and Townsend, B. W. (1982). Designing Total Energy Systems.

Google Scholar

Oliveira, C. M., Pavão, L. V., Ravagnani, M. A. S. S., Cruz, A. J. G., and Costa, C. B. B. (2018). Process integration of a multiperiod sugarcane biorefinery. Appl. Energy 213, 520–539. doi: 10.1016/j.apenergy.2017.11.020

CrossRef Full Text | Google Scholar

Paterson, W. R. (1984). A replacement for the logarithmic mean. Chem. Eng. Sci. 39, 1635–1636. doi: 10.1016/0009-2509(84)80090-1

CrossRef Full Text | Google Scholar

Pavão, L. V., Costa, C. B. B., and Ravagnani, M. A. D. S. S. (2016). Automated heat exchanger network synthesis by using hybrid natural algorithms and parallel processing. Comput. Chem. Eng. 94, 370–386. doi: 10.1016/j.compchemeng.2016.08.009

CrossRef Full Text | Google Scholar

Pavão, L. V., Costa, C. B. B., Ravagnani, M. A. D. S. S., and Jiménez, L. (2017b). Large-scale heat exchanger networks synthesis using simulated annealing and the novel rocket fireworks optimization. AIChE J. 63, 1582–1601. doi: 10.1002/aic.15524

CrossRef Full Text | Google Scholar

Pavão, L. V., Costa, C. B. B., and Ravagnani, M. A. S. S. (2017a). Heat exchanger network synthesis without stream splits using parallelized and simplified simulated annealing and particle swarm optimization. Chem. Eng. Sci. 158, 96–107. doi: 10.1016/j.ces.2016.09.030

CrossRef Full Text | Google Scholar

Pavão, L. V., Costa, C. B. B., and Ravagnani, M. A. S. S. (2018a). An enhanced stage-wise superstructure for heat exchanger networks synthesis with new options for heaters and coolers placement. Ind. Eng. Chem. Res. 57, 2560–2573. doi: 10.1021/acs.iecr.7b03336

CrossRef Full Text | Google Scholar

Pavão, L. V., Costa, C. B. B., and Ravagnani, M. A. S. S. (2018b). A new stage-wise superstructure for heat exchanger network synthesis considering substages, sub-splits and cross flows. Appl. Therm. Eng. 143, 719–735. doi: 10.1016/j.applthermaleng.2018.07.075

CrossRef Full Text | Google Scholar

Peng, F., and Cui, G. (2015). Efficient simultaneous synthesis for heat exchanger network with simulated annealing algorithm. Appl. Therm. Eng. 78, 136–149. doi: 10.1016/j.applthermaleng.2014.12.031

CrossRef Full Text | Google Scholar

Pourfarhady Myankooh, Y., and Shafiei, S. (2016). A specific strategy for determination of feasible domain of heat exchanger networks with no stream splitting and its assessment by application of ACOR Algorithm. Appl. Therm. Eng. 104, 791–803. doi: 10.1016/j.applthermaleng.2016.05.115

CrossRef Full Text | Google Scholar

Rathjens, M., and Fieg, G. (2020). A novel hybrid strategy for cost-optimal heat exchanger network synthesis suited for large-scale problems. Appl. Therm. Eng. 167, 114771. doi: 10.1016/j.applthermaleng.2019.114771

CrossRef Full Text | Google Scholar

Ravagnani, M. A. S. S., Silva, A. P., Arroyo, P. A., and Constantino, A. A. (2005). Heat exchanger network synthesis and optimisation using genetic algorithm. Appl. Therm. Eng. 25, 1003–1017. doi: 10.1016/j.applthermaleng.2004.06.024

CrossRef Full Text | Google Scholar

Storn, R., and Price, K. (1997). Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11, 341–359. doi: 10.1023/A:1008202821328

CrossRef Full Text | Google Scholar

Toffolo, A. (2009). The synthesis of cost optimal heat exchanger networks with unconstrained topology, Appl. Therm. Eng. 29, 3518–3528. doi: 10.1016/j.applthermaleng.2009.06.009

CrossRef Full Text | Google Scholar

Xu, Y., Cui, G., Deng, W., Xiao, Y., and Ambonisye, H. K. (2019). Relaxation strategy for heat exchanger network synthesis with fixed capital cost. Appl. Therm. Eng. 152, 184–195. doi: 10.1016/j.applthermaleng.2019.02.054

CrossRef Full Text | Google Scholar

Yee, T. F., and Grossmann, I. E. (1990). Simultaneous optimization models for heat integration—II. Heat exchanger network synthesis. Comput. Chem. Eng. 14, 1165–1184. doi: 10.1016/0098-1354(90)85010-8

CrossRef Full Text | Google Scholar

Yerramsetty, K. M., and Murty, C. V. S. (2008). Synthesis of cost-optimal heat exchanger networks using differential evolution. Comput. Chem. Eng. 32, 1861–1876. doi: 10.1016/j.compchemeng.2007.10.005

CrossRef Full Text | Google Scholar

Zamora, J. M., and Grossmann, I. E. (1997). A comprehensive global optimization approach for the synthesis of heat exchanger networks with no stream splits. Comput. Chem. Eng. 21, S65–S70. doi: 10.1016/S0098-1354(97)87480-7

CrossRef Full Text | Google Scholar

Zhang, H., Cui, G., Xiao, Y., and Chen, J. (2017). A novel simultaneous optimization model with efficient stream arrangement for heat exchanger network synthesis. Appl. Therm. Eng. 110, 1659–1673. doi: 10.1016/j.applthermaleng.2016.09.045

CrossRef Full Text | Google Scholar

Zhaoyi, H., Liang, Z., Hongchao, Y., and Jianxiong, Y. (2013). Simultaneous synthesis of structural-constrained heat exchanger networks with and without stream splits. Can. J. Chem. Eng. 91, 830–842. doi: 10.1002/cjce.21702

CrossRef Full Text | Google Scholar

Nomenclature

Mathematical model

Greek Symbols

β The area exponent.

Subscripts