An Integer Linear Programming Model for an Ecovat Buffer

An increase in the number of volatile renewables in the electricity grid enhances the imbalance of supply and demand. One promising candidate to solve this problem is to improve the energy storage. The Ecovat system is a new seasonal thermal energy storage system currently under development. In this paper, an integer linear programming model is developed to describe the behaviour and potential of this system. Furthermore, it is compared with a previously developed model, which is simplifying the behaviour of the Ecovat system much more, but is much less computationally expensive. It is shown that the new approach performs significantly better for several cases. For controlling a real Ecovat system in the future we may incorporate a number of improvements identified by our comparison analysis into the previously developed approach, which may help increase the quality of the obtained results without increasing the computational effort too much.


Introduction
To reduce greenhouse gas emissions, renewable energy sources such as sun and wind are more and more favoured over fossil fuels.Aside from that, the amount of cheaply accessible fossil fuels is limited, thereby increasing the need for sustainable alternatives even more.However, due to the unpredictable nature of renewables it is harder to match supply and demand.One of the solutions to this problem is (thermal) energy storage (e.g., [1,2]).For example, short term thermal and electrical storage can be used to shift excess supply to times of the day with a high demand.Likewise, seasonal thermal storage can be used to shift surplus thermal energy from high-supply periods (summer) to high-demand periods (winter).
In this paper the focus is on (seasonal) thermal storage.A considerable amount of research has been done in the past decades on such thermal storage.In [3] an overview of thermal storage technologies is presented, while [4] focusses specifically on seasonal thermal storage.Hereby thermal energy storage systems can be divided into three classes: sensible, latent and chemical storage.In sensible thermal storage the energy of the storage is increased by increasing the temperature of the storage medium.The most used medium is water, because of its high availability, low cost and high specific heat.Examples of sensible thermal storage using water as the storage medium include water tanks [5,6], solar ponds [7] and aquifers [8].An important phenomenon in such water-based heat storage systems is thermal stratification (e.g., [9,10]); this phenomenon is discussed extensively in e.g., [11].Another common storage medium for sensible storage is rock, for example in rock beds [12,13].One of the advantages of using rock is that it can work for temperatures over 100 • C, while a disadvantage is that a larger volume is needed to store the same amount of energy when compared to water.In latent thermal storage energy is stored due to a phase change in the storage medium, so called phase change materials.Advantages of latent thermal storage include a smaller temperature range due to the phase change of the medium happening at a certain nearly constant temperature and requiring smaller storage volumes, due to a higher energy density, when compared to sensible storage.Reviews on thermal storage using phase change materials are given by e.g., [14,15].In chemical storage energy is stored by reversible chemical reactions or sorption processes.In general these systems have an even higher energy density, however, they are currently more expensive.A review on chemical thermal storage is given in e.g., [16].The Ecovat system falls in the sensible thermal storage class, more specifically it is a water tank system.There are however some differences compared to traditional water tank systems, as is discussed later in this section.To the authors' knowledge, the optimization of the control of a system like the Ecovat has not been researched before aside from a heuristical model for the Ecovat buffer described in [17].
The Ecovat system [18] is a new seasonal thermal storage solution that is currently under development.The most important component of the Ecovat system is the Ecovat buffer, which is a large subterranean buffer filled with ground water (Figure 1).The buffer is divided into a number of segments, N seg , which can each be charged or discharged separately using heat exchangers integrated into the buffer walls.Note that these segments are not physically separated, but specify a certain region within the buffer.The charging is done through a set of devices accompanying the buffer, a resistance heater, photovoltaic thermal (PVT) panels and a number of heat pumps.These devices are used as follows.In summer the heat of the PVT panels can be used to charge the buffer.In addition to thermal energy the PVT panels also generate electrical energy that can be used to power the other devices in the Ecovat system or be sold on the energy market.The resistance heater is exclusively used when the energy price is negative, due to its low efficiency.During those times the resistance heater can be used to charge the buffer with a large amount of energy in a short time.This can be done all year round.Aside from charging the buffer with external energy, the energy quality of the buffer can also be increased by using the heat pumps to transfer energy from a cold to a warm segment within the buffer.Using the heat exchangers in a smart way ensures that the buffer remains thermally stratified, which has been shown to increase the efficiency of the buffer [9,10].Construction of an Ecovat buffer in Uden, the Netherlands [18].This Ecovat buffer is 16 meter deep with a diameter of 11 meter.However, in the future construction of even larger buffers is planned.
The buffer is very well insulated to prevent heat losses to the surroundings.A thermal analysis from an internal report revealed heat losses of about 10% or less, depending on the size of the buffer, over a period of 6 months.The discerning feature of the Ecovat buffer compared to other thermal storages is that there is no water being pumped in or out of the buffer.Charging and discharging of the buffer is done solely through the heat exchangers in the buffer walls, which greatly reduces the amount of water movement in the buffer compared to pumped storage.
To make efficient use of an Ecovat system a method is needed that determines the optimal strategy for the charging and discharging of the buffer given an initial temperature distribution for the buffer segments, as well as predicted weather conditions and energy prices during the considered period.For this, first a proper model of the Ecovat itself is needed.This model mainly has to specify the temperature (energy) distribution within the buffer, which depends in reality on all three spatial dimensions as well as on time.However, due to the complexity of the problem a three-dimensional model of the buffer is computationally very expensive.To simplify the model the buffer is described in fewer dimensions.As such the model proposed in this paper is a one-dimensional model; i.e., each segment is modelled as having a single temperature.This model can form the basis for an extended model in more dimensions if needed in future work.
This paper presents an integer linear programming (ILP) model to describe the optimal charging and discharging of the Ecovat system.Additionally, results obtained from the ILP are compared with results from a heuristical approach described in [17], which runs much faster but gives potentially suboptimal charging strategies.The goal is to use the insight of our model to develop a fast heuristic, that is comparable in quality to the ILP model, but which allows for much faster simulation times compared to the ILP model.Such an approach is especially important when controlling an Ecovat system in real time or incorporating the Ecovat model into a larger simulation, for example a demand side management simulation of a neighbourhood including an Ecovat buffer.To this end the ILP model developed in this work is compared with the previously developed heuristic model to identify whether the results of both models are of comparable quality and which improvements (if any) should be made to the heuristic model.
The structure of this paper is given below: in the following section a detailed description of the ILP model of the Ecovat system and its charging/discharging is presented as well as a short description of the heuristic model described in [17].The input profiles and parameters needed for the ILP model are discussed in Section 3. The results of the ILP model and heuristic model are compared and discussed in Section 4. Finally, in Section 5 conclusions based upon the current work are presented and plans for future work are outlined.

