Simultaneous integration of heat pumps and different thermal energy storages into a tightened multi-period MILP HENS superstructure formulation for industrial applications

Design optimization of industrial installations with help of mathematical programs offers high potential for energy savings and cost reduction, which steadily gains importance. Heat exchange network synthesis (HENS) is one of the commonly used and promising methods. This paper deals with simultaneous integration of heat pumps and two different storage types into mixed integer linear programming (MILP) HENS by extending a multi-period superstructure formulation proposed by Beck and Hofmann (2018c). Further measures for tightening the solution space for the coefficient of performance (COP) are introduced. This improves the accuracy of the heat pump (HP) approximation, reduces computational effort and prevents solutions with uneconomical low COP. The obtained approach allows to design cost efficient heat recovery systems with storages (ST) and HP for the improvement of the energy recovery capabilities. A constructed test case was used to analyze the performance of the method and to show its possible potential. © 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
In the last decades extensive research in the area of waste heat recovery has been carried out and it still gains scientific popularity. This is caused by the complexity of the topic and the considerable potential for energy savings which provides economic benefits with short payback times ( Eichhammer and Rohde, 2016 ). The problem gets even more complex by extending classical HENS by additional installations like HP or ST. Simplified approaches are needed to find suitable solutions with acceptable computational effort. The complexity is caused by the nature of the mathematical formulation of HENS. It is proven to be an NP-hard problem, which means there is no possibility for the existence of a polynomial exact solution algorithm for the problem ( Furman and Sahinidis, 2001 ). The advantage of linear programming approaches is that they always deliver global optimal solutions without the need for initial feasible solutions. Even for small problems suitable initial solutions are often hard to find ( Escobar and Trierweiler, 2013 ). As proposed by Nemet et al. (2019) for basic HENS without HP, solutions obtained by MILP approaches could also be used as near optimum starting solutions for a mixed integer nonlinear programming (MINLP) optimization. As already stated by Ciric and Floudas (1991) , the simultaneous optimization with a minimum cost target is superior to sequential performed targeting, because every decomposition adds an element of uncertainty to the solution. Other researchers already approached the problem of integrating HP into heat exchange networks. For example Stampfli et al. (2019) used the principle of heat recovery loops to integrate ST and HP into non-continuous processes by utilizing insight based and nonlinear programming techniques sequential with the target of optimal HP operation. Miah et al. (2015) developed an methodological framework for heat integration which consists of several analysis steps that divide the problem into different zones to decide where HP have to be located. Becker and Maréchal (2012) utilize heat cascade formulations and calculate the optimal solution by setting up an multi objective optimization algorithm that solves the thermodynamic calculations and the energy integration serial as sub problems. In contrast, this work is based on the linearization of the widely used superstructure formulation initially introduced by Yee and Grossmann (1990) , which was proposed by Beck and Hofmann (2018c) .
Differing from the sequential approach for integration of multiple thermal ST used by Beck and Hofmann (2019) , here an simultaneous optimization that considers different thermal ST and HP in one step was developed. The linear integration of the HP and ST allows to find a global optimal solution for the given parameters without the need for initial feasible solutions. All adaptions that have been made compared to previous superstructure approaches are shown in the following Chapter 2 .

Proposed extension of the HENS superstructure
A previous publication Prendl and Hofmann (2020) extended the linearized superstructure proposed by Beck and Hofmann (2018c) to consider the integration of HP. In continuation of this work a further extension has been made by introducing additional storage possibilities to enhance the capability of shifting energy on different temperature levels over time periods. A real vapor compression HP was modeled and simulated to obtain a realistic HP characteristic as explained in Chapter 4 . Also the solution space for the COP has been tightened to improve the accuracy of the approximation and to reduce the computational effort. As can be seen in Fig. 1 , possible installations are a one tank ST (1T ST), a two tank ST (2T ST), heat pumps between streams and the 2T ST, HEX between streams, HEX between streams and ST as well as hot and cold utilities. The 1T ST has a preset fixed mass and variable temperature. Because the initial temperature of the 1T ST is also an optimization parameter the storage gets integrated on the economically optimal temperature interval. The 2T ST operates at two preset temperature levels. The size of the ST gets optimized simultaneously with the solution of the superstructure. For both ST the heat transfer fluid is used as storage medium so no HEX inside the ST are needed. The simultaneous integration of ST with fixed mass and variable temperature as described by Beck and Hofmann (2019) and ST with fixed temperatures and variable mass allows to utilize the benefits of both types simultaneously ( Walmsley et al., 2014 ). The course of the charging states of the different thermal energy storages over the operation periods resulting from the solving of the optimization problem allows to gain better insights into the system behavior. HEX, HP and ST are possible in every stage, while utilities are only possible before the first and after the last stage. The utilities provide the energy needed to fulfill the temperature boundaries of the streams if it is not possible to satisfy them within the stages.

