A Framework for Flexible and Cost-E ﬃ cient Retroﬁt Measures of Heat Exchanger Networks

: Retroﬁtting of industrial heat recovery systems can contribute signiﬁcantly to meeting energy e ﬃ ciency targets for industrial plants. One issue to consider when screening retroﬁt design proposals is that industrial heat recovery systems must be able to handle variations, e.g., in inlet temperatures or heat capacity ﬂow rates, in such a way that operational targets are reached. Consequently, there is a need for systematic retroﬁtting methodologies that are applicable to multi-period heat exchanger networks (HENs). In this study, a framework was developed to achieve ﬂexible and cost-e ﬃ cient retroﬁt measures of (industrial) HENs. The main idea is to split the retroﬁtting processes into several sub-steps. This splitting allows well-proven (single period) retroﬁt methodologies to be used to generate di ﬀ erent design proposals, which are collected in a superstructure. By means of structural feasibility assessment, structurally infeasible design proposals can be discarded from further analysis, yielding a reduced superstructure. Additionally, critical point analysis is applied to identify those operating points within the uncertainty span that determine necessary overdesign of heat exchangers. In the ﬁnal step, the most cost-e ﬃ cient design proposal within the reduced superstructure is identiﬁed. The proposed framework was applied to a HEN retroﬁt case study to illustrate the proposed framework.


Introduction
After the Paris Agreement 2015, decarbonization of industry and energy-intensive industrial processes has been proposed as one of the key measures to limit the increase of the global average temperature well below 2 • C. Within decarbonization of industry, waste heat recovery is dedicated to play a major role.
Systematic approaches to design and synthesize heat recovery systems, especially heat exchanger networks (HENs), have been subject to research since the early 1980s with the introduction of the graphical Pinch analysis method [1]. In addition to the graphical Pinch analysis method, mathematical programming has been well established as a complementary approach to HEN synthesis. Both sequential approaches [2] and simultaneous approaches [3], where targeting and network synthesis are performed separately and simultaneously, respectively, have been developed and reported in the literature. In addition, hybrid methods, which combine graphical methods and mathematical programming [4], have been suggested. Among others, Axelsson et al. [5], Bengtsson et al. [6], as well as Escobar and Trierweiler [7] have applied graphical and mathematical methods to case studies, including industrial applications, and reported results and comparisons. Furthermore, Papalexandri et al. [32] proposed an iterative scheme between the developed multi-period MI(N)LP and a flexibility analysis subproblem to identify the retrofit alternative, which is operable over a priori defined variations and yields the minimum total annualized cost. Another retrofitting methodology was reported by Kang and Liu [31], who introduced a two-step method to achieve retrofit design proposals that can operate cost efficiently at multiple periods. In a first step, the multi-period HEN synthesis model as it is used in greenfield design problems is employed. In the second step, existing exchangers are relocated in order to meet the required area demands identified in step one.
Although the above-mentioned retrofit approaches account for variations and uncertainty in operating data, there are several reasons why these approaches are not easily applied to large and complex industrial applications. Commonly, complexities, such as splitting or recirculation of streams, are present in industrial heat recovery systems, which cause nonlinearities in the mathematical formulation of these systems. Therefore, the resulting multi-period MINLP, including retrofit design proposals incorporated as hyper structures as proposed in [32], becomes complex, which inevitably leads to difficultness in finding feasible solutions. Additionally, the flexibility analysis subproblem for large and complex industrial applications is difficult to solve even with state-of-the-art global MINLP solvers. Moreover, depending on the problem complexity itself, the solution produced by common multi-period HEN synthesis formulations and the current network layout may differ essentially. Consequently, relocating existing exchangers to meet identified area demands as proposed in [31] can easily result in inefficient trial and error procedures.
The above reported methodologies are based on mathematical programming, which coincides with the benefits of mathematical programming when dealing with variations in operating data. However, in comparison to these methodologies, graphical insights-based methods have great potential to be applied to large and complex applications due to their beneficial user interaction. Recently, different strategies have been reported in the literature to utilize approaches based on graphical insights in the design process of multi-period HEN retrofitting case studies. A common strategy to deal with variations in operating conditions from a Pinch analysis perspective is to develop different retrofit proposals for a number of selected sets of operating conditions (e.g., annual, seasonal, and monthly average values) [8]. The different design proposals are then evaluated and may be combined to achieve an, over all considered operating points, operable and energy-efficient retrofit proposal. Another strategy is to develop different retrofit proposals by employing graphical insights-based retrofitting methods (e.g., Bridge analyses [22] or advanced composite curves [33]) for a specific nominal point and analyze the network's response to variations to obtain insights and identify the best performing retrofit design proposal [34]. Additionally, Langner et al. [35] proposed to decouple the design and analysis steps in retrofitting processes. By means of this decoupling, well-proven (single period) retrofit design methodologies can be utilized to generate different design proposals, which thereafter are evaluated and adjusted in a trial and error manner with respect to flexibility and energy efficiency [35].
In conclusion, approaching a retrofit problem subject to variations in operating data from a Pinch analysis perspective relies to a large extent on trial and error procedures as different design proposals must be evaluated and manually combined. Good results can be obtained, but the design process itself is often very inefficient due to the trial and error character of the respective strategies.
As pointed out, systematic methodologies for retrofitting industrial HENs subject to variation in operating data are needed, but the available methodologies fulfill this demand only partly. Therefore, we propose a new framework in this paper to achieve flexible and cost-efficient retrofit measures by combining the beneficial designer interaction of graphical approaches (e.g., Pinch based) at an early stage in the design process with the efficiency of mathematical programming. By means of this combination, inefficient trial and error procedures can be avoided. The proposed framework is based on single-period (e.g., Pinch based) retrofit methodologies, (structural) flexibility analysis, critical point analysis, and multi-period optimization. The proposed framework is outlined in Section 3 of this paper. In the following section, the theoretical background of (structural) flexibility analysis and critical point analysis is provided.

