Dynamic optimization of control setpoints for an integrated heating and cooling system with thermal energy storages

Energy systems for buildings and neighborhoods are expected to become more complex and flexible. Advanced control strategies are required to exploit the full potential of this flexibility and are especially important for systems with storages. In this study, the control of an integrated heating and cooling system for a building complex in Oslo, Norway, was analyzed. Focus was on the control setpoints for the main heat pumps, which had a total heating capacity of about 1 MW and were connected to thermal storage tanks. Previously developed simulation models of the system and its main components were made suitable for dynamic optimization with long time horizons. JModelica.org was used to find optimal control trajectories for the system with two different objectives. The first objective was to reduce the electricity use of the system and the second objective was to reduce the electricity costs of the system. The results showed that the electricity use of the system could be reduced by about 5% for the analyzed year. The electricity costs could be reduced by about 5e11% for the three analyzed winter months, depending on the variability of the electricity price and the size of the storage tanks. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Heating and cooling demands of buildings account for a large share of the world's energy use [1,2]. The development of new buildings and integrated energy systems aims at reducing the environmental impact of these energy demands. Such future systems are expected to be more complex and flexible due to the inclusion of fluctuating energy sources [3] and thermal energy storage (TES) [4]. Simulation and optimization are key methods for the analysis of these complex systems and their operation [5,6]. In particular, optimized control is essential to unlock the full potential of TES [7,8]. However, TES increases the optimization problem complexity, especially when short-and long-term thermal storage are combined [9,10].
Several studies on the optimization of TES operation can be found. Liu et al. optimized the charging of a hot and a cold storage tank with a dual-mode transcritical CO 2 heat pump using Dymola and a genetic algorithm (GA) and report energy savings of almost 20% [11]. Kamal et al. optimized TES control in a large office building using EnergyPlus and a GA and report cost savings of 10e17% [12]. Urbanucci et al. optimized the component sizes and the operation of a cogeneration system with TES by coupling a mixed integer linear programming formulation with a GA [13]. Li et al. optimized the thermal management for industrial waste heat recovery with phase change material TES using a biogeographybased optimization algorithm and report daily fuel savings of 6.9% [14]. All of these studies used derivative-free optimization algorithms, which are not the best choice for the optimization of TES system operation [15]. On the contrary, Knudsen et al. optimized TES operation for surplus-heat exchange in an industry cluster using a gradient-based optimization algorithm and report an increase in surplus-heat utilization from 77% to 98% [16].
In a previous study, an integrated heating and cooling system (IHCS) with TES for a small neighborhood in Oslo, Norway has been analyzed by means of dynamic simulations with Modelica models [17]. It was shown that the control setpoints influenced the system performance. Therefore, the aim of this study was to optimize these setpoints with two main objectives. The first objective was to analyze the potential reduction of the system's electricity use with the currently installed components. The second objective was to analyze if the installation of larger storage tanks could lead to a reduction of the system's electricity costs. To this end, dynamic optimizations with JModelica.org were performed. JModelica.org is a Modelica-based open-source optimization platform, which essentially uses a gradient-based algorithm to solve a nonlinear programming (NLP) problem [18]. The optimization models in this study contained both short-and long-term TES and a heat pump model. This led to a challenging optimization problem due to the different time-scales of the storages and the part-load operation of the heat pump. In addition, relatively long time horizons were necessary to avoid suboptimal usage of the long-term storage during optimization. This led to relatively large optimization problem sizes. The resulting NLP sizes in this study were larger than the reported sizes in other studies using JModelica.org, e.g. Refs. [16,19e24].
The remainder of this paper is structured as follows: a short description of the IHCS and brief results from the previous study are given in Section 2. In Section 3, the optimization approach is described in detail. Results from the optimizations are presented in Section 4, followed by a discussion in Section 5 and concluding remarks in Section 6.

System description
The IHCS was located in the Norwegian capital Oslo and supplied thermal energy for space heating (SH), domestic hot water (DHW) heating, snow melting (SM), space cooling (SC), and product cooling to several different building types. Snow melting was applied to the walkways between the buildings and product cooling was supplied to a food court. A simplified schematic of the IHCS is shown in Fig. 1 The main components of the IHCS were heat pumps (HPs), heat exchangers, solar collectors, storage tanks, ice thermal energy storage (ITES), and borehole thermal energy storage (BTES). The ITES was used to reduce space cooling peak demands during summer days and was charged during the night. The BTES was used as heat source during heating season and as heat sink during cooling season. The IHCS was designed to supply heat at 50e55 C and was also connected to the city's district heating (DH) system. This DH system supplied heat at a temperature of 85e120 C and was used as backup (in case of system failure), for peak demand coverage, and as temperature lift for DHW heating, which had a supply temperature setpoint of 70 C. A more detailed description of the IHCS can be found in Refs. [17,25].
The IHCS was equipped with a control and monitoring platform.