Mathematical model
The objective for the extended stagewise superstructure shown in Fig. 1 is a minimization of the total annual costs (TAC) to provide the needed energy to change the temperatures of the streams to match the given process parameters. The cost function ( Eq. (1) ) consists of three different types of costs:

•
Step fixed investment costs for the installation of streamstream HEX, hot-and cold utilities (Hu,Cu), Storages (ST), streamstorage HEX and heat pumps (HP).
• Variable investment costs for the surface area and thus the size of all HEX as well as variable costs for the size of the two tank storage.
• Energy costs for the external energy demand of the obtained solution which consists of the costs for the hot and cold utilities and the costs for the electrical energy consumed by the HP.

Energy balances
In the following ( Eq. (2) ), the energy balances for all process streams are given with the extensions necessary for considering the heat flows to the storages and heat pumps.
T in and T out , the inlet and outlet temperatures of the streams, where i are the hot streams and j are the cold streams, with their heat capacities cp ip and cp jp as well as their massflows ˙ m ip and ˙ m jp are the physical requirements that have to be fulfilled for a given problem. The balances for every stage k in every time period p which are given in the following Eq. (3) show that isothermal mixing is assumed after every stage. If more than one installation occurs at one one position i jkp, the stream is split up and they are arranged in parallel configuration.
For simplification of the formulation, utilities are only allowed directly before the outlet of streams as in the basic superstructure proposed by Yee and Grossmann (1990) . The utility heat loads in Eq. (4) are calculated as the heat flows that are needed to change the temperatures of the streams to the demanded output temperatures required by the given processes. In the proposed superstructure, utilities are modeled as streams with fixed input and output temperatures that are provided in the needed quantities. They have corresponding HEX to interact with the hotstreams and coldstreams.
For the 1T ST, the accumulated energy content is calculated through the temperature change of the storage medium ( Eq. (5) ), where a perfect mixing inside the storage tank is assumed. The temperature differences for the heat transfers ( Eq. (11) ) are calculated with the temperatures at the start and end of the time periods according to the formulation by Beck and Hofmann (2018a) . For taking the periodic operation into account, the temperature at the start of the first time period has to be the same as the temperature at the end of the last time period and is defined as T shi f t . To consider the physical constraint of thermal stability of the storage medium, a maximum storage temperature T max 1 T is defined. The requirement to keep the program linear prevents a simultaneous optimization of the temperature and sizing of the one tank storage. However, the resulting temperature curve of the storage can be used to gain information of the system and to adapt the storage size or material.
The 2T ST operates at fixed temperatures and the masses in the two tanks change. The energy balance of the ST ( Eq. (6) ) is an adaption of the storage model by Beck and Hofmann (2018a) . The constant temperatures are one of the enablers for the linear integration of heat pumps. This storage can be charged or discharged either directly over HEX from the streams or over heat pumps.
The variable CH 2 T represents the charging state of the ST and has the physical limitation to be between zero (empty) and one (full charged). Similar to the 1T ST, CH shi f t gets used to ensure that the ST has the same state of charge at the beginning and end of the cycle. The storage size S st 2 T gets calculated by multiplying the available storage mass m 2 T with the maximum charging difference that occurs during one cycle. It is assumed as a linear component of the cost function ( Eq. (1) ).
The heat pumps are formulated as power to heat device without losses. In Eq. (7) it can be seen that HP that charge the ST can only occur between hot streams and ST. HP that discharge the ST are only possible between cold stream and ST.

Additional constraints
The constraints given in the following Eq. (8) make sure that the binary variables for the existence of installations are only non zero when the corresponding heat flows lie within their boundaries and thus establish the needed connection between them. For the stream -stream HEX and the stream -storage HEX, the minimum heat flow ˙ Q min has the aim of tightening the solution space for the HEX area. Thus the accuracy of the linear approximations can be improved and the computational time can be reduced ( Beck and Hofmann, 2018b ).
The constraints for the temperature differences for all heat exchanges in the system ( Eq. (11) ) are set up by using BIG-M formulations, where is a sufficient large number to deactivate the constraint if no installation exists at a given position. Additional, a lower boundary for all temperature differences is set as a tightening measure ( Eq. (9) ).
T ≥ T min (9) With decreasing logarithmic mean temperature difference (LMTD) values, the heat exchange areas and thus the variable costs increase, which causes that the optimization wants to maximize the temperature differences in HEX. Thus it is possible to formulate the constraints for the LMTD as given in Eq. (10) . CLMT D are the piecewise linear approximated solution spaces for the logarithmic mean temperature difference as introduced by Beck and Hofmann (2018c) . This approach benefits from the strict convexity of the LMTD that is proven in Mistry and Misener (2016) .