Model
The developed ILP model of the Ecovat system considers the relevant components of an Ecovat system as illustrated in Figure 2.For the case shown, the number of segments is N seg = 5.The water inside the Ecovat buffer will be stratified, i.e., different layers are formed due to a difference in density, in this case caused by differences in temperature.The layers are ordered from least dense (highest temperature) at the top of the buffer to most dense (lowest temperature) at the bottom of the buffer.For the Ecovat buffer this means that a given segment will always have a temperature equal to or greater than the segment below it.The input for the model consists of predictions for the energy prices, weather data (ambient temperature and solar radiation) and heat demand for the optimization period.The model returns the optimal charging strategy, which is determined by the decision variables of the model, as well as the temperature evolution inside the buffer during the considered time horizon.The set of segments is defined as S = {1, ...., N seg }, which are numbered from top to bottom as shown in Figure 2. The time horizon is discretized into a set of time intervals and defined as T = {1, 2, ...., N time + 1}, where N time is the number of time intervals in the time horizon.The temperature of segment s at time t is given by T t,s , where the simplifying assumption has been made that an entire segment has a single temperature value at any given time step.Finally, the set of devices that charge/discharge the Ecovat buffer, D, is defined as D = {pvt, aw, ww1, ww2, res, dem}, where the labels designate the different devices in the system.The labels correspond to the PVT panels, the air/water heat pump, the high temperature and low temperature water/water heat pumps, the resistance heater and the heat demand respectively.We assume that each device is in the same state during a complete time interval.

Decision Variables
All the considered devices are on/off devices, this means all the decision variables can be modelled as binary variables that describe whether a device is connected to segment s at time t, i.e., the decision variables are written of the form: x d t,s and represent the state of device d during time interval t with respect to segment s.A decision variable with a value of 1 means the given device is turned on and is connected to segment s at time t, while a decision variable with a value of 0 means the given device is not connected to segment s at time t.If at time t the variable x d t,s = 0 for all s ∈ S then the device is turned off at that time.
For every time step any device can only be connected to a single segment, i.e.,: Additionally, an extra optional constraint can be added which states that any segment can only have a single device connected to it: This constraint is merely left optional in the model to allow comparison with the heuristic model, which does not take this limitation into account yet.For the real system this is a mandatory restriction.
As mentioned there are N time time intervals within the time horizon.However, within the model the time index t runs from 1 to N time + 1, whereby T t,s is the temperature of segment s of the Ecovat buffer before time interval t and T N time +1,s is the temperature of segment s at the end of the time horizon.As the decision variables, x d t,s give the decisions made during time interval t and the time interval N time + 1 lies behind the optimization horizon, the decision variables x d N time +1,s are not relevant.
The temperature distribution is constrained by the fact that the temperature in the buffer should not exceed a certain maximum allowed temperature in the buffer, T max , i.e.,: Furthermore, a constraint is added to ensure the temperature gradient from the top to the bottom of the buffer is always decreasing:

PVT Panels
Photovoltaic thermal (PVT) panels generate electrical as well as thermal energy from incident solar radiation.The efficiency of PVT panels depends on the temperature of the panels [19]; the lower the temperature the higher the efficiency.Therefore, the PVT panels are cooled by connecting them to the lowest segment of the buffer when turned on.Note that the maximum temperature allowed in the bottom segment will be low, e.g., 5 • C, to thus increase the electrical efficiency of the PVT panels.
Whether the thermal output of the PVT panels is connected to the buffer at a given time step depends on the output temperature of the water coming from the PVT panels, which should be higher than T t,N seg to be able to charge the buffer.Moreover, the model should also have the option of choosing not to connect the PVT panels if that is better.A binary variable, w t , is defined that indicates whether the output temperature of the PVT panels is high enough to be able to add thermal energy to the buffer.The value of this binary is set by the following constraints: To allow for the option of not connecting the PVT panels to the buffer even if the temperature is sufficient the following constraint is added: where is the decision variable indicating whether the PVT panels are connected to the bottom segment of the buffer at time t.
If the PVT panels are connected to the bottom segment of the buffer the thermal and electrical energy generated depends on the thermal and electrical efficiencies of the PVT panels, which will be derived below.The steady state thermal efficiency of a PVT panel is defined by dividing the yield of the panel by the total solar radiation incident on the panel [20]: where η th is the thermal efficiency of the PVT panels, ṁ is the mass flow rate of the water through the PVT panels, c p is the specific heat coefficient of water, T in and T out are the temperatures of the water entering and exiting the PVT panels respectively, G is the solar radiation and A pvt is the surface area of the PVT panel.The temperature T in is the same as the temperature of the bottom segment of the buffer.The water exiting the PVT panels with temperature T out is fed back to bottom segment of the buffer, thus charging the buffer and cooling the water back down to T in to keep cooling the PVT panels.Due to this the bottom segment warms up slightly over time, meaning one of the heat pumps needs to cool the bottom segment during some time intervals in which the PVT panels are not connected to the buffer e.g., in the evening.Alternatively, the thermal efficiency can be approximated as a linear function of the reduced temperature (see e.g., [19,21]): where η th 0 is the thermal efficiency at a reduced temperature of zero ( m 2 • • C W ), a th is the thermal loss coefficient of the PVT panels and T red is the reduced temperature which is defined as: where T amb is the ambient temperature.Combining ( 8), ( 9) and ( 10), leads to the following constraint to be used in the ILP: Figure 3 shows the output temperature of the PVT panels as a function of the ambient temperature for different global radiation levels.As mentioned before, if the output temperature is lower than the temperature of the bottom segment of the buffer the PVT panels will not be connected to the buffer.As one would expect, this only occurs for low levels of global radiation and/or low ambient temperatures.The constraint specifying the reduced temperature is then defined by: The constraints for the thermal efficiency are based on (9), with the restriction that the range for the thermal efficiency is between 0 and η th max , where η th max is the maximum thermal efficiency of the PVT panels.The following constraints give the thermal efficiency: where y th t is a binary indicator variable that determines which of the constraints ( 13) or ( 14) is active and As a consequence the value of η th t will be determined by which of the two constraints gives the tightest bound in this case.The objective function makes sure the program always picks y th t = 0 if possible and that it takes the upper bounds on the thermal efficiency in the last two cases.
The constraints on the electrical efficiency of the PVT panels are similar: where η el t is the electrical efficiency of the PVT panels, η el 0 is the electrical efficiency at a reduced temperature of zero, η el max is the maximum electrical efficiency of the PVT panels, a el is the electrical loss coefficient of the PVT panels and y el t is a binary indicator variable that determines which of the constraints ( 16) or ( 17) is active.