Results from the previous study
The IHCS was analyzed in a previous study by means of dynamic simulations [17]. Modelica models were developed for the main components with focus on a good trade-off between accuracy and simulation time. The system model was calibrated with one-year measurement data and the results showed good agreement between simulated and measured values for electricity use and DH import.
The simulation results showed that more heat was extracted from the BTES during winter than was injected during summer. This can be seen in Fig. 3, which shows the daily heat balances for the BTES and the solar collectors. This imbalance was shown to reduce the system's long-term performance.
A sensitivity study was performed and the results showed that the supply temperature setpoints for heating, T supply_heat , and space cooling, T supply_cold , were important for the system performance. These setpoints were used for the control of the main heat pumps as well as the BTES circulation pumps, which were responsible for 62% and 5% of the total electricity use, respectively. The IHCS had two operation modes: heating mode and cooling mode. T supply_heat was set to 55 C during heating mode and 51 C during cooling mode. T supply_cold was set to 6 C during both operation modes. This control approach is called business as usual (BAU) throughout this paper.

Methodology
The IHCS was modeled in Modelica and simulated with Dymola as described in the previous section. The aim of this study was to optimize the operation of the IHCS, in particular the setpoint temperatures for heating and cooling supply, T supply_heat and T supply_cold , respectively. To this end, dynamic optimizations were performed with JModelica.org. The resulting setpoints from the optimizations were implemented into the system simulations in Dymola. The workflow of the entire analysis is shown in Fig. 4.
All elements of the optimization procedure (Part 2 in Fig. 4) are explained in detail in the next subsections. First, the optimization platform JModelica.org is described in Section 3.1. The development of seasonal models, which were suitable for numerical optimization, is described in Section 3.2. Finally, in Section 3.3, the control variables, constraints, and objective functions of the optimization problems are explained.

Optimization platform
JModelica.org is an open-source platform for simulation and optimization of complex dynamic systems [18]. JModelica.org is based on the modeling language Modelica and the Functional Mock-up Interface standard which enables coupling to different software packages. It uses the language extension Optimica, which enables high-level formulation of optimization problems [26]. All the optimizations in this work were performed with JModelica.org version 2.2 via 64-bit Python scripting. The main steps of the optimization procedure are described in this section and are shown in Fig. 5.
Step 1: An initial simulation was required to obtain variable trajectory data for initialization and scaling of the NLP variables in Step 5, see Fig. 5. To this end, the Modelica model for initialization was compiled into a Functional Mock-Up Unit and simulated using the CVode solver from the SUNDIALS suite [27], which is included in JModelica.org.
Step 2: The Modelica model for optimization and the problem formulation (Optimica code) were compiled and transferred to the CasADi interface of JModelica.org. CasADi was used for the computation of derivatives using algorithmic differentiation [28].
Step 3: Symbolic elimination based on block-triangular ordering was applied to reduce the number of algebraic variables as explained in Ref. [29]. This step was found to be crucial for successful converge as it significantly reduced the size of the resulting NLP.
Step 4: Code for orthogonal collocation on finite elements [30] is included in JModelica.org and was used to transform the infinitedimensional dynamic optimization problem into a finitedimensional NLP. The number of collocation elements and the number of collocation points in each element has a strong influence on the size of the resulting NLP. Initial testing with a prediction horizon of one-week was performed to compare the resulting setpoint temperatures with two different collocation configurations. The fine discretization had an element length of 15 min and  two collocation points per element. The coarse discretization had an element length of 30 min and one collocation point per element. Applying the coarse discretization led to a reduction of the NLP variables from 598•10 3 to 178•10 3 and a reduction of the solution time from 103.2 s to 30.4 s for the one-week prediction horizon. However, the average absolute difference between the optimized setpoints for the fine and coarse discretization was less than 0.1 K, which was regarded as insignificant. Therefore, the coarse collocation discretization was used for all the optimization cases in this study.
Step 5: Variable trajectory data obtained during the initial simulation (Step 1) was used for automatic initialization and scaling of the NLP variables.
Step 6: The resulting NLP was solved using version 3.12.4 of the primal-dual interior-point solver IPOPT [31] with linear solver MA57 from HSL [32]. All optimizations were performed with an Intel® Core™ i7-6700 K processor (4 GHz) and 64 GB RAM.