(Structural) Flexibility Analysis
Flexibility analysis of HENs has been subject to research since the early 1980s. In 1982, Marselle et al. [36] first introduced the concept of resilient HENs with respect to a certain disturbance range in the inlet conditions. In 1985, Saboo et al. [37] introduced a resilience index to quantify the resilience of HENs. In the same year, Swaney and Grossmann [38] extended the concept of the resilience index to a flexibility index, which is applicable not only to HENs but also to chemical processes in general. Both the resilience and the flexibility index indicate the maximum disturbance range in which inlet conditions may vary while at the same time achieve feasible operation. This maximum disturbance range can be interpreted as a hyperrectangle in the space of the varying inlet conditions. In this context, both indices are defined as the ratio between the largest scaled hyperrectangle within the feasible region and the hyperrectangle defined by an expected disturbance range. Therefore, feasibility is achieved if the respective index is larger or equal to 1. For two varying inlet conditions, this can be visualized as in Figure 1. The largest scaled rectangle within the feasible region can mathematically be expressed by the following set of equations (δ corresponds to flexibility/resilience index): at an early stage in the design process with the efficiency of mathematical programming. By means of this combination, inefficient trial and error procedures can be avoided. The proposed framework is based on single-period (e.g., Pinch based) retrofit methodologies, (structural) flexibility analysis, critical point analysis, and multi-period optimization. The proposed framework is outlined in Section 3 of this paper. In the following section, the theoretical background of (structural) flexibility analysis and critical point analysis is provided.