Air/Water Heat Pump
Heat pumps are devices that transfer thermal energy from a heat source to a heat sink, where the heat source has a lower temperature than the heat sink.To accomplish this transfer of thermal energy the heat pump requires electrical energy to operate.The air/water heat pump uses the ambient air as heat source and one of the buffer segments as heat sink.
The air/water heat pump can only charge segments that have a temperature within a given temperature range: T aw min T t,s T aw max , where T aw min and T aw max are the minimum and maximum temperatures respectively the air/water heat pump can supply.
The air/water heat pump can be modelled by the following constraints: These constraints ensure that the chosen segment lies within the allowed temperature range.Constraint (20) contains a multiplication of a binary variable, x aw t,s , with a continuous variable, T t,s , which means that in this form the constraint is not linear and can not be used in an ILP model directly.However, such quadratic terms can be linearised by using the well know method of defining the multiplication as a new variable and adding some constraints to that new variable (see e.g., [22] for more details).Note that all non-linearities in the developed model are of this form and can thus be linearised using the same method.

Water/Water Heat Pumps
The water/water heat pumps work similarly to the air/water heat pump, except that instead of using the ambient air one of the other buffer segments is used as the heat source.The configuration of the system shown in Figure 2 has two heat pumps, which are modelled exactly in the same way, only using different parameters.The water/water heat pumps have a limited temperature range T ww1 min T t,s T ww1 max similar to the air/water heat pump, where T ww1 min and T ww1 max are the minimum and maximum temperatures respectively the water/water heat pumps can supply.For the water/water heat pump both the heat source and heat sink temperatures need to be in this range.To accommodate for the fact that the water/water heat pump has to be connected to two segments the device label ww1 will be split into labels ww1w and ww1c, designating the heat sink and heat source of the water/water heat pump respectively.The decision variables that indicate whether segment s is used as heat sink or heat source during time interval t are then defined as x ww1w t,s and x ww1c t,s respectively.To model a water/water heat pump the following constraints are specified: x ww1w t,s To ensure the chosen segment lies within the allowed temperature range, To make sure the heat sink has a higher temperature than the heat source and To ensure that either both a heat source and heat sink are selected, or neither is.

Resistance Heater
The resistance heater is a device that simply converts electrical energy to thermal energy on a one-to-one basis.The modelling of the resistance heater is very straightforward, it only needs to be connected to a single buffer segment and has no temperature restrictions, except for the fact that it has to respect (3).The only constraint needed to model the resistance heater itself is (1).

Heat Demand
The predicted heat demand consists of an amount of energy demanded, D t , during every time interval t = 1, ..., N time as well as the temperature this demand should be supplied at, T dem t .The constraint that makes sure that the segment selected to supply the demand has a sufficiently high temperature is given by: T t,s T dem t where x dem t,s is the decision variable that describes whether segment s is supplying the heat demand during time interval t.The buffer is discharged to satisfy the heat demand by running (low temperature) water through the heat exchangers in the buffer walls of a segment with high enough temperature to heat this water to the demand temperature.Instead of (1) the heat demand is subject to a slightly different constraint: which makes sure that only one segment is selected to supply the demand if there is any demand, or no segment is selected if there is zero heat demand at time t.

Heat Losses
In an internal report a calculation for the heat losses of an Ecovat buffer over a period of 6 months has been performed.Assuming that the relative heat loss in any given hour within the time horizon is constant the heat lost by a buffer segment to the groundwater can be written as: where Q loss t,s is the heat lost by segment s at time t, β is the heat loss factor over 6 months, T gw is the temperature of the groundwater, m s is the mass of segment s and the factor 4380 in the exponent is the number of hours in half a year.Alternatively, if a value for the heat loss coefficient through the walls has been calculated or measured that can also be used to determine the heat losses instead: where U is the heat loss coefficient and A s is the surface of the walls through which segment s is losing heat to the surroundings.For the top segment A s also includes the surface of the lid of the buffer.At the time of this research the first method is used to compare the results with [17], however, most likely (28) will be used instead once measurements on the first Ecovat buffer are available.

Coefficients of Performance
The coefficient of performance (COP) of a heat pump determines how much thermal energy in watts will be transferred for every watt of electrical energy consumed.The resistance heater simply converts electrical energy into thermal energy on a one-to-one basis, which means it has a COP of 1.The COP of a heat pump depends on the temperatures of both the heat source and heat sink in a non-linear way.Since it is not possible to model this in an ILP the COPs have been modelled as constant for the results presented in this paper.A model incorporating COPs that depend linearly on the heat source and heat sink temperatures is a topic for future research.