Optimization models
Initial testing showed that the complexity of the full system model developed in Ref. [17] impeded its applicability for dynamic optimization. To enable optimization of the system model over relevant prediction horizons, seasonal models were developed as a means of reducing model complexity. This section describes the reduction of the full system model and its decomposition into seasonal models.

Reduction of the full system model
The full system model could not be used for dynamic optimization due to the large number of components and their interconnections as mentioned above. Therefore, certain parts of the system had to be removed to reduce the complexity and the size of the resulting NLP. The full system model and the removed parts (covered with gray) are shown in Fig. 6.
It can be seen in Fig. 6 that the solar collector loop, the DHW heating substation, the product cooling and ITES charging loop, and the recovery loop were removed. The solar collector loop was removed because it played a minor role for the system performance, see Fig. 3. The DHW heating substation, the product cooling, and the ITES charging loop were removed because the recovered heat from HP3 was similar to the supplied heat in the DHW heating substation (620 MWh and 682 MWh, respectively). Removing these parts therefore caused insignificant mismatch in the total heat balance. The simulated electricity use of the removed parts accounted for 18% of the total electricity use in the previous study, which showed that the key components of the system were kept. For clarity, the reduced system model is shown in Fig. 7.

Modifications of the component models
The component models described in Ref. [17] were developed for stable and fast dynamic simulations. However, due to the    The substation model for simulation, see Fig. 6, received a demanded heat flow rate as input signal, which was based on the corresponding demand type. This input signal was used to control the mass flow rate of the circulation pumps within the substation model. Initial testing showed that the implemented PI-controller model led to convergence issues. Therefore, a different approach was chosen for the optimization model. The mass flow rate of the circulation pumps was used as input signal, see the yellow boxes in Fig. 7, and the required heat flow rate was formulated as a constraint in the optimization problem, see Section 3.3.
Similar to the substation model, a new approach was also chosen for the HP model. The HP model for simulation received one of the outlet temperatures on the secondary side as input signal (not shown in Fig. 6). The HP model for optimization received the heat pump power (P HP ) as input signal, see the yellow box in Fig. 7. Initial testing showed that this modification increased the convergence rate significantly.
The calculation of the Lorentz temperature in the heat pump model in Ref. [17] impeded successful convergence and was therefore approximated in the optimization model as: The difference in Lorentz temperature due to this modification was less than 0.1 K for all relevant operating conditions, which was regarded as insignificant.
The numerical discretization of the BTES and the storage tanks had strong influence on the number of NLP variables. A one-week test optimization was performed to compare the resulting setpoint temperatures with high and low discretization values. The horizontal and vertical discretization of the BTES was set to 30 and 4 for the high discretization case and 10 and 2 for the low discretization case, respectively. The discretization of the heating and cooling tanks was set to 15 and 5 for the high discretization case and 5 and 2 for the low discretization case, respectively. Reducing the discretization led to a reduction of the NLP variables from 178•10 3 to 53•10 3 and a reduction of the solution time from 30.4 s to 2.8 s. However, the average absolute difference between the optimized setpoints for the high and low discretization case was less than 0.1 K, which was regarded as insignificant. Therefore, the low discretization values were used for all the optimizations in this study.

Seasonal models
Some parts of the IHCS were only used during certain periods of the year, because the heating and cooling demands varied from season to season, see Fig. 2. Optimizing unused parts would unnecessarily increase the optimization problem size. Therefore, three seasonal models were created based on the reduced system model shown in Fig. 7 and the unused parts of each model were removed. For the winter model, the BTES charging heat exchanger and respective circulation pumps could be removed because no heat was sent to the BTES during winter, see Fig. 3. For the transition model used for spring and fall season, the snow melting substation could be removed since there was no snow melting demand, see Fig. 2. For the summer model, the snow melting substation and the BTES discharging pump could be removed because there was no   Table 1 together with the unused parts. Two versions of each seasonal model were required: one for the initial simulation and one for the optimization, see Fig. 5. In the initialization models, the component models for simulation were used. The component models for optimization, explained in Section 3.2.2, were used during the optimizations.