(Structural) Flexibility Analysis
Flexibility analysis of HENs has been subject to research since the early 1980s. In 1982, Marselle et al. [36] first introduced the concept of resilient HENs with respect to a certain disturbance range in the inlet conditions. In 1985, Saboo et al. [37] introduced a resilience index to quantify the resilience of HENs. In the same year, Swaney and Grossmann [38] extended the concept of the resilience index to a flexibility index, which is applicable not only to HENs but also to chemical processes in general. Both the resilience and the flexibility index indicate the maximum disturbance range in which inlet conditions may vary while at the same time achieve feasible operation. This maximum disturbance range can be interpreted as a hyperrectangle in the space of the varying inlet conditions. In this context, both indices are defined as the ratio between the largest scaled hyperrectangle within the feasible region and the hyperrectangle defined by an expected disturbance range. Therefore, feasibility is achieved if the respective index is larger or equal to 1. For two varying inlet conditions, this can be visualized as in Figure 1. The largest scaled rectangle within the feasible region can mathematically be expressed by the following set of equations (δ corresponds to flexibility/resilience index): In Figure 1, the largest scaled rectangle within the feasible region is smaller than the rectangle according to the expected variations. Thus, feasibility with respect to the expected variations is not achieved (flexibility or resilience index is smaller than 1). In both index formulations, the physical performance of the HEN or the chemical process is described by the following set of constraints: In Figure 1, the largest scaled rectangle within the feasible region is smaller than the rectangle according to the expected variations. Thus, feasibility with respect to the expected variations is not achieved (flexibility or resilience index is smaller than 1). In both index formulations, the physical performance of the HEN or the chemical process is described by the following set of constraints: where d is the vector of design variables, x corresponds to the state variables, z is used for the control variables, and the varying inlet conditions or uncertain parameters are depicted by Θ (see [37,38] constraints ensure a minimum temperature difference as well as that heat transfer is only possible from hot to cold streams. Additionally, design constraints of HEXs (in form of heat transport equations) may also be considered as operational inequality constraints. Feasibility is achieved when all constraints i ∈ I and j ∈ J are satisfied at the point of operation. With both formulations, it is possible to describe the resilience or flexibility for convex problems in which the solution lies at a vertex point of the largest scaled (hyper-)rectangle within the feasible region (see Figure 1). In 1987, Grossmann and Floudas [39] developed an active set approach to guarantee a global solution of the flexibility index problem also for some non-convex problems. More recently, Li et al. [40] suggested a framework to calculate the flexibility index by means of an alternating direction matrix embedded in a simulated annealing algorithm. Furthermore, Zhao and Chen [41] proposed to explicitly calculate the shape of the uncertainty space via cylindrical algebraic decomposition and quantifier elimination.
In the context of flexibility analysis of HENs, the term structural flexibility is used to define the set of constraints that are included in the flexibility analysis. Often, it is distinguished between structural constraints (e.g., heat and mass balances) and design constraints (e.g., heat transport equations of HEX). Marselle et al. [36] distinguished, for example, between a resilient network structure and a resilient network itself. Marselle et al. [36] further suggest that a network structure is resilient if it remains feasible for the specified disturbance range independent of the HEX areas (i.e., HEX areas are not specified). A resilient network structure is, thus, the premise for a resilient network, which remains feasible for the specified disturbance range for specified HEX areas [36]. This definition of structural resilience was later used by Li et al. [42] to describe the structural flexibility of HENs. In accordance with the literature, our work distinguishes between the structural feasibility of a design and the (general) feasibility of a design. Both can be assessed by solving the flexibility index problem. In the case of structural feasibility assessment, design constraints are discarded, and only structural constraints are considered, i.e., the heat transferred by HEXs is not limited by design characteristics. In the case of (general) feasibility assessment, design constraints are included, i.e., the heat transferred by HEXs is limited by design characteristics (i.e., the installed surface area and the overall heat transfer coefficient of exchangers).

Multi-Period Design Problem and Critical Point Analysis
In a multi-period design problem, independently of the chosen approach, sufficiently many sets of operating points need to be considered to guarantee the flexibility and cost efficiency of the achieved solution. However, with an increasing number of sets of operating points, the complexity of the problem increases. Especially for large-scale industrial applications, the computational capacity can be a limiting factor. To overcome these difficulties, different works in the literature deal with the reduction of the sets of operating points, e.g., by sensitivity analyses [43] or by the identification of critical operating points [44] of a given design proposal. The main idea of these approaches is to identify those variations (within a given uncertainty span) of the uncertain parameters, which require the largest equipment size and are, thus, critical for the design and operation of the HEN. This implies that the HEN itself must be structurally feasible and the feasibility may only be limited by design characteristics, i.e., HEX surface area and heat transfer coefficients (see Section 2.1). The identified critical operating points can be utilized to find the necessary overdesign of HEXs in a fixed HEN structure. By including the identified critical points in the design problem, feasibility for the entire uncertainty span can be achieved. Additionally, representative operating points/scenarios may be identified, e.g., by the screening of historical data, and included in the design problem to achieve a cost-efficient solution (cost optimal for the set of critical points and representative operating points/scenarios). In conclusion, the design problem subject to critical and representative operating points/scenarios can be expressed as the following multi-period MI(N)LP optimization problem: In (5), TAC represents the total annualized cost of the HEN design, which is described by the annual operating cost c operating in €/year (utility cost) of a number of defined representative operating points/scenarios s ∈ S, which are weighted with their normalized duration factor w s (time period represented by operating point divided by entire time period), and the investment cost c investment in € (e.g., HEX investment cost), which are annualized with a given capital recovery factor CRF. The set of equality constraints is depicted by h i with i ∈ I (heat and mass balances) and the set of inequality constraints is represented by g j with j ∈ J (temperature and other operational restrictions). The set of design constraints g d depends on the set of the non-negative design variables d with d ∈ DV. In h i , g j , and g d , x is the vector of the state variables, z corresponds to the control variables, and the varying inlet conditions or uncertain parameters are depicted by θ. Consequently, in order to find the optimal value for TAC (for the respective operating points), the degrees of freedom of the optimization problem are the non-negative design variables d with d ∈ DV and the control variables (of the HEN) z. The set of constraints must be satisfied at all representative operating points/scenarios s ∈ S and at all identified critical operating points c ∈ CP, i.e., for certain, previously identified, fixed combinations of values for the uncertain parameters. In (5), the set of constraints must be satisfied at each operating point present in the union of the two sets S and CP: op ∈ (S ∪ CP). This way, a flexible and cost-efficient design is achieved.
To identify those variations (within a given uncertainty span) that require the largest equipment size, sensitivity analyses can be used if monotonic correlations between uncertain parameters and design variables exist. One example is the correlation between the uncertain heat transfer coefficient and the surface area of a HEX-the maximum area of the HEX is obtained at the smallest value of the heat transfer coefficient within the uncertainty span [45]. However, it can be assumed that not many parameters follow these clear correlations and critical operating points can be determined following the procedure introduced by Pintarič and Kravanja in [44], which is based on the following idea: "The main idea is to identify points with the largest values of design variables at optimum objective function. This may be achieved by maximizing design variables one by one, while allowing uncertain parameters to obtain any value between the specified bounds, and simultaneously minimizing cost function." [44] (p. 1607) This implies a max-min problem for each design variable d ∈ DV, i.e., maximizing the design variable of interest d i while minimizing the cost function C(x, z, d, θ), which is stated in problem (6). Compared to the cost function in (5) (TAC), in the cost function of (6) (C(x, z, d, θ)), the operating cost of the representative operating points is not included: To solve the max-min problem, Pintarič and Kravanja [44] developed different formulations: The Karush-Kuhn-Tucker (KKT-) formulation; • The two-level formulation; and • The approximate one-level formulation.
By means of these formulations, possible candidates for critical points are identified. In a second step, a set covering algorithm is applied to the identified candidates to merge them into a final set of critical points. Both the two-level and the approximate one-level formulation approximate the solution of the KKT-formulation. Moreover, the approximate one-level formulation is dependent on a heuristically chosen Big-M parameter and solutions can be very sensitive to this parameter. Additionally, if the system is convex, all critical points are vertices of the uncertain space. These critical vertices can be identified by solving problem (5), excluding the representative operating points, sequentially at all vertices. Those vertices at which the maximum values for each design variable are obtained are critical [44].

Methodology
In order to achieve retrofit design proposals that can operate cost efficiently at multiple operating points, a framework was developed. The framework is shown in Figure 2 and will be outlined in the following subsections. The proposed framework utilizes single-period and well-proven retrofit design methodologies but also "experience-based" retrofit design proposals can be considered. Design proposals that are structurally infeasible are discarded after a structural feasibility assessment. To avoid suboptimal solutions due to unnecessary overdesign of equipment, the proposed framework is based on the identification of critical points. Cost efficiency for the entire operating period is ensured, by considering representative operating points besides the identified critical points when solving the eventual design problem. A final feasibility check is performed to ensure that all critical points have been identified (i.e., the set of critical points is complete); otherwise, the new critical point, which is identified during this check, is added to the set of critical point(s) and the design problem is solved again.

HEN Retrofit Design Proposals for Single Operating Point-Generation of Superstructure
In a first step, retrofit design measures are derived using single-period and well-proven retrofit design methodologies. Additionally, "experience-based" design proposals may be considered, which can be based on, e.g., already known operability issues. By means of this, the designer is actively involved in the design process, which can be beneficial in industrial applications. The derived design changes are collected in a superstructure, which eventually represents all design proposals that are considered.

Structural Feasibility Assessment-Reduction of Superstructure
In the next step, the structural feasibility of the identified proposals is analyzed. As outlined in Section 2.1, this can be done by calculating the flexibility index discarding all design constraints and considering only structural constraints. The flexibility index can be calculated using the active set approach developed by Floudas and Grossman [39]. However, the active set approach may result in impractical large problem formulations as industrial applications are often complex. Alternative formulations exist, such as the calculation of the flexibility index via the direction matrix [40] or via cylindrical algebraic decomposition and quantifier elimination [41]. However, the first alternative relies on the simulated annealing algorithm, which means that the global solution may not be found in a reasonable period of time. Moreover, applying cylindrical decomposition to industrial applications (which typically involve more than four HEXs and uncertain parameters) requires excessive computational capacity, which may not be available. An alternative is Monte Carlo network simulations with random variations within the uncertainty span while degrees of freedom are utilized to control operational targets and minimize utility cost in a similar way as it was proposed by Kachacha et al. in [46]. With an increasing number of uncertain parameters and complexity of the problem, computational capacity is the limiting factor for this approach. Although all listed methodologies for (structural) feasibility assessment have known drawbacks, it is assumed that structural feasibility assessment of HENs (also in more complex industrial applications) is possible. To reduce the probability of faulty results, the authors suggest applying several of the proposed methods and compare the achieved results. If design proposals can be identified that are structurally infeasible for at least some operating points within the uncertainty span (i.e., flexibility index smaller than 1), these design proposals are excluded from the superstructure, yielding a reduced superstructure.

HEN Retrofit Design Proposals for Single Operating Point-Generation of Superstructure
In a first step, retrofit design measures are derived using single-period and well-proven retrofit design methodologies. Additionally, "experience-based" design proposals may be considered, which can be based on, e.g., already known operability issues. By means of this, the designer is actively involved in the design process, which can be beneficial in industrial applications. The derived design changes are collected in a superstructure, which eventually represents all design proposals that are considered.

Structural Feasibility Assessment-Reduction of Superstructure
In the next step, the structural feasibility of the identified proposals is analyzed. As outlined in Section 2.1, this can be done by calculating the flexibility index discarding all design constraints and considering only structural constraints. The flexibility index can be calculated using the active set approach developed by Floudas and Grossman [39]. However, the active set approach may result in impractical large problem formulations as industrial applications are often complex. Alternative formulations exist, such as the calculation of the flexibility index via the direction matrix [40] or via cylindrical algebraic decomposition and quantifier elimination [41]. However, the first alternative relies on the simulated annealing algorithm, which means that the global solution may not be found in a reasonable period of time. Moreover, applying cylindrical decomposition to industrial applications (which typically involve more than four HEXs and uncertain parameters) requires

Determination of Critical Points
In the next step, critical operating points for each design proposal are identified using the different strategies for the identification of candidates for critical points and the set covering algorithm introduced by Pintarič and Kravanja [44]. The determination of critical operating points has been implemented in a methodology for synthesis of HENs with large numbers of uncertain parameters [45]. However, it has not yet been implemented in a retrofitting framework. The main difference in a retrofitting framework is that besides the investment cost of new HEXs, investment in existing equipment can also be considered, e.g., heat transfer enhancement or an increase of the area of existing HEXs, as well as structural changes, such as resequencing and repiping of existing HEXs. However, structurally, problems (5) and (6) are similar for a HEN synthesis and a HEN retrofitting problem.
It is worth mentioning that for the critical point analysis, each design proposal (which is present in the reduced superstructure) is considered individually, i.e., for each design proposal, this step yields an individual set of critical points. To allow a fast evaluation of the different design proposals, the different strategies for the identification of candidates for critical points and the set covering algorithm proposed Energies 2020, 13, 1472 9 of 24 in [44] were automated in this work in order to be applicable to a superstructure-based approach and industrial applications (see Supplementary Material).

Multi-Period Design Problem for Reduced Superstructure
With all critical points obtained, the multi-period design problem for the reduced superstructure and the representative operating points can be formulated. This problem is very similar to problem (5), with the exception that each design proposal p ∈ P may include a not generalizable number of individual equality constraints i ∈ I p , inequality constraints j ∈ J p , and design constraints/variables d ∈ DV p as well as an individual set of critical points c ∈ CP p , which results in an individual set of operating points op p ∈ S ∪ CP p at which the set of constraints of the respective design proposal must be satisfied. Therefore, problem (5) is reformulated: In problem (7), the set of individual design proposals of the reduced superstructure is represented by p ∈ P. Each design proposal is expressed via a binary variable y p and one and only one of these binary variables is forced to be 1, which will be the design proposal, which performs most cost efficiently at the representative operating points (and critical operating points). Instead of solving problem (7) simultaneously, it may also be solved proposal-wise, setting one y p to 1 and eventually comparing the TAC of all proposals to identify the most cost-efficient proposal of the reduced superstructure. In this context, it should be noted that all costs that can be associated with the respective retrofit design proposal must be included in the respective cost functions c operating,p and c investment,p . This is especially important for retrofitting projects of industrial HENs, which are often very interconnected, and which may cause increased cost in other operational units.

Feasibility Check
In a final step, the feasibility is checked by means of the methodologies presented in Section 3.2. This step is necessary to ensure that all critical points have been identified (i.e., the set of critical points is complete). As outlined in Section 2.1, the previously derived HEX design characteristics are included, i.e., design constraints and solution obtained by solving the multi-period design problem. Since including design characteristics is likely to increase the complexity of the problem, solutions obtained by the methodologies presented in Section 3.2 should be analyzed and it is advisable to compare the achieved results of different methodologies to reduce the probability of failures. If feasibility for the respective design cannot be guaranteed, the result of the feasibility check can provide useful insights. Besides the flexibility index, the feasibility check provides a point in the space of the uncertain parameters, which limits the flexibility. The vector between this point and the average operating point (see Figure 1) can be used to generate candidates of critical points, which can be added to the multi-period design problem. If the newly generated design is feasible, the critical point has been successfully identified; otherwise, the new flexibility limit is analyzed, and the step is repeated until a feasible design is obtained.

Illustrative Example
The methodology described in Section 3 is illustrated using a four-stream example. The example was adapted from Lal et al. [34] with some modifications for demonstration purposes and was used as an initial network in a retrofitting process. In Figure 3, the initial network, the average stream data, and six different retrofit proposals are presented. These retrofit proposals are partly inspired by design proposals obtained by applying the Bridge analysis methodology presented in [34]. The main objective of the different retrofit proposals was to reduce the utility consumption in the heaters of the HEN. Table 1 shows the hot utility demand of the initial network as well as the hot utility demand and the hot utility savings of the six retrofit proposals.
feasibility for the respective design cannot be guaranteed, the result of the feasibility check can provide useful insights. Besides the flexibility index, the feasibility check provides a point in the space of the uncertain parameters, which limits the flexibility. The vector between this point and the average operating point (see Figure 1) can be used to generate candidates of critical points, which can be added to the multi-period design problem. If the newly generated design is feasible, the critical point has been successfully identified; otherwise, the new flexibility limit is analyzed, and the step is repeated until a feasible design is obtained.

Illustrative Example
The methodology described in Section 3 is illustrated using a four-stream example. The example was adapted from Lal et al. [34] with some modifications for demonstration purposes and was used as an initial network in a retrofitting process. In Figure 3, the initial network, the average stream data, and six different retrofit proposals are presented. These retrofit proposals are partly inspired by design proposals obtained by applying the Bridge analysis methodology presented in [34]. The main objective of the different retrofit proposals was to reduce the utility consumption in the heaters of the HEN. Table 1 shows the hot utility demand of the initial network as well as the hot utility demand and the hot utility savings of the six retrofit proposals.   The initial network consists of two process-to-process HEX and four utility HEXs. In Table 2, the design characteristics of the existing process-to-process HEX are listed. For simplicity, it was assumed that the utility HEXs are not limited by size, i.e., any possible duty can be satisfied. For process-to-process HEX, a constant overall heat transfer coefficient (i.e., U-value) of 0.523 kW/(m 2 • C) was assumed.

Structural Feasibility Assessment
The proposed retrofit designs can be considered as different pathways in a superstructure, and in a next step, the structural feasibility of the different proposals was evaluated. For the evaluation, Energies 2020, 13, 1472 11 of 24 variations in the four inlet temperatures and heat capacity flow rates were assumed. The variations are shown in Table 3. The structural feasibility was evaluated by calculating the flexibility index using the active set approach developed by Floudas and Grossman [39]. Since all constraints i ∈ I and j ∈ J (see Section 2.1) are linear, the active set approach guarantees the globality of the solution [39]. In Table 4, the structural flexibility index for each retrofit proposal is listed. The flexibility analysis revealed that retrofit proposal 5 is structurally infeasible, i.e., the structural flexibility index is smaller than 1 (see Table 4), meaning that for some operating points within the uncertainty span, the target temperature of at least one stream cannot be reached. This retrofit proposal was therefore removed from the initial superstructure of retrofit proposals, yielding a reduced superstructure.

Economic Data and Identification of Critical Points
In order to identify the critical points of the different retrofit design proposals, design variables were defined. These design variables could be modified in order to achieve feasibility for the entire range of operating conditions. It was assumed that the area of the two existing units could be increased while new units could be designed freely. Furthermore, the structure of the objective function of the overall design problem was specified. All in all, the total annualized cost (TAC) of the design should be minimized. Consequently, different cost functions were specified: Additionally, representative operating points s ∈ S, their normalized duration factor w s , the capital recovery factor CRF, and the annual operating time t operating were specified, in order to achieve cost efficiency for the entire operating period. The structure of the objective function for the retrofit design proposals is given in Equation (8): minTAC = s∈S w s p CU Q CU,s + p HU Q HU,s * t operating + (Inv.cos t increase + Inv.cos t new ) * CRF. (8) Energies 2020, 13, 1472

of 24
As mentioned in Section 2.2, the cost function, which is minimized in problem (6) to identify critical points, does not include the operating cost of the representative operating points. The cost function that was used for determining the critical points is given in Equation (9): minC = p CU Q CU + p HU Q HU * t operating + (Inv.cos t increase + Inv.cos t new ) * CRF.
For simplicity, it was assumed that the increase of the two existing HEX units could be described by the same cost function. The cost for increasing a HEX in € is given in Equation (10) For new HEX units, it was also assumed that the same cost function is applicable for all possible new HEX units. The cost of a new HEX in € is given in Equation (11) A capital recovery factor of 0.1 (e.g., 15 years lifetime and an interest rate of 7%) and an annual operating time of 8200 h were assumed. Furthermore, for the operating cost, the data in Table 5 were used. It was further assumed that no additional costs are associated with the retrofit proposals. For each of the five retrofit design proposals, which together constitute the reduced superstructure, one set of critical points was identified. These sets of critical points are listed in Table A5 in Appendix B. The sets of critical points listed in Table A5 are complete, i.e., the designs based on these sets were proven to be feasible by means of flexibility analysis (see Section 4.3). To avoid numerical difficulties, the logarithmic mean temperature difference in the design constraints of the HEXs was approximated using Paterson's approximation (see, e.g., [47]). As the determination of the critical points was connected to difficulties that were neither reported in the literature nor have been experienced when the results of the available literature examples of critical point analysis of HENs (see [44,45]) were reproduced, a description and a discussion on the determination of the critical points can be found in Appendix A. In this context, modifications on the existing methodologies were suggested to be able to determine critical points for more complex HEN structures but also for retrofit studies of HENs. For more information, see Appendix A.

Solution of the Multi-Period Design Problem and Final Feasibility Check
In the final step, the multi-period design problem was solved. In order to calculate the TAC of the different design proposals, representative operating points and their respective normalized duration factors were defined. For the illustrative example, 11 representative operating points were assumed. Furthermore, it was assumed that each operating point had the same normalized duration factor. The different representative operating points and their normalized duration factors are given in Table A4 in Appendix B. As can be seen in Table A4, extreme values (which can be calculated with the values given in Table 3) were not considered as representative operating points since these values represent extreme situations, which are usually not representative for longer operating periods.
To reduce the computational complexity, the multi-period design problem was solved individually for each retrofit proposal. The problems were solved using the global solver BARON. Table 6 shows the TAC for the proposals in the reduced superstructure, the annual net savings (difference between the annual operating cost of the initial HEN and the TAC of the proposals), and the flexibility index. The annual operating cost of the initial HEN was used as a benchmark for the TAC of the different retrofit proposals and was calculated considering the cost data presented in Section 4.2. The annual operating cost of the initial network is 314,600 €/y. Additionally, the derived design values of the process-to-process HEXs of the different retrofit design proposals are shown in Tables A6-A10 in Appendix B. The results in Table 6 indicate that proposal MER is the most cost efficient. Table 6. Total annualized cost, annual net savings, and flexibility index of the different retrofit design proposals in the reduced superstructure. The flexibility index was calculated using the active set approach developed by Floudas and Grossman [39] and the results were verified by searching for the minimum direction matrix using simulated annealing as it is reported in [40]. In contrast to [40], the subproblem to find the maximum value for δ (maximum feasible variation of the uncertain parameters) in a given direction was obtained by utilizing BARON. The parameters that were used for the simulated annealing algorithm are shown in Table A11 in Appendix B.

Proposal
Introducing different representative operating points to the problem increases the problem size and complexity, which may cause longer CPU times to guarantee the globality of the obtained solution. The solutions presented above were obtained on an Intel i7-6600 2.6 GHz processor with 16.0 GB RAM in 1800 s. It is worth mentioning, that globality (within default EpsA range of BARON: 1e-6) could be achieved in the case of proposals 2, 3, and 4. However, it should be mentioned that the optimality gap for proposal 1 was less than 500 €/y while the largest remaining optimality gap was 37,000 €/y in case of proposal MER. Although the remaining optimality gap of the proposal MER is considerable in comparison to the possible savings, the found local solution was considered satisfactory as it indicates that proposal MER is the most cost efficient.

Discussion of the Results of the Illustrative Example
As mentioned previously, a discussion around the determination of the set of critical points for the different retrofit proposals of the illustrative example can be found in Appendix A.
Based on the results presented in Table 6, the proposal MER performs the most cost efficiently. This, however, relates to the used cost functions. Different cost functions will certainly lead to different results and another proposal may be most cost efficient. The influence of the cost functions can be illustrated by a comparison of the solutions of proposals 1 and MER. An analysis of the results of the multi-period design problem of proposals 1 and MER revealed that the TAC of proposal 1 consists of 23% annualized capital cost (and 77% operating cost) while the TAC of proposal MER consists of 56% annualized capital cost (and 44% operating cost). Thus, the used cost functions favor design proposals that allow more utilization of the HEX surface area. However, if the cost for the area is increased (e.g., increased parameters in investment cost functions), it can be assumed that design proposals that allow more utilization of utilities will be more cost efficient.
From Table 6, it can be seen that design proposal 2 has a significantly higher flexibility index than 1, which indicates overdesign of the HEX units. This overdesign results from the consideration of the chosen representative operating points. When the multi-period design problem for proposal 2 was solved considering only the identified set of critical points (and discard all operating cost), the flexibility index of the achieved design was calculated as 1. To check the influence of the representative operating periods and thereby the completeness of the identified sets of critical points, similar calculations were Energies 2020, 13, 1472 14 of 24 also made for the other design proposals in the reduced superstructure and the results are shown in Table 7. The results in Table 7 indicate that the sets of critical points shown in Table A5 are complete as it was possible to achieve feasible design characteristics (flexibility index is 1) respecting only the sets of critical points. As all operating costs were discarded, the objective function of these design problems represented only the total annualized capital cost of a certain design.

Conclusions
In this work, a novel framework for retrofitting multi-period HENs was developed and presented. The proposed framework divides the retrofit design process into five sub-steps, which allows for the combining of the beneficial designer interaction of graphical approaches (e.g., Pinch based) at an early stage in the design process with the efficiency of mathematical programming to derive flexible and cost-efficient retrofit measures. In this context, inefficient trial and error procedures are avoided. Additionally, by means of splitting the design process in five different sub-steps, the complexity of the sub-problems is decreased compared to the overall problem.
The proposed framework utilizes the designer interaction of graphical approaches to derive different design proposals, which can be very beneficial for large and complex HENs usually found in industrial heat recovery systems. Well-proven, single-period retrofit design methods, such as Bridge analysis, but also "experience-based" retrofit design proposals may be employed. As the single-period character of graphical retrofitting methodologies does not ensure the flexibility and cost-efficiency of the generated proposals for the entire operating period, different mathematical evaluation strategies are incorporated in the proposed framework. By means of structural feasibility assessment based on the calculation of the flexibility index, structurally infeasible design proposals can be identified and discarded from further analysis. Additionally, critical point analysis was suggested to identify those operating points within the uncertainty span that determine necessary overdesign of the HEXs to ensure feasibility of the HEN retrofit proposals for the entire operating range. In order to identify the most cost-efficient proposal, a multi-period MI(N)LP optimization problem was formulated, which considerers the identified critical points (for feasibility) as well as representative operating points (for cost-efficiency). A final feasibility check was suggested to ensure that the complete set of critical points has been identified.
The determination of critical points based on [44] was automated and the achieved results were discussed. Compared to the available literature examples, the HEN examples included in this study are more complex, which caused difficulties when determining critical points. Modifications of the existing methodologies (see Appendix A) were suggested in order to identify the complete set of critical points for the HEN retrofit example presented in this study. It is worth mentioning that by means of the above-mentioned final feasibility check, an incomplete set of critical points is recognized. Additionally, by means of the obtained results of the final feasibility check, the missing critical point(s) can be identified to be included in the multi-period MI(N)LP design problem.
The framework presented in this paper yields opportunities to increase heat recovery in applications operating at multiple periods, e.g., industrial applications. As decarbonization of industry is essential to limit the increase of the global average temperature well below 2 • C, systematic approaches to successfully retrofit industrial heat recovery systems are in demand more than ever. In comparison to existing approaches for retrofitting multi-period HENs, the presented framework splits the design process in five different sub-steps, which decreases the complexity of the sub-problems. The decreased complexity can be the decisive factor for a successful retrofit project when large and complex industrial applications are addressed. The novel contributions of this paper are summarized in the following list: • Development of a novel framework to achieve flexible and cost-efficient retrofit measures for HENs operating at multiple periods; • Implementation of critical point analysis in a retrofitting framework to achieve flexible retrofit solutions of HENs; • Possibility to employ well-proven, single-period retrofit design methods, e.g., advanced composite curves, Bridge analysis, etc., but also "experience-based" retrofit design proposals in a framework that guarantees flexibility and cost efficiency (with respect to investment and operating cost) for the entire operating period; • Possibility to compare design proposals of different single-period retrofit design methods in a superstructure-based mathematical program ensuring flexibility and cost efficiency with respect to investment and operating cost; • Automatization of the KKT-formulation, the two-level formulation, and the set-covering algorithm initially suggested by Pintarič and Kravanja [44] to identify critical points of different design proposals (see Supplementary Material); and • Extension of the set covering algorithm suggested by Pintarič and Kravanja [44] in order to handle more complex network structures (see Appendix A).

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Discussion on the Determination of Critical Points
As mentioned in Section 4, the determination of critical points for the illustrative example was connected to difficulties. These difficulties are described in the following subsections and explanation attempts, and future work is proposed. This section of the appendix is intended for the interested reader and concerns specific issues related to the mathematical formulation of the critical point analysis.
For the determination of critical points, the three different strategies and the set covering algorithm presented by Pintarič and Kravanja in [44] were considered. As previously mentioned, the approximate one-level formulation depends on a heuristically chosen Big-M parameter and solutions can be very sensitive to this parameter [44]. For this reason, the approximate one-level formulation was discarded from this study. It was identified that depending on the chosen strategy, different aspects need to be considered, which may relate to the structure of the HEN under investigation.

Appendix A.1. Two-Level Formulation
The two-level formulation approximates the solution of the KKT-formulation and may be impractical for large industrial applications since an iterative scheme is involved. The main idea of the two-level formulation is to decompose the max-min problem ((6) in Section 2.2) into two subproblems, which can be solved separately in an iterative procedure. The two subproblems are: • The lower (control-design) level problem; and • The upper (uncertainty) level problem.
At iteration k, in the lower (control-design) level problem, the solution of the minimization problem (min x,z,d C(x, z, d, θ)) is approximated for fixed values of the uncertain parameters (obtained from upper (uncertainty) level problem in iteration k-1). In the upper (uncertainty) level problem, the solution of the maximization problem (maxd i θ ) is approximated by fixing the control variables to the values obtained by the lower (control-design) level problem while the uncertain parameters are relaxed to continuous variables within their respective bounds. The iteration procedure is stopped if the values obtained for the design variable d i at the lower and upper level problem are within a previously defined range. By means of the two-level formulation, the operating point in the uncertainty span is obtained at which the maximum value for the design variable d i is achieved while the cost function C(x, z, d, θ) is minimized (i.e., a potential critical point). For further details concerning the iterative procedure, see [44].
The iterative scheme can be tedious for large-scale industrial problems with many design variables, control variables, and uncertain parameters. To overcome this problem, the iterative scheme was automated as part of this study (see Supplementary Material). In this automation, the lower (control-design) level problem was solved using BARON and the upper (uncertainty) level problem was solved using the GAMS interface of CONOPT or IPOPT to obtain marginal values. In order to apply the two-level formulation, control variables need to be determined. Industrial systems may be complex and very interconnected, which results in many different and complex control strategies.
A mathematical analysis of the model equations (h i with i ∈ I) of the design proposals, which form the reduced superstructure of the illustrative example, revealed that the proposals demand either two (proposals 3 and 4) or three control variables (proposals 1, 2, and MER). This allows different combinations of control variables and it is worth mentioning that the obtained sets of critical points were not equivalent with respect to the solution of the final design problem. This will be explained at the example of design proposal 3. If the control variables are limited to the duties of the utility exchangers (Q3, Q4, and Q5 in Figure 3), three control variable combinations are possible. When using these different combinations to obtain the critical points through the two-level formulation and the set covering algorithm described in [44], different sets of critical points were obtained. These different sets, the TAC of the designs based on the different sets, and the corresponding flexibility index are shown in Table A1. As can be seen in Table A1, with only one of the tested control variable combinations, it was possible to achieve a feasible design. Table A1. Results of the two-level formulation for design proposal 3. It was identified that in order to be applicable to HEN retrofitting problems, including the possibility to increase the size of existing HEX units, the lower (control-design) level problem of iteration k (A1) of the two-level formulation (compare [44]) needs to be modified. For a selected design variable d i , the lower (control-design) problem is relaxed the following way: Physically, this reformulation implies that existing HEX units can be bypassed, which may be necessary to ensure the feasibility of the lower (control-design) level problem for given values of the uncertain parameters (obtained by the upper (uncertainty) level problem of the previous iteration).
As mentioned previously, the two-level formulation approximates the solution of the KKT-formulation, which would explain discrepancies between the sets of critical points identified with the two-level and the KKT-formulation. However, to the authors' knowledge, the influence of differently chosen control variables on the sets of critical points (and thereby on the flexibility index of the corresponding design) when applying the two-level formulation was not yet reported in the literature. It is worth mentioning that similar observations were made for proposal 4, although the values for the TAC and the flexibility index were different compared to Table A1.
Furthermore, it was observed that with more complex network structures, additional difficulties may occur if the two-level formulation is applied. This will be explained at the example of proposal 1. In the case of proposal 1, three control variables need to be defined, which are fixed when solving the upper (uncertainty) level problem. Again, different combinations were tested, and the flexibility indices of the derived designs were calculated (see Table A2). For explanation of the variable names of the control variables, refer to Figure A1 in Appendix B. As can be seen from the flexibility indices in Table A2, with none of the tested combinations of control variables, a feasible design could be derived. It is worth mentioning that depending on the choice of the control variables, convergence between the lower (control-design) level problem and the upper (uncertainty) level problem could not be achieved for at least some of the design variables of proposal 1. Similar observations were made when the two-level formulation was used to identify the sets of critical points for proposals 2 and MER. Table A2. Results of the two-level formulation for design proposal 1. These observations could have several explanations: Another, not tested, combination of control variables must be chosen, the calculation of the marginal values of GAMS is not sufficiently accurate, the set covering algorithm is not globally valid, or the approximative character of the solution of the two-level formulation is too strong. When analyzing the complete set of critical points of proposal 1, which is shown in Table A5 in Appendix B, it was found that the actual critical point is [240.0, 190.0, 19.0, 120.0, 10.0, 20.0, 20.8, 32.0] (i.e., this point is sufficient to derive a feasible design for proposal 1). An indication for a not valid set covering algorithm would be if this critical point can be found as a combination of one of the above presented points, i.e., if a "0" element in one of the points could be replaced with another value to achieve the critical point. This, however, is not possible. This is not proof that the set covering algorithm is globally valid. It indicates, however, that in the case of proposal 1, the set covering algorithm is not responsible for missing the critical point and one of the other above listed reasons must be responsible. Future work should focus on identifying strategies or guidelines for the application of the two-level formulation to HENs.
As mentioned previously, the difficulties experienced with the two-level formulation are not reported in the literature. These difficulties did not appear when the results reported in the literature (application of two-level formulation to the HEN example; see [44]) were reproduced. It was assumed that the experienced difficulties are related to structural differences. Compared to the structure of the HEN example in the literature, the structure of the examples presented in Section 4 differs. In the literature example, the respective HEN consists of more streams (seven) while the number of process-to-process HEX (four) is similar to (some of) the presented examples in Section 4. Consequently, it was observed that the distribution of HEXs is different in the literature example. For example, in the literature example, only one process stream is connected to two process-to-process HEXs while all other streams are connected to one process-to-process HEX, only. Additionally, the literature example demands only one control variable. Therefore, it can be assumed that the control structure of the literature example is simpler as the HEN itself is less interconnected compared to the HEN examples in this study. A preliminary conclusion drawn from this is that with increasing structural complexity of the HEN of interest, the risk for failure of the two-level formulation increases. Possible assessment criteria for the structural complexity of a HEN are: • Ratio between the number of process streams and the number of process-to-process HEX: The higher this number, the less complex the HEN; and • Number of process-to-process HEX connected to a single process stream: The higher this number, the more complex the HEN.
It this context, it is worth mentioning that it was possible to derive the complete set of critical points by means of the two-level formulation for the structurally less complex design proposals 3 and 4 while it was not possible for the more complex design proposals 1, 2, and MER.

Appendix A.2. A Karush-Kuhn-Tucker (KKT-)Formulation
The main idea of the KKT-formulation is to apply the KKT conditions to the minimization problem of problem (6) (min x,z,d C(x, z, d, θ)) and then solve the remaining single-level optimization problem (maxd i θ ) for each design variable d i ∈ DV individually. Similar to the result of the two-level formulation, the result of the KKT-formulation is the operating point in the uncertainty span at which the maximum value for d i is obtained while the cost function C(x, z, d, θ) is minimized (i.e., a potential critical point) [44]. If one is able to solve the KKT-formulation globally, the solution represents the global solution of problem (6). In contrast to the two-level formulation, the solution process of the KKT-formulation does not involve an iterative procedure. However, a global NLP solver is necessary to solve the resulting problem, which may cause difficulties for large industrial applications. Additionally, as global solvers usually do not provide marginal values, an (local) NLP solver needs to be applied afterwards to obtain marginal values, which are necessary to identify the influencing uncertain parameters (compare [44]).
In order to obtain marginal values, in this study, the KKT-formulation was first solved using BARON and the found solution was used as an initialization for IPOPT in GAMS. However, the calculation of BARON was interrupted after 300 s and the best (local) obtained solution was used for the initialization of IPOPT. In all tested cases, IPOPT returned the same value as BARON and marginal values were provided. The KKT-formulation was solved for the five different design proposals. For design proposals 1, 2, 3, and 4, the sets of critical points which are shown in Table A5 in Appendix B (i.e., the complete sets of critical points) were obtained by solving the KKT-formulation and applying the set covering algorithm as it was reported by Pintarič and Kravanja in [44]. However, the situation was different for the design proposal MER. The set of critical points that was obtained by the KKT-formulation and the set covering algorithm for proposal MER together with the corresponding TAC and flexibility index are shown in Table A3. Table A3. Results of the Karush-Kuhn-Tucker (KKT-)formulation of design proposal MER. Comparing the TAC and the flexibility index shown in Table A3 to the results for the proposal MER shown in Table 6 indicates that the set of critical points (shown in Table A3) is incomplete. A first solution attempt was to investigate the influence of the uncertain parameters that have "0" elements in the set of critical points (e.g., Fcp H2 in the first identified critical point of the set shown in Table A3). These "0" elements result from the calculation of the marginal values of the uncertain parameters, meaning that the corresponding uncertain parameter was discovered to have no influence on at least one of the design variables when solving the KKT formulation (compare [44]). When solving the design problem, the uncertain parameters corresponding to these "0" elements are defined as variables for the respective critical point, i.e., these uncertain parameters can obtain any value within the uncertainty span at the solution of the design problem. As the identified sets were incomplete, it seemed that the influence of one or several uncertain parameters was not correctly detected. This was investigated by extending the set covering algorithm. In addition to the three steps reported in the appendix of [44], a fourth step was added. In this step, each "0" element present in the set of critical points, obtained after the first three steps, was replaced with all possible combinations with respect to the other points in the set and the set covering formulation (AP6 in [44]) is solved.
With this additional step in the set covering algorithm, the set of critical points shown in Table A3 was transformed to the set of critical points that is shown in Table A5 in Appendix B. Eventually, with the additional step in the set covering algorithm, the complete set of critical points was found, and the feasible design for the proposal MER was achieved (see Table 6).
Obviously, the KKT-formulation is easier to solve for less complex HENs as the different HEX units are less dependent on each other. This is reflected in the achieved results as the extended set covering algorithm was only necessary to identify the complete set of critical points for proposal MER. In all other cases, the extra step in the extended set covering algorithm was not necessary. It was further investigated if increasing the maximum calculation time of BARON would result in different solutions for proposal MER. However, even an increase to 1000 s did not have any impact on the results. Actually, for several design variables, the best (local) solution was found before BARON started iterating, i.e., during the local search before the iteration scheme of BARON is executed. This implies that different starting points could have an impact. However, the influence of different starting points was not studied.