Temperature Evolution
To describe the temperature evolution in the buffer a differential equation similar to the one used in [5] for their one-dimensional model is used: where m s is the mass of segment s and ∑ Q is a summation of all the heating and cooling terms from the different devices and the heat losses to the surroundings.Additionally, a term that describes heat exchange/mixing between the different segments may be added, however, for the time being this term is omitted for the sake of simplicity.Discretizing the derivative in (29) and rewriting gives the constraint: where ∆t is the length of the time intervals used in the optimization.The Q d t,s terms differ for the different devices: where N pvt is the number of PVT panels in the system, COP aw and COP ww are the coefficients of performance of the air/water and water/water heat pumps respectively and C aw , C ww and C res are the amounts of electrical energy consumed by the air/water heat pump, water/water heat pump and resistance heater respectively when they are turned on.The factor (COP ww − 1) in ( 34) is due to the fact that even though C aw • COP aw watt of thermal energy is transferred to the heat sink, only C aw • (COP aw − 1) watt is thermal energy coming from the heat source, the remaining C aw watt is the electrical energy the heat pump consumed to transfer the thermal energy.This assumes there are no losses in heat pump operation.Combining (28), ( 30) and ( 31) to (36) gives the final constraint that describes the temperature evolution inside the Ecovat buffer for all segments except the bottom segment: For the bottom segment the final constraint is given by: For the bottom segment a number of terms can be dropped from the final constraint.The heating terms by the heat pumps can be dropped because there is no lower situated segment to use as heat source.The terms from the air/water heat pump and resistance heater can be dropped because the higher situated segments will always be favoured over the bottom segment due to the objective function (see next subsection).To use equations (37) and (38) an initial temperature distribution, I s in the Ecovat buffer needs to be specified:

Objective Function
The objective of the Ecovat model is to minimize costs under the restriction that the heat demand should always be satisfied.As the gain from supplying heat to satisfy the demand is equal for all charging strategies, adding this gain to the objective function is unnecessary.The cost to use devices during a given time interval t depends on the electricity price during that time interval (which may be negative).The cost made during time interval t is given by: where K t is cost of the system at time t and E t is the electricity price at time t.It is assumed that the unused electrical energy produced by the PVT panels during time interval t can be sold for the energy price at that time.However, this could easily be changed to a different price.The total cost of the system is to be minimized.However, a few more terms are added to the objective function to ensure the system shows the desired behaviour.First a term is added that incentives the program to charge segments situated higher in the buffer over lower ones.This is desirable because only segments with temperatures higher than the demand temperature are able to satisfy the heat demand.By incorporating such an incentive the system is able to supply the heat demand also at times after the optimization horizon.This is important for the real world system, even if it gives slightly higher costs, since the system is operating continuously.The second term added is to ensure the upper bound for the thermal efficiency constraints is chosen (as discussed in Section 2.2).In general, the highest efficiencies will be selected even without this term.However, there can be cases where it would in theory be more beneficial to take a thermal efficiency of zero while having a high electrical efficiency, but this would be an unphysical situation where the panels would be cooled to achieve a high electrical efficiency without the panels heating up the water that is cooling them.This term prevents such cases from occurring in the model.The full objective function is then given by: min where c 1 and c 2 are small constants.The total cost of the system is defined as K tot = ∑ t∈T K t , which is the quantity of interest and is the one compared in Section 4 for different scenarios.

Summary of the Heuristic Model
The heuristic model described in [17] is based upon a few simple rules to determine the charging/discharging of the Ecovat buffer.For every time interval t the model determines whether a device is turned on or off.A device is turned on if the energy price during that time interval is negative, with higher prices being accepted for the heat pumps (but not the resistance heater) if the current energy content of the buffer drops below some threshold value.The PVT panels are turned on whenever they generate energy.Assignment of the segment to be charged by a device is done by simply taking the segment with the highest temperature that the device can charge.For the heat pumps the heat source is always taken to be the segment with the lowest temperature the heat pump can operate on.Due to these simple rules constraint (2), which states that a segment can only have one device connected to it, is not enforced.The other major difference compared to the ILP model is that the heuristic model only considers what happens at time interval t and bases the charging/discharging strategy on that, instead of looking ahead for future opportunities.