Optimal control problem formulation
The seasonal models described in the previous section were used to find optimal heating and cooling supply temperature setpoints for simulations with the full system model as shown in Fig. 4. The optimization problems for the different models were formulated as continuous-time optimal control problems. The control variables, constraints, and objective functions of the optimization problems are explained in the following subsections.

Control variables
The control variables in the optimal-control problems were the heat pump power, P HP , and the mass flow rates for the circulation pumps. These are marked yellow in Fig. 7 and are written as a vector: i2fSH; SM; SC; BTES_Evap; Cond_BTESg (2) The temperatures T supply_heat and T supply_cold were not included in the vector u(t). This was due to the fact that the optimization models did not contain PI-controllers, as explained in Section 3.2.2, and thus could not receive a setpoint temperature. The temperatures T supply_heat and T supply_cold depended on the control variables and were calculated during the optimizations. The resulting values were then used as input for the new simulations (Part 3 in Fig. 4).

Constraints
Lower and upper bounds were defined for the control variables based on their operational limits, yielding the following linear inequality constraints: The maximum mass flow rates in Equation (4) were set based on manufacturer specifications as follows:  the supply temperatures for heating and cooling were constrained by: T supply heat ðtÞ 65 C T supply cold ðtÞ ! À 5 C Constraints were also added to ensure that the correct amount of energy was supplied from the IHCS to the connected buildings. Enforcing this demand satisfaction as an equality constraint led to convergence issues. Therefore, the following upper and lower bounds were defined for the heat flow rates in the substations based on the heating and cooling demands: This formulation improved the numerical performance. The parameter ε was set to 1.005 so that only a small slack in energy supply was allowed.

Definition of the objective function for reduction of electricity use
The calculated electricity use of the system consisted of three parts: the electricity use of the heat pumps, the electricity use of the circulation pumps, and a constant term from auxiliary systems [17]. The first aim of this study was to analyze how much this electricity use could be reduced. Therefore, the following objective function was defined in order to minimize total electricity use: Note that the constant term of the electricity use had no influence on the optimal solution and was therefore removed from the objective function. The year was divided into seasonal periods and each period was optimized separately with the corresponding values for the period's beginning (t start ) and end (t final ). The length of each season and the resulting NLP problem size of the respective optimization are shown in Table 2.
The initial state of the BTES and storage tank models for each season were chosen based on the result of the previous season.