LM T D i jkp ≤ CLM T D i jkp
As physical constraints Eq. (12) are given which prevent that the heat flows or the HEX areas get negative.
The lower boundary for the HEX area is necessary because in the constraints for the heat exchange areas ( Eq. (13) ), the BIG-M coefficient A has to be big enough to make the right side of the equation smaller than zero to deactivate the constraint ( Beck and Hofmann, 2018c ). The boundary for the heat flows enforces the monotonic temperature decrease over the stages.
The decision whether a storage is introduced or not is derived from the existence of installations that connect streams to them as given in the following Eq. (14) .
For the proposed formulation it is assumed, that if a HEX, HP or utility exits at one position in more than one time step p, the biggest heat exchange area gets installed and at the other time steps bypasses balance out the differences ( Eq. (15) ).
The following Eq. (16) formulate the tightening of the solution space for the COP which is described in Chapter 4 and Eq. (17) are the HP formulations derived from Eqs. (8) and (13)

Linearized HP model
Due to the nonlinearities of the process-specific thermodynamic relations as well as material properties, the heat pump characteristics between Pel, ˙ Q and T are of nonlinear nature. As stated in Prendl and Hofmann (2020) , it is necessary to linearize this characteristics for a linear integration of heatpumps into the superstructure. The approach explained in the following is based on the linearization of the nonlinear relationship between COP and T . While in the previous publication ( Prendl and Hofmann, 2020 ), a typical heat pump characteristic was assumed, now a real vapor compression HP was modeled and simulated with a given refrigerant and given temperature ranges to obtain the corresponding HP characteristic.
The T-h-diagram in Fig. 2 shows the thermodynamic cycle of the refrigerant (yellow) in a heat pump. This cycle consists of compression of the superheated vapor (1-2), condensation and subcooling (2-3), expansion (3-4) and the vaporization (4-1). In the HP model these single processes were modeled according to their corresponding thermodynamic relations and connected to a closed thermodynamic cycle. As stated in Eq. (7) any external losses were neglected. The refrigerant R1234 ZE was chosen in the model. With the assumption of exact saturation of vapor and condensation a characteristic for the relationship between COP and temperature difference was obtained. In Fig. 3 this characteristic is shown  alongside with its linear approximation given in Eq. (18) . The maximum deviation between the calculated COP values and the approximation for this case lies at 7.5%.
With this relation the electrical power consumption can be expressed as a function of the supplied heat flow (in this case Q hpst ) and the temperature lift ( Eq. (19) ).
Due to the similarity of the nonlinear equation of Pel to the calculation of the heat exchanger area, a similar linearization approach as Beck and Hofmann (2018c) proposed was used. The nonlinear and nonconvex solution space for Pel gets split up into three regions along the T hp axis according to Eq. (20) with the index x . In this regions Pel as formulated in Eq. (19) gets approximated with linear polynomials as expressed in Eq. (21) with least squares procedure. The maximum of the obtained planes ( Eq. (22) ), which are shown in Fig. 4 , is then used as piecewise linear and convex Constraint CPel in Eq. (17) with means of BIG-M formulations.
Fixed boundaries for parameters of the HP, a minimum and maximum power consumption ( Pel min = 400 kw, Pel max = 20 0 0 kw ) as well as boundaries for the temperature lift T hp min = 20 K and T hp max = 50 K were used. Also the HP approach temperature gets set to fixed value ( T hp ap = 5 K). As further tightening measure, a lower boundary for the COP is chosen ( COP min = 3 ), which cuts of the solution space as displayed in Fig. 5 . The solution space above the red surface gets discarded which reduces the possible solutions and thus the computational effort considerably. The limitation of the range for the COP is justified by the fact, that heat pumps are often only considered as an economically reasonable option over an certain COP. Also the accuracy of the approximation improves and the computational times sink if the approximation domain gets tightened.
In Fig. 6 it is visible, that the percentual deviation of the approximation lies well under 10% in broad areas of the solution space. It can also be seen that the tightening cuts off areas with higher deviations, which are caused by the type and shape of the approximation surface. The lower boundary of the power consumption of the heat pump cuts of areas where even small absolute deviations cause high relative deviations where the calculated power consumption approaches zero.