Optimization Setup
The ILP model described in the previous section can be used to determine an optimal steering of the charging/discharging of the Ecovat system over a longer time period.It requires a number of input parameters and profiles.The first is the overall heat demand profile of the buildings connected to the Ecovat buffer.For the results presented in the next section a case is considered which is described in [17] where an Ecovat buffer is connected to 78 apartments divided over five apartment buildings.From historical data a typical heat demand profile for a winter, a spring, an autumn and a summer day is constructed.As the spring and autumn days are quite similar for these periods the same profile is chosen.Figure 4 shows the resulting daily heat demand profiles for the different seasons.As the profiles include both space heating and tap water demand there is still a significant demand during the summer.The optimizations performed were all for a period of one year.The optimizations start on January 1, the first 60 days are winter days, followed by 91 spring days, 91 summer days, 91 autumn days and finally 32 winter days for a total of 365 days.Every day uses the heat demand profile from Figure 4 corresponding to that season.The heat demand temperature T dem t was taken constant throughout the year.Two different situations were considered with two distinct heating systems having different demand temperatures, one where T dem t = 40 • C and the other where T dem t = 60 • C. Next to the heat demand profiles, also profiles for the weather data and energy prices are required.For these profiles the ones used by [17] are used to be able to compare results between the models.For the weather data, historical data from 2005 to 2011 for the city of Eindhoven, where the apartment buildings are located, was averaged to obtain profiles for the ambient temperature and solar radiation during the year.For the energy price profile the prices on the Dutch imbalance market of 2014 were used.The resulting profiles are shown in Figure 5.It should be noted that due to the averaging of weather data over a number of years any extreme weather conditions will have been averaged out.Consequently the results presented in the next section will not include such extreme weather conditions.In future work we plan to use weather data of only a single year, which not only will include such conditions, but will also allow us to compare years with different weather conditions.Additionally, we will use the energy prices from the same year as it is expected that these are, among others, influenced by the weather conditions.
Version July 1, 2016 submitted to Energies 12 of 21 described in [17] where an Ecovat buffer is connected to 78 apartments divided over five apartment 332 buildings.From historical data a typical heat demand profile for a winter, a spring, an autumn and 333 a summer day is constructed.As the spring and autumn days are quite similar for these periods 334 the same profile is chosen.Figure 4 shows the resulting daily heat demand profiles for the different 335 seasons.As the profiles include both space heating and tap water demand there is still a significant 336 demand during the summer.

337
The optimizations performed were all for a period of one year.

343
Next to the heat demand profiles, also profiles for the weather data and energy prices are 344 required.For these profiles the ones used by [17] are used to be able to compare results between    The used input parameters for the performed optimizations are given in Table 1.Furthermore, the system configuration shown in Figure 2 is used, i.e., a configuration with two water/water heat pumps accompanying the buffer and N seg = 5 segments is used.The optimization horizon of one year is divided in time intervals of 15 minutes, resulting in N time = 35040 total time steps.
Table 1.Listing of the input parameters for the Ecovat system.* The top 3 segments have a length of 3.3 m and the bottom 2 segments a length of 2.9 m, the diameter of the buffer is 20 m and the density of water is taken as 1000 kg m 3 .+ Except for the bottom segment where the maximum temperature is set to 5 • C to be able to cool the PVT panels.

Input Parameters Ecovat
The values for η th 0 , a th , el 0 and a el are obtained from fits to measurements given in [21].Values for the maximum efficiencies, η th max and η el max have been selected slightly higher than those at zero reduced temperature, because of the positive effect cooling has on the PVT panels.
As noted before, the COP values in the ILP model have been taken as constants in this research.To determine these values a year long simulation using the heuristic model, which uses a non-linear formula to determine the COPs, was executed and the average COPs of that simulation over the entire year were taken as the COP values for the ILP.For the results in the following section these constant COP values are also used in the heuristic model to be able to make a better comparison between the ILP and heuristic model.It is worth mentioning that the difference in objective value between using a non-linear COP and using a constant COP in the heuristic model is less than one percent.
Due to the very large number of variables in the ILP model it is not possible to optimize over a whole year at once.For this reason the optimization is divided into smaller overlapping parts whereby from an optimization of each of these parts only the first part is executed.This process is depicted in Figure 6.The initial conditions of the buffer are taken and combined with predictions for the first period to optimize the charging strategy.From this solution, only the realization for the first part of the period is used.The temperature distribution in the buffer at the end of this first part is then used as the initial conditions for the next period which starts at the end of that first part, and so on.For the case shown in Figure 6 a time horizon of two days is used for every step while only the first day is executed.The effect of different horizons has been investigated and is presented in the following section.It is clear that the total objective value for the entire optimization will be higher using this method compared to an optimization over the whole year at once.However, the iterative approach simulates a more realistic use case, where good weather predictions are not available for more than a few days ahead.The stopping criteria for the ILP in every step of this procedure are: • the gap between the best bound and the best solution found by the solver is less than 0.2 % or • the solver takes longer than an hour to reach a gap smaller than this percentage. .Scheme, how the solution is build for an optimization of a year; in this case the horizon for every optimization step is two days.