Definition of the objective function for reduction of electricity costs
The second aim of this study was to analyze if the electricity costs of the system could be reduced with improved control setpoints. In Norway, electricity prices are much higher during winter than during summer due to the market based electricity price model and the high amount of electricity used for space heating. The first three months of the year accounted for 44% of the total electricity costs for the simulated year. Therefore, these three months were chosen for the cost-reduction analysis to limit the number of required optimization runs. This way, all the optimizations could be performed with the winter model.
The electricity spot prices for the location of the IHCS for the first three months of the previous four years are shown in Fig. 8. At the time of writing, the exchange rate from Norwegian Krone to Euro is 1 NOK ¼ 0.0982 EUR (XE [34]. Note that the prices in Fig. 8 are market spot prices. Customers also have to pay electricity grid prices and additional fees, which were not considered in this study. It can be seen from Fig. 8 that the electricity price showed relatively little variation in 2015. Therefore, additional price signals were defined with different fluctuations to analyze the influence of the variability of the electricity price, v, on the cost saving potential. The price signals were based on the average price of the first three months of 2015 (239 NOK/MWh) and the original price signal (e Oslo2015 ). Values of 0, 1, 2, and 3 were chosen for v and the price signals were calculated as follows: The four resulting price signals were used for the optimizations and are shown in Fig. 9.
This approach, similar to the one presented in Ref. [35], was chosen instead of using electricity prices from other years to maintain the correlation between the electricity price and the climate conditions. Note that this correlation is not kept for v ¼ 0,  which corresponds to a constant and thus unrealistic electricity price.
The storage tanks of the IHCS were relatively small and only used as buffer to even out the supply temperatures of the heating and cooling loop. Storage tanks are a relatively cheap component, so the installation of larger tanks may be considered as a realistic retrofitting option. To investigate the effect of larger storage tanks on the cost saving potential, three different tank size combinations were chosen: the installed 10 m 3 and 2 m 3 for the heating and cooling tank, respectively, as well as 100 m 3 and 500 m 3 for both tanks.
The four different price signals and the three different tank size combinations led to the twelve optimization cases listed in Table 3.
All the cases listed in Table 3 were optimized separately with the winter model. According to the prediction horizon of three months, t start and t final were set to 0 and 7.776•10 6 in the objective function, Equation 10, respectively. Optimal operation over this prediction horizon would lead to emptied short-term storages at t final , i.e. the average temperature (T avg ) in the hot storage tank would be as low as possible and the average temperature in the cold storage tank would be as high as possible. This would lead to an unfair comparison, especially when different tank sizes were compared. Therefore, the following constraints were added for these twelve optimizations to avoid this effect and thus ensure a fair comparison:

Results
In this section, the results from the simulations and the optimizations are presented. First, the results leading to reduced Table 3 Case IDs of optimizations for electricity cost reduction. electricity use of the system are shown. For these results, the five optimizations listed in Table 2 were performed, which included all the seasonal models and an analyzed period of one year. Afterwards, the results leading to reduced electricity costs of the system are presented. For these, the twelve optimizations listed in Table 3 were performed with the winter model and an analyzed period of three months. Note that perfect prediction of the heating and cooling demands as well as the electricity price was assumed for all the simulations and optimizations.

Reduction of annual electricity use with optimization-based control
A validation was performed to confirm that the energy demand constraints, Equations (7) and (8), were not violated during the optimizations. To this end, a one-year simulation with the reduced system model shown in Fig. 7 was performed. The optimized values for the control variables were used as input and the resulting heat flow rates in the substations' heat exchangers were compared to the demanded heat flow rates. Daily values for demanded and supplied energy are shown in Fig. 10.
It can be seen from Fig. 10 that there was no mismatch between the demanded and the supplied energy. Slight deviations were observed on hourly basis. This was due to the slack formulation in Equation (8) and the fact that the constraints were only enforced at the collocation points and not during the entire width of the collocation element. However, the R 2 -values for all three demand types were above 0.99 on hourly basis, showing that the deviations were insignificant.
The optimized values for heating and cooling supply temperature are shown in Fig. 11 and Fig. 12, respectively. The former setpoints from the previous study, used during Part 1 in Fig. 4, are also shown for comparison.
It can be seen from Fig. 11 that the optimal heating supply temperature was lower than the BAU setpoint for most of the year, except for a few winter days. It can also be seen that the optimal temperature was lower during summer than during winter, so the reduction of the setpoint during cooling mode was a good choice for the BAU case. Fig. 12 shows that the optimal cooling supply temperature varied around the BAU setpoint throughout the year. It was higher during summer than during winter except for a short period in the beginning of the year. This period might be affected by the initialization of the BTES. In general, the optimal cooling supply temperature showed larger variations than the optimal heating supply temperature.
The optimized values for T supply_heat and T supply_cold were implemented into the full system model, where they were used as replacement for the mode-based setpoints of the BAU case (Part 3 in Fig. 4). The resulting energy amounts for the simulated year are shown in Fig. 13.
It can be seen from Fig. 13 that the electricity use for the heat pumps and the circulation pumps decreased by 5 and 14%, respectively, with the optimized setpoints compared to the BAU case. Due to the circulation pump's low share of electricity use, this corresponded to a minor total reduction of electricity use. The amount of heat imported from DH increased by 12% for the simulated year. The amount of heat taken from the long-term storage decreased by 7%.

Reduction of electricity costs during winter with optimizationbased control
Selected result values from the optimizations leading to minimized electricity costs are shown in this section. February 14 th and February 3 rd were days with very different variations in electricity spot price. The price signals for these two days are shown in Figs. 14 and 15, respectively.
It can be seen from Fig. 14   constant on February 14 th . On the contrary, the electricity price varied significantly on February 3 rd as shown in Fig. 15. Peak hours were in the morning and afternoon, which is typical for Norway [36]. Detailed results for the optimal heat pump power and temperature setpoints are presented for these two days for selected cases from Table 3. The results for February 14 th for the cases with the original electricity price and different tank size combinations are shown in Fig. 16. It can be seen from Fig. 16 that the different tank size combinations yielded very similar results for February 14 th . This was expected due to the relatively constant electricity price during that day. The results for February 3 rd for the same cases are shown in Fig. 17.
It can be seen from Fig. 17 that the optimal control trajectories for February 3 rd depended highly on the size of the storage tanks. Larger tanks led to larger variations, due to the possibility to shift electricity use from peak hours (with high prices) to off-peak hours (with low prices) and thus decrease the total electricity costs. Fig. 17 clearly shows that the installed tanks (Case 10e2_e 1 ) were too small to take advantage of the electricity price variations. The heat pump power only varied between 150 kW and 270 kW for this case and the temperatures setpoints were relatively constant as well, expect for two short peaks of T supply_heat . For Case 100-100_e 1 , the heat pump power varied across nearly the entire allowed range from 0 to 300 kW. It was higher during off-peak hours to charge the storage tanks, corresponding to high values for T supply_heat and low values for T supply_cold . On the contrary, the heat pump power was low during peak hours and the energy demands of the buildings were to a large extent covered by discharging the tanks. For Case 500-500_e 1 , this effect was even more pronounced, leading to the largest variations in the optimal values for T supply_heat and T supply_cold .
The results for February 3 rd for the cases with the largest tanks and different variability of the electricity price are shown in Fig. 18.
It can be seen from Fig. 18 that there were large differences between the results with a constant electricity price, Case 500-500_e 0 , and the cases with price variations. Although the costs were optimized for all the cases, the constant price led to a minimization of the total electricity use for Case 500-500_e 0 (i.e. Equations (9) and (10) yielded equal results). The control of the heat pump and the circulation pumps were therefore optimized depending on the energy demands of the buildings. For the other three cases, the electricity use was significantly higher during offpeak hours. The cases with different variability showed very similar results for February 3 rd . The optimal control trajectories became slightly more pronounced for larger values of variability, but only T supply_cold showed significant differences. This showed that even larger tanks would be required to take advantage of the variations during that day. However, other days showed larger differences between these cases.
The optimized setpoints were implemented into the full system model and a simulation for the first three months was performed for all the cases listed in Table 3. The simulated total electricity costs for this period are shown in Fig. 19. The simulated costs with BAU control were included to show the potential savings. All the results are shown relative to the BAU case, because the different price signals led to different costs for the BAU case.
It can be seen from Fig. 19 that all the optimized cases led to lower electricity costs compared to the BAU case. The relative savings were in the range of 5e11%. The relative savings increased with larger variability of the electricity price signal. Larger tanks also led increased relative savings, except for the cases with constant electricity price (e 0 ). Peak hours and off-peak hours were defined to quantify how much of the total electricity use was shifted from the former to the latter. For each day, the minimum and maximum electricity price were used to define the daily price range. Hours with a price in the top 25% of this range were defined as peak hours. Similarly, hours with a price in the bottom 25% were defined as off-peak hours. Fig. 20 shows the electricity use during different pricing hours for the BAU case and the optimized cases.
It can be seen from Fig. 20 that the installed tanks (10/2) were too small to shift electricity use from peak hours to off-peak hours. On the contrary, the percentage of electricity use during off-peak hours was increased for the cases with larger tanks. This load shift was more pronounced for the cases with higher variability of the electricity price.

Discussion
In this section, general matters regarding the applied methodology are discussed first. Afterwards, the reduction of electricity use and electricity costs are discussed in detail.
All the optimizations were performed with reduced system models. Minor components, responsible for 18% of the annual electricity use, were removed from the full system model to avoid convergence issues of the NLP solver. In addition, the year was divided into five seasons and each season was optimized separately. From this point of view, the results can be seen as a lower bound, since the control of the removed parts were not optimized. This means that potential improvements were disregarded and a oneyear optimization of the full system model would be desirable. However, numerical optimization is significantly more challenging than simulation for this type of integrated systems and the complexity of the full system model impeded a one-year optimization. The authors find it worth noting that JModelica.org version 2.0 was used initially, which only supported 32-bit Python. The memory usage of a 32-bit Python process is limited to about 2 GB, which was insufficient for the optimizations in this study and led to frequent memory allocation errors. JModelica.org version 2.2 was released in March 2018 and was the first version to support 64-bit Python. The upgrade to version 2.2 was crucial for this study and the same results could therefore not have been produced before March 2018. The input data for the whole year were used as input in this study. The optimizations were thus performed with perfect prediction. From this point of view, the results can be seen as an upper bound, since perfect prediction is not a realistic scenario. The energy demands of buildings and the electricity price in Norway both depend on ambient conditions. In practice, the uncertainty of the weather forecast thus makes detailed optimizations over a long prediction horizon obsolete. Shorter periods are therefore chosen in practical applications such as Model Predictive Control [37]. An advantage of shorter prediction horizons are that more detailed models can be optimized. A disadvantage is that the use of longterm storages needs special attention. For a short prediction 500-500_e 1 100-100_e 1   horizon, the optimization of a long-term storage is fundamentally difficult. Unsustainable usage may result unless a sufficiently high cost is put in the objective function or constraints are imposed. The implementation of such measures was outside the scope of this study and long prediction horizons were chosen to ensure optimal operation of the BTES. The simulated electricity use of the IHCS was reduced by improving the control setpoints T supply_heat and T supply_cold . The coefficient of performance (COP) of the heat pump model was calculated based on a constant Lorentz efficiency, see Ref. [17]. Thus, the COP depended on the temperature lift of the heat pump, which varied significantly with the optimized setpoint trajectories. Part-load operation was therefore included in the model, but depended only on the temperature levels and not on the heat flow rates. The optimized setpoints led to reduced electricity use of the heat pumps and the circulation pumps as shown in Fig. 13. However, the electricity used by the heat pumps was converted to useful heat and thus partly covered the heating demands of the buildings. Reducing this electricity use thus led to more heat being imported from DH. The increase in DH import depended on the amount of heat taken from the BTES. Too high heat extraction from the BTES can lead to unsustainable operation as shown in the previous study [17]. However, importing heat from DH is more expensive in the short run. An economic analysis is therefore required to find the optimal operation strategy that balances short-and long-term cost considerations. This was outside the scope of this article, because the DH import was not included in the optimization models. This should be added in future studies, but as DH is an essential technology for decarbonization, its use should generally be preferred over electricity use [38]. Charging the BTES with low-grade heat from the DH return line could thus also be an interesting option to investigate in future work.
The electricity costs shown in Fig. 19 were calculated by multiplying the electricity use of the system by the local electricity spot price, the 25% taxes that have to be paid were neglected. However, this is only a part of the actual costs that large customers have to pay in Norway. The electricity grid in Norway is stressed significantly more during the winter than during the rest of the year due to the high use of electricity for space heating. Therefore, the electricity grid prices include additional costs to consider the electricity grid stress. For business customers, this may induce peak-load tariffs and charging for their peak electricity use of each calendar month. This was not taken into account in this study as the measurement data showed that the peak use of the IHCS was almost the same for all the winter months. This cost was therefore assumed fixed and not included in the optimizations. The trend towards more flexible systems and higher incentives for peak load reduction may require including peak-load tariffs in future  optimization studies. Such studies could also influence the sizing of the heat pump and the storages during the design phase by taking flexible operation into account. However, it is important to ensure that the heat pump can handle the large variations in operation.
The relative savings shown in Fig. 19 were obtained by comparing three-month simulations with different tank size combinations and different temperature control setpoints. Larger tanks were shown to lead to reduced electricity costs. However, the difference between the BAU case (10/2 -BAU) and the case with the currently installed tanks and optimized setpoints (10/2 -Optimized) was larger than the difference between the cases with different tank sizes and optimized setpoints (10/2 e Optimized vs. 500/500 -Optimized). This means that the optimized control led to higher relative savings than the installation of larger tanks. However, these savings only included the electricity costs and not the costs for DH import. Since the DH import increased for the cases with the optimized setpoints compared to the BAU case, an economic analysis including the calculation of the total operating costs is required to decide if larger storage tanks should be installed. The costs for the advanced control system should be taken into account in such an analysis since the installation of larger tanks would not lead to savings with the current control strategy.

Conclusions
The simulated performance of an integrated heating and cooling system with thermal energy storages was analyzed in this study. Dynamic optimizations were applied to find optimal control trajectories for operation leading to reduced electricity use and reduced electricity costs. The results showed that the electricity use of the system could be reduced by about 5%. However, this led to increased import of heat from the district heating grid. Possible savings therefore depended on the electricity and district heating prices.
The installation of larger storage tanks was shown to decrease electricity costs when the optimized control setpoints were implemented. However, the savings depended on the variability of the electricity price and could only be achieved with a more advanced control system than the one currently implemented. During the analyzed period, the variability of the electricity price was too low to make the installation of larger storage tanks seem profitable in practice. Higher peak-load tariffs and/or an increased variability of the electricity price might change this conclusion in the near future. Further work should therefore include more detailed optimization models and more advanced cost calculations.