Test case
For the evaluation of the proposed method, a previous used test case from Prendl and Hofmann (2020) has been extended. The assumed process has a cycle time of four hours and is split into four operational periods of one hour each. It consists of three hot and three cold process streams which have varying mass flows in the different time slices. The annual operation time of the process is assumed as 8600 h. The extended superstructure model was set up with two stages with the stream data and cost coefficients given in Table 1 . The two storages operate with thermo-oil as storage medium, where both oils have a heat transfer coefficient of h oil = 0 . 5 kW m −2 K −1 , the oil for the one tank ST has a heat capacity of cp oil 1 T = 1 . 5 kJk g −1 K −1 and the oil for the two tank ST has a heat capacity of cp oil 2 T = 2 kJk g −1 K −1 . The fixed size of the one tank ST is 10 0 0 0 0 kg and the two tank ST operates at 70 • C and 100 • C. The parameters and boundaries of the chosen HP are given in Chapter 4 except the assumed heat transfer coefficient of h hp = 5 kW m −2 K −1 . Gurobi 8.1.0 was used as solver for the MILP. The test case was chosen rather small to keep the results more traceable and the computational effort low. It has to be mentioned that even in the linear cost function, small changes of the coefficients can have huge impacts on the solution of the system. For example a modification of the variable HEX costs has an effect on all installation options and the changes in the results are hard to trace because of the complexity of the problem. Some parameters have been varied to test the behavior of the optimization for plausibility. When the electrical power costs increase, the system tends to chose HP that operate at points with higher COP, until the costs surpass a critical point where no more HP get chosen for the system. The same behavior is obtained for sinking utility costs. Also when the costs for the two tank ST increase, at some point the ST and also the HP are not chosen anymore because of the connection between them. These reactions are expected from the optimization, and are only traceable for small problems.

Results
The test case described in Chapter 5 was optimized in two different setups with the same parameters. First the optimization was done without the option for ST or HP to be chosen in the HEN to be able to compare the new method to traditional HENS. Then the proposed integration of HP and ST was applied to the problem.

Test case optimized without HP and ST
The optimized HEN without options for HP or ST is shown in Fig. 7 and has total annual costs of T AC = 3 , 132 , 700 € y −1 . It consists of seven stream -stream HEX and six utility HEX and the heat flows of the installations for each time period are given in Table 2 . The total utility energy demand adds up 25.413 GWh y −1 .

Test case optimized with integrated HP and ST
The obtained extended HEN after the application of the proposed method is shown in Fig. 8 and has total annual costs of T AC = 1 , 214 , 400 € y −1 . The network consists of six streamstream HEX, five stream -ST HEX, four utility HEX, two HP, one 1T ST and one 2T ST and its calculated heat flows and electrical power demands are given in Table 3 . For the extended HEN the total utility energy demand is 5.7873 GWh y −1 and the electrical energy demands for powering the heat pumps is 4.7904 GWh y −1 which adds up to a total external energy demand of 10.578 GWh y −1 . The optimized storage size for the 2T ST is m St = 247684 kg thermo-oil with a storage capacity of 4.1280 MWh. Its charging state displayed over the cycle time is given in Fig. 9 . The 1T ST operates between 132 • C and 161 • C and its temperature profile over the cycle time is given in Fig. 10 .

Comparison
The TAC of the HEN obtained with the proposed extended approach are 61.2% lower than the TAC for the network optimized without HP and ST. The total external energy demand is 58.4% lower compared to the total external energy demand of the basic HEN which can also be seen from the values in Tables 2 and 3 . Much less energy is needed because the storages allow to shift the energy over the time periods as visible in Figs. 9 and 10 . The 2T ST gets charged in the first and second period and discharged at the third and fourth period. Similar the 1T ST gets charged during the first time period and discharged in the third and fourth.

Conclusion
A multi period MILP HENS approach that allows the integration of HP has been extended for the option to integrate ST with variable mass and fixed temperatures as well as ST with fixed mass and variable temperatures. In order to get more realistic HP characteristic curves, a simulation of a vapor compression HP with a chosen refrigerant was deployed. Further tightening of the COP solution space is used to improve the accuracy of the approximation and to reduce the computational effort. An previously used test case has been adapted to demonstrate the proposed method and to verify its behavior. The HEN resulting from the optimization with the new approach was able to reduced the TAC by 61.2% and the total external energy demand by 58.4% compared to the classic HEN obtained without HP and ST options. It was found that the computational times of the optimization are very sensitive to changes of physical parameters or cost coefficients. This sensitivity can also be seen in the very different network solutions that result when parameters are varied, which is caused by the nature of mixed integer programming. For comparisons with real applications, additional factors have to be considered. Costs like piping or instrumentation of installations are not included and the costs coefficients are highly dependent on specific geographical and eco-nomical factors. Although, the high reduction of cost and energy demand of the used test case is not directly comparable with real applications, even after inclusion of additional costs considerable savings can be expected trough the utilization of the developed mathematical approach.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.