Results
In this section the results of the ILP model and the heuristic model are compared for a number of scenarios, which differ only in the value of the demand temperature T dem t .The comparison of the two models is based on the achieved objective values.It should be noted, though, that the objective values alone do not tell the whole story, since for a system like the Ecovat buffer the useful energy content (i.e., energy at temperatures higher than the demand temperature) of the buffer at the end of the planning horizon is also an important quantity if the system continues to operate after the planning horizon, which is the case for a real life system.To take into account also the energy content at the end of the optimization we define the state of charge of the Ecovat buffer as the ratio of the useful energy content and the useful energy content at the start of the optimization period.Note that this means the state of charge can have a value larger than 1.
For the first scenario, the demand temperature, T dem t , is set to 40 • C. The results for this scenario are shown in Figure 7.In the figures the different lines give the temperature evolution inside the buffer.The objective value is also given in the figures.These results are presented for simulations using horizons of two, three and five days as well as for the heuristic model given in [17].
Comparing the ILP model for the different horizons shows that for a horizon of three days the optimization has a lower objective value than for a horizon of two days, however, it also has a lower state of charge at the end of the optimization as seen in Figure 8.When using a horizon of five days the objective value is a little lower than for a three day horizon while also having an almost equivalent state of charge as for a two day horizon.This means that a horizon of five days gives the best results of the three considered horizon lengths, as may have been expected.However, a five day horizon is also computationally the most expensive.The heuristic model gives charging strategies similar to those of the ILP model except that it tends to have a higher state of charge at any given time.Also the state of charge at the end of the optimization is higher than in any of the ILP solutions.However, looking at the objective value, the heuristic model gives a much higher value than the ILP optimizations.This is caused by the fact that the heuristic model simply considers the state of the buffer and the inputs at a given time without looking ahead for future opportunities.This showcases the first important improvement that should be made to the heuristic model, namely incorporating a method that allows the heuristic model to look ahead for potential future opportunities.
The second scenario uses a demand temperature, T dem t , of 60 • C. It is immediately clear from the results given in Figure 9 that the solutions of both the ILP and heuristic are significantly worse than for the case with a demand temperature of 40 • C. In general the objective values are higher, while the state of charge (see Figure 10) at the end of the optimization is much lower compared to Figures 7 and 8.This was to be expected as for a higher demand temperature the buffer has a lower useful energy content for the same temperature distribution in the buffer.When comparing the ILP solutions for different horizons the two and three day horizon solutions look very similar.However, the solution for a five day horizon has a significantly higher objective value, which may seem surprising.A closer look shows that this result is due to the stopping criteria (Section 3), in particular due to the maximum run time per step.For the case using a five day horizon one hour is simply not enough time to find good solutions for some of the steps in the scenario with a high demand temperature.Note, that for this case it takes already around 150 hours to solve the optimization over the whole year.As a consequence simply increasing the maximum run time per step by a meaningful amount is not a realistic solution to this problem.Comparing the heuristic to the ILP solution for two and three day horizons shows that the difference in objective values is a lot larger than for a lower demand temperature scenario (Figure 7).However, this is somewhat mitigated by the fact that the difference in state of charge is also larger.In this scenario the heuristic model is hurt even more by the fact that it can not look ahead for future opportunities.
A last thing to note is that in the solution given by the heuristic model two of the temperature lines cross each other.This is of course physically not possible, in reality the two segments would mix, breaking the thermal stratification, which is not desirable.As a consequence making sure that the heuristic model never exhibits this behaviour is the second important improvement that should be made to the heuristic model.
In the last scenario, we consider a demand temperature of 40 • C, but this time we enforce that the optional constraint (2) is respected.This constraint states that every segment can only have a single device connected to it for any given time interval t.The achieved results are presented in Figure 11, whereby the heuristic model is not given in this figure since it can not take this restriction into account yet.The results for this optimization are very similar to the results shown in Figure 7 except with slightly higher objective values and slightly lower states of charge (see Figure 12) at the end of the optimization due to the extra restriction placed upon the system.In this case the differences between different horizon lengths are very small both in terms of the state of charge at the end of the optimization as well as the objective value.As this constraint (2) has to be respected in practice, the third improvement that should be made to the heuristic model is to include constraint (2).

Conclusions and Future Work
In this paper an ILP model was proposed to describe the charging and discharging of the Ecovat system.This model was applied in a simulation to control an Ecovat for a scenario of a complete year.Hereby the optimization was executed in a rolling horizon fashion, whereby different time horizons per optimization step were used.When comparing the achieved results of the ILP model using different time horizons it was found that the longer horizons performed slightly better in general.However, for the scenario closest to the real system (i.e., including constraint (2)) the differences between the different horizons are almost negligible.
The results achieved with this model were compared with a heuristical approach described in [17].Furthermore, the results obtained using the heuristic model were found to be significantly worse than those obtained using the ILP model, in particular for the scenario with a demand temperature of 60 • C. Especially in this last case, the heuristic model is hurt by the fact that it only considers the current time interval and is not able to look ahead for future charging opportunities.
Based on the achieved results a number of improvements to the heuristic model have been identified and will be addressed in future work.First, in the heuristic model future time steps should be taken into account instead of just looking at the current time step.Second, the restriction that every segment in the buffer has a higher temperature than the segment below it should be incorporated into the heuristic model, to keep segments from mixing.Third, the restriction that every segment in the buffer can only have a single device connected to it should be added to the heuristic model.Furthermore, in future work COPs that are linear functions of the heat source and heat sink temperatures instead of constant COPs should be incorporated into the model.Furthermore, the model should be extended to more dimensions if measurements to the real Ecovat buffer, which is currently under construction, indicate that the one-dimensional model of the Ecovat buffer is not sufficient.Total cost of the system over the optimization period

Figure 1 .
Figure 1.Construction of an Ecovat buffer in Uden, the Netherlands[18].This Ecovat buffer is 16 meter deep with a diameter of 11 meter.However, in the future construction of even larger buffers is planned.

Figure 2 .
Figure 2. A schematic overview of the way the Ecovat system is modelled, showing the relation between the different components in the model.

Figure 3 .
Figure 3.The output temperature of the PVT panels as a function of the ambient temperature.The input temperature of the PVT panels is set at 5 • C. The solid line corresponds to a global radiation of 100 W m 2 , the dashed line to a global radiation of 300 W m 2 and the dotted line to a global radiation of 500 W m 2 .

1 Figure 4 .
Figure 4.The seasonal heat demand profiles used in the optimizations for the Ecovat system.The solid line shows the heat demand profile for winter days, the dashed line for spring and autumn days and the dotted line for summer days.

345
the models.For the weather data, historical data from 2005 to 2011 for the city of Eindhoven, where 346 the apartment buildings are located, was averaged to obtain profiles for the ambient temperature and 347 solar radiation during the year.For the energy price profile the prices on the Dutch imbalance market 348 of 2014 were used.The resulting profiles are shown in Figure 5.

Figure 4 .Figure 5 .
Figure 4.The seasonal heat demand profiles used in the optimizations for the Ecovat system.The solid line shows the heat demand profile for winter days, the dashed line for spring and autumn days and the dotted line for summer days.

Figure 5 .
Figure 5.The weather and energy prices input data that are used in the optimizations for the Ecovat system.

1 Figure 6
Figure 6.Scheme, how the solution is build for an optimization of a year; in this case the horizon for every optimization step is two days.

KT t, 1 T t, 2 T t, 3 T t, 4 T t, 5 1Figure 7 . 1 Figure 8 .
Figure 7. Temperature evolution inside the Ecovat buffer for different horizons in the ILP model and the heuristic model for a demand temperature of 40 • C.

K 1 Figure 9 .
Figure 9. Temperature evolution inside the Ecovat buffer for different horizons in the model and the heuristic model for a demand temperature of 60 • C.

Figure 10 .
Figure 10.State of charge for different horizons in the ILP model as well as for the heuristic model for a demand temperature of 60 • C.

Figure 11 . 1 Figure 12 .
Figure 11.Temperature evolution inside the Ecovat buffer for different horizons in ILP model for a demand temperature of 40 • C while adding optional constraint (2) to the ILP model.
The optimizations start on 338 January 1, the first 60 days are winter days, followed by 91 spring days, 91 summer days, 91 autumn 339 days and finally 32 winter days for a total of 365 days.Every day uses the heat demand profile 341 throughout the year.Two different situations were considered with two distinct heating systems 342 having different demand temperatures, one where T dem t = 40 • C and the other where T dem t = 60 • C.
Amount of electrical energy consumed by the air/water heat pumped if turned on C ww1 , C ww2 Amount of electrical energy consumed by the water/water heat pumps if turned on C res Amount of electrical energy consumed by the resistance heater if turned on COP aw Coefficient of performance of the air/water heat pump COP ww1 , COP ww2 Coefficients of performance of the water/water heat pumps c 1 , c 2 Small constants used in the objective function M Large constant used in some constraints T t,s Temperature in segment s at time t x pvt t,N seg Decision variable indicating whether the PVT panels are connected to the bottom segment of the buffer at time t Decision variables that indicate whether segment s is used as heat sink for one of the water/water heat pumps at time t Decision variable indicating whether segment s is supplying the heat demand at time t w t Indicates whether the output temperature of the PVT panels has a high enough temperature to charge the buffer T in Energy lost to the surroundings by segment s at time t Energy added to segment s at time t by the air/water heat pump Energy added to segment s at time t by one of the water/water Energy added to segment s at time t by the resistance heater Energy extracted from segment s to supply the heat demand at time t K t Cost of the system at time t K tot