Development of a data-driven model to characterize the heat storage of the building thermal mass in energy planning tools

In the one hand, energy planning tools compute the cost-optimal investment in the energy system minimizing life cycle costs (LCC). These tools often consider optimal control. The building (or cluster of buildings) is represented by a node where the time profiles of energy demands are given as inputs. The indoor temperate in buildings is typically not considered and may even be difficult to define for a cluster of buildings. Secondly, to perform optimization, the model of the energy system is often linear (e.g. using MILP). In the other hand, the building thermal mass has proven to be a cheap and large source of energy flexibility. Therefore, there is a need for a linear model of the building thermal dynamics when there is a limited knowledge of the indoor temperature. Consequently, the paper proposes a model that tracks the change of indoor temperature and space-heating power rather than their absolute values: the model focuses on the deviations from the reference energy profiles given as input. This framework gives a simple model that is less dependent on the boundary conditions (i.e. the weather, user behaviour and internal gains). In addition, a second-order model is proposed to characterize the transfer function. The model has only four parameters, which simplifies its identification. The model is validated using detailed building performance simulation, namely IDA ICE, on a Norwegian wooden detached house during demand response (DR).


Introduction
According to IEA EBC Annex 67 [1], the building energy flexibility is its ability to manage its demand and generation according to local climate conditions, user needs and grid requirements. Energy flexibility of buildings will allow for demand side management, load control and demand response (DR) based on the requirements of the surrounding grids and on availability of RES to minimize the CO2 emissions. For example, energy flexibility can be used for peak shaving to prevent congestion in the distribution electricity grids and avoid (or postpone) investments to reinforce them. For local energy systems, like district heating, reducing the peaks enables to minimize the size of the heat generation systems and thus the investment, or to maximize the energy use of the baseload heating generation system [2].
It has been shown that the thermal mass of buildings can be a significant short-term heat storage to perform DR [3,4]. The exploitation of such thermal storage requires the indoor temperature to fluctuate within limits acceptable for the occupants. Therefore, it is crucial to model the relation between space-heating power and the resulting indoor temperature. However, the requirement for modeling accuracy differs between applications or the specific stakeholder involved. * Corresponding author : laurent.georges@ntnu.no Energy planning tools are optimization tools aiming to determine the cost-optimal investment and operation in an energy system. One example is the energy hub approach initially proposed in [5]. Such models are used to evaluate the costs over the economic lifetime of the project, considering optimal control of the system. To limit the computational time, the optimization is often performed using mixed-integer linear programming (MILP). Energy planning models generally have a linear structure and include the time profile of energy loads, electrical and thermal, as inputs. Buildings are represented as a node in the energy system with a given energy use profile that often represents several buildings (i.e., aggregated profiles). The corresponding indoor temperature is not represented. In other words, the buildings are usually represented as an inflexible load that needs to be served. Finally, to limit the computational complexity, energy planning tools typically have a time step of one hour. The time step duration can be decreased to address problems related to energy flexibility but should be kept large enough to limit the computational time. The relatively long time step means that the short-term dynamics, like local controllers, are not modeled. Energy planning tools evaluate the optimal operation at the supervisory control level, and it is assumed the local control is done perfectly. In this framework, energy planning tools perform a direct control, imposing the power directly to the building while introducing constraints on the indoor temperature (if modeled). Since building flexibility can impact the energy system's optimal operation and expansion plan, it is important to represent this aspect in the planning tools.

Literature review
A few authors have investigated how to integrate flexible assets into energy planning tools. Junker et al. have [6] proposed a linear function to characterize energy flexibility. The method has recently been generalized to nonlinear systems [7]. In this last paper, the energy flexibility is evaluated with stochastic differential equations describing the change of state-ofcharge (SOC) as a function of the change of power compared to a reference scenario. In this paper, this framework has been applied to the specific case of the thermal mass energy storage. Other authors have instead developed a methodology mainly dedicated to the thermal energy storage. In this context, it should then be distinguished between methods that evaluate the absolute value of the indoor temperature (category A) and methods that evaluate the change of indoor temperature compared to a reference scenario (category B). Thermal comfort criteria are typically expressed as a range of indoor temperature bounded by a minimum and maximum threshold, in absolute value [8]. To be able to model the indoor temperature, it requires models that incorporate many complex physical phenomena. Essential aspects to model the underlying physics of a building include the influence of solar irradiation, the heating and ventilation system dynamics. More critically, the user behaviour influences the internal gains and the selection of the set-point indoor temperature. For instance, Patteuw et al. introduce a methodology [9,10] to model the building thermal storage in energy planning tools. A recent contribution of Hedegaard et al. [11] proposes a clustering to reduce the number of buildings in the optimization and thus the computational time.
However, reformulating the problem in relative terms (category B) leads to important simplifications. The change of indoor temperature compared to the reference scenario (i.e., the SOC) is evaluated as a function of the shift in space-heating power. This is based on the idea that the building dynamics is complicated but that the thermal mass energy storage only corresponds to a limited indoor temperature change. Therefore, the building dynamics are linearized around the reference scenario, and a linear model should be accurate enough for small deviations from the reference. This idea has been proposed in the paper of Kensby et al. [2] and Romanchenko et al. [12] in the context of district heating. The energy flexibility is also defined as a change in energy demand compared to a reference operation scenario in [4,6,13].

Contribution and scope
This paper proposes a fit-to-purpose linear model to represent the dynamics of the building thermal mass in energy planning tools, with sufficient accuracy and with compatible computationally cost. In addition, the datadriven model of category B enables to characterize the flexibility potential of building thermal mass.
The paper is organized as follows. Firstly, the concept of linearizing changes of temperature and power in relative terms is formalized. The linearization is a simple and known concept; however, its consequences on the thermal mass energy storage have not been adequately discussed in the literature. Secondly, the paper proposes and compares simulation procedures to evaluate the transfer function between space-heating power change and the change of indoor temperature. Thirdly, the paper proposes to characterize the energy flexibility using a second-order model stepresponse. The step response is used as a signature of the storage performance and is parametrized using the overall heat transfer coefficient as well the two main times constants of the building envelope. Finally, the simulation accuracy of the proposed model is compared with detailed dynamic simulations in IDA ICE [14] on a single-family house.

Linearization of the building dynamics
As a starting point, it is crucial to define the system boundaries. If the energy flexibility is activated, part of the energy system should be explicitly modeled in the energy planning tool. Regarding the energy storage in building thermal mass, the modeling includes the building envelope and the heat emitter, excluding its local control (e.g., a thermostatic valve). The objective is to obtain a linear model to represent the relation between the increase of power delivered to the heat emitter and the increase of indoor temperature.
With these boundaries, the building dynamics can be described by a nonlinear state-space model: (2) where X is the state vector, Y is the scalar output, namely the indoor temperature, Tin, and U is the input vector with the outdoor temperature (Tout), the wind speed, two components of the solar irradiation, such as the global irradiation on a horizontal plane (Ig,h) and the beam irradiation on a normal plane (Ib,n), the power of the internal gains (Pint) and the power delivered to the heat emission subsystem (Pe). The flexible operation can be expressed based on the reference scenario as: The nonlinear state-space model can then be linearized around the reference scenario: w X X X X (6) ' ' w X X X (7) This is a linear state-space model with a single input and a single output (SISO), ΔPe and ΔTi, respectively. It has the corresponding transfer function H(s): (8) Based on this, these important characteristics can be expressed: x Compared to the original nonlinear system of equations (Eqs 1 and 2), the linearized model only models the influence of ΔPe on ΔTi and phenomena that are complex to model, like solar gains and internal gains, are excluded, meaning that these phenomena are already taken into account in the reference scenario.
x Even though the deviations are evaluated with a linear model, the reference scenario can still result from a nonlinear model. x The matrix A, B, C are expected to be weakly dependant on the reference scenario considered. Consequently, the matrices A, B and C can be identified using a first reference scenario and then used in the energy planning tool for another reference scenario. However, this hypothesis is questionable when highly nonlinear phenomena occur in the building, like active solar shading controlled with indoor temperature or variable effectiveness of the ventilation heat recovery.
x An abundant literature has demonstrated that a firstorder model is not accurate enough to capture the thermal dynamics of the building envelope, see, e.g., [2,12]. A second-order model with fast and slow dynamics have proved to give good simulation performance.
We can then distinguish between two application cases: x Tin,ref is known. The constraint on thermal comfort can be imposed on Tin using Eq. 4. The advantage of the method is that the reference scenario can be obtained from measurements or detailed simulations. As the reference scenario is not optimized (e.g. using MILP), the model complexity to evaluate the reference scenario can be significantly higher than linear models, like Urban Energy Modelling tools (UBEM) [15]. This gives a better chance to capture the influence of solar and internal gains correctly than using models only relying on linear equations (category A). x Tin,ref is unknown. This is typically the case for energy planning tools. The constraint on thermal comfort should then be expressed in terms of temperature deviation (ΔTi). If this framework is acceptable, the modeling is greatly simplified. Also, the method is relatively easy to scale to a group of buildings. If several buildings are identical (meaning with a same H(s)), but the user behaviour and weather conditions different, the method considers all the buildings equivalent and they can can be aggregated directly.

Simulation-based identification procedures
A significant drawback of the method is related to the model identification. It requires to have Tin and Pe measurements for the reference scenario and the deviation from this scenario when the exact same boundary conditions (i.e., input vector, U) are applied to the building. This information is difficult to obtain using field measurements (see discussion in next Section 2.3). However, this is straightforward to obtain using simulation. Two procedures can be proposed. The first procedure (method M1) is the most elaborated. It is divided into two main steps: one closedloop and one open-loop simulation (see Figure 1).  The second procedure (method M2) is easier to implement. The building is simulated without solar and internal gains and with a constant outdoor temperature. The power delivered to the space-heating emission subsystem (Pe) is directly imposed (i.e., closed-loop). The simulation is started with a constant Pe,ref until the indoor temperature reaches a steady-state, Tin,ref(t). Then, an excitation signal ΔPe(t) is added to Pe. The resulting change of indoor temperature is then registered, ΔTi(t). Compared to method M1, the reference scenario of method M1 is much more specific and can be far from the usual operating conditions. Using these procedures, time series of ΔPe(t) and ΔTi(t) are then available. The transfer function of equation (8) can be identified using linear black-box or grey-box models, in deterministic or stochastic forms.

Characterization using step response
As a complement to the methodology, we propose to characterize energy flexibility using a step response. It means that ΔPe in Eq. 6 is taken as a step function. The step response ΔTi(t) will represent the signature of the building energy flexibility potential. Several authors have stressed the close relationship between the thermal mass energy flexibility and its time constants [2,6,12,16]. In the same way, the step response is parametrized as a function of the building's time constants, meaning two constants for a second-order model: Utot is the overall heat transfer coefficient, τ1 and τ2, the two building time constants, and α is a weighting factor for the relative importance of each time constant in the step response. This parameterization of the step response is expected to be more universal and ideally to not vary much within a given building category or archetype. Firstly, Utot is a steady-state performance indicator and can be evaluated with simple evaluation methods. Secondly, the characteristics of the thermal dynamics are here translated by the time constants. This is believed to be a better choice than thermal capacitances (C). C, like heat transfer coefficients (U), depends strongly on the building size. The time constant is a ratio of C and U and is thus expected to be less dependent on the building size. The step response can be converted into a linear state-space model using different methods. A 2R2C grey-box model with four parameters can be used. For simplicity, the conversion will be done directly using the transfer function H(s). The time constants are directly related to the eigenvalues of the linear system: To simplify the notation, a weighted eigenvalue can be defined: The transfer function (Eq. 8) corresponding to the step function (Eq. 9) has then following expression: The transfer function can be converted back into the physical domain, firstly into a second-order ordinary differential equation and, afterwards, into a secondorder linear state-space model: The system can be identified in different manners. Firstly, ΔPe can be applied, and the resulting step response can be directly used to determine the four parameters using Eqs. 9 and 10. Then, the linear statespace model can be obtained with Eqs. 14 and 15. On the contrary, more sophisticated excitation signals can be used for ΔPe, like a Pseudo-Random Binary Sequence [17]. The parameter of Eqs. 14 and 15 can then be identified, like a grey-box model. The characterization curve (Eq. 9) is obtained as the parameters are known.
It has been explained in Section 2.2 that the method is difficult to apply using field measurements rather than simulation. However, using the energy flexibility characterization in Section 2.3, several existing approaches can be used to determine the model parameters using field measurements: x The recent paper of Palmer Real et al. [16] proposes to identify the two time constants during the spaceheating temperature setback during nighttime, when internal and solar gains should be negligible and the space-heating power equal to zero. They develop a stochastic approach to identify the time constants based on the indoor temperature decay curve.
x In Kensby et al. [2], the authors create an average profile by analysing several days. This profile will then serve as the reference scenario. x A standard second-order grey-box model can be identified with all the terms of the U vector; see, e.g., Madsen et al. [18,19]. The A matrix is taken, but only the term of the B matrix corresponding to the contribution of Pe is included (while the other elements of B are disregarded).
Finally, the aggregation of several buildings is straightforward. If N buildings have the same parameters (Utot, τ1, τ2 and α) but different user behaviour or weather conditions, a same increase of indoor temperature will take place in each of N buildings if the same increase of power is applied. In the energy planning tool, an optimization constraint is set on ΔTi and the change of power at aggregated level, ΔPe,tot, is divided equally between the N buildings:

Description of the building
The case of a single-family detached house is analysed using detailed dynamic simulation in IDA ICE [14]. The building has been initially defined and implemented by Rønneseth et al. [20] and is available for different insulation levels. Some adjustments of the model have been applied by Elin Storlien [21]. The building is heated by electric radiators, here represented by ideal heaters controlled by a PI controller. A real electric radiator has most likely an on-off control. However, energy planning tools consider large time steps and assume the local controller to be perfect. Besides, if many similar buildings are aggregated, each radiator's short-term dynamics related to the on-off control will be averaged. Consequently, a PI control is assumed to be reasonable. As the ideal heater has no thermal mass, the heat delivered to the ideal heater (Pe) equals the heat delivered to the room. The ratio between convection and total heat emitted is fixed to 0.4. The building is divided into three thermal zones: the ground floor with the living room, the bedrooms, and the bathroom on the first floor. Each zone has one heat emitter. The nominal power (Pe,n) of each emitter has been sized using a design heating load simulation. The internal doors between zones have been assumed closed.
The building is constructed in wood, meaning a lightweight structure. Two different performances of the building envelope are considered. Firstly, a superinsulated version where the building complies with the Norwegian definition of the passive house standard [22]. The building then has balanced mechanical ventilation equipped with heat recovery with constant effectiveness taken at 85%. Second, a poorly insulated version corresponding to the Norwegian building regulation requirements of 1987, TEK87 [23]. In this case, the building has natural ventilation, here modeled as balanced mechanical ventilation without heat recovery. Some comments should be given regarding the modeling of the building physics in IDA ICE. The tool automatically integrates a ventilation network model [24]. Pressure coefficients (Cp) have been defined on each external wall. Consequently, wind-and buoyancydriven air infiltrations are computed. The heat conduction in walls assumes constant thermal properties, which makes it linear. However, surface convection coefficients and thermal radiation between surfaces are nonlinear. As the TEK87 is less insulated, the influence of the surface heat transfer is more important and the infiltrations are larger. Consequently, the thermal dynamics of the TEK87 building is expected to be more nonlinear than the passive house. The building has no active solar shading but is surrounded by obstacles representing neighbouring buildings. The influence of these obstacles on the direct solar irradiation is evaluated in detail at every time step by IDA ICE. In all simulations, the data import and export in IDA ICE are performed using a sampling time (Δt) of 6 minutes.

Building identification procedure
The building step response (Eq. 9) will be evaluated using the two methods described in Section 2.2. In the present paper, the building is directly excited by a step response.
In method M1, the building has the usual operating conditions (reference scenario 1). Internal gains are taken from the technical standard SN:TS3031 [25] and the weather file from a typical meteorological year (TMY) of Oslo. For the closed-loop simulation (step 1), the living room has a constant set-point temperature (Tday) at 22°C for the passive house and 21°C for the TEK87 building, the bedroom (Tbed) at 18°C and the bathroom 24°C (Tbath). The building is simulated for three months, from the beginning of December to the end of February. The period is equally divided into two: one initialization period and the analysis period. The Pe during the initialization period is the same for the closed-loop and open-loop simulations. However, during the analysis period, a step increase of the emitted power ΔPe in the living room is applied to the open-loop simulation. The magnitude of ΔPe is taken so that the steady-state ΔTi,∞ is about 2°C. As shown in Figure 3, there is a considerable variation of the outdoor temperature during the analysis period. In method M2, the building has no internal and solar gains, and the outdoor temperature is fixed at -5°C (reference scenario 2). The building is only simulated in open-loop. During the training period, a constant Pe,ref is imposed on the three heat emitters. A step increase in the radiator's power in the living room is applied in the analysis period. The baseline ΔPe is taken so that the steady-state ΔTi,∞ is about 2°C. A more severe increase (called a high scenario) is also tested. ΔPe is then doubled so that ΔTi,∞ is about 4°C.
The measured step response is used to identify the parameters. The entire period is used to identify Utot using Eq. 10, while the first 6 hours are used to determine the time constants, τ1 and τ2, with Eq. 9. The nonlinear regression function of MATLAB, nlinfit, is used, minimizing the squared error.

Demand response scenario
To validate the modeling approach, the indoor temperature computed with IDA ICE and the secondorder model will be compared during a DR event. A third reference scenario demonstrates that the reference scenario can be taken differently during the application of the model than during its identification. In reference scenario 3, the building has the same internal gains and weather file as reference scenario 1. However, a night set-point temperature setback is applied in the living room. The set-point (Tday) is at 22/21°C during daytime (here 7 AM to 11 PM), 16°C during nighttime (Tnight). The bedroom set-point temperature (Tbed) is reduced from 18°C to 16°C. The DR case is defined as follows. Like in Clauss et al. [26], typical peak hours in the Norwegian electricity grid are between 7 and 10 AM and between 5 and 8 PM. The peak-shaving is performed with a simple rule-based control. The building thermal mass is pre-heating three hours before the morning peaks and two hours before the evening peaks. The intensity of the pre-heating is done to decrease the electric space-heating during the pre-defined peak hours. The amplitude of the decrease has been selected to give daily fluctuations of the indoor temperature in the range of 5°C. With this approach, more space-heating energy is shifted relatively in the passive house case than in the TEK87 case. This DRcase can be seen in Fig. 4 as the black solid line. The demand response is applied during a period equivalent to the analysis period (defined in Section 3.2). The space-heating energy use during the pre-defined peak hours is defined as Eref while the amount of energy shifted from peak hours is called ΔE. ΔE can be expressed in absolute (ΔE) or relative terms (ΔErel):

Linearity hypothesis
The overall procedure is based on the idea that a linearization of the nonlinear state-space model is enough to analyse thermal mass activation. It will be checked if the step response is dependent on the reference scenario, the amplitude of the step or if the building is unloaded or loaded. If the thermal dynamics is linear, the step function should be symmetric for loading and unloading events. To make a comparison, all the step responses are scaled to the same ΔPe,norm here arbitrarily taken at 100W for the passive house and 250 W for the TEK87 building:  Step response for the passive house using identification methods 1 and 2: scaled for the same ΔPe,norm.
The data series of the step response is first analysed for the passive house. Figure 5 shows the step response for the first 6 hours. All the curves are very similar. It means that the transfer function H(s) (Eq. 8) is not much dependant on the reference scenario or the amplitude of the step function ΔPe. The cases for loading and their corresponding case in unloading are almost symmetric. Hence, it can be argued that the linearity hypothesis is reasonably valid. The identification method 2 did not lead to different results than method 1 even though the reference scenario 1 is significantly different than the reference scenario 2.

Figure 6.
Step response for the building with low insulated building (TEK87) using identification methods M1 and M2: scaled for the same ΔPe,norm. The analysis of the TEK87 building should be more demanding as the physics is expected to be more nonlinear. In Figure 6, the step function is not much influence by the reference scenario (i.e., method M1 vs. M2). However, a difference appears between loading and unloading cases, which are not fully symmetric anymore. A possible explanation for this is the nonlinearity of the convection coefficient for the indoor surface of the walls.
As proposed by Palmer Real et al. [16], the temperature decay during nighttime can be used to evaluate the two main time constants of a building. Using a similar approach, the temperature decay from the unloading case of Method M2 (using step function) can be compared to the temperature decay taking place during the 34 nights in the reference scenario 3, see Figure 7(a). The curves are normalized using the spaceheating power taken before the sudden reduction of the set-point temperature from 22/21°C to 16°C at 11 PM. In reality, the outdoor temperature cannot be assumed constant during the night temperature setback. Also, internal gains contribute to the space-heating. This makes the ΔPe at the beginning of the temperature setback challenging to define. Consequently, a different curve is obtained for each temperature decay. The mean value and standard deviation are then evaluated and shown in Figure 7(b). It can be seen that the mean profile of the temperature decay can be taken to identify the step-response (if the change of internal gains during the night-setback is properly taken into account).  The step-response in Eq. 9 is now fitted to the time series. The case of the passive house is shown in Figure  8. It can be seen that the second-order response is able to match well with the virtual experiments. Unlike a first-order model (obtained by imposing coefficient α to be 1), the second-order model can capture the short-term dynamics and long-term behaviour. Results are thus in line with this well-known conclusion in building physics, see, e.g. [2,11,12,27].

Demand response (DR)
The performance of the identified models to predict the indoor temperature during DR can now be compared to IDA ICE. The linear state-space model of Eqs. (14) and (15) is simulated using the lsim function of MATLAB [28]. IDA ICE assumes that variables are linear between discrete values at each time step. To be consistent with this, the identified models are simulated in lsim using a first-order hold for the B term. The identified model captures the change of indoor temperature. Nonetheless, it has been decided to show the absolute value of the temperature to facilitate the physical interpretation of the results. The case of the passive house is analysed first in Figure 9. This case is expected to give the best simulation performance. As mentioned, the change of power has been designed to limit daily indoor temperature fluctuations to 5°C. For the passive house case, this is enough to remove the space-heating needs during the pre-defined peak hours, see Table 3. It is questionable if such large temperature fluctuations are acceptable for building occupants. The objective with the validation case was rather to challenge the model using large variations. It can be seen that the second-order linear model can reasonably reproduce the predictions of IDA ICE. There is not a significant difference in simulation performance during loading and unloading events. The error is well below the expected accuracy, as the uncertainty on the temperature limits acceptable by occupants is much more uncertain. The case of the TEK87 house is shown in Figure 10.
The second-order model is also able to fairly reproduce the predictions of IDA ICE. As shown in Tables 4, only a fraction of the space-heating needs is reduced during peak hours. month of the space-heating season. The maximum absolute error (MAX), the mean bias error (MBE), the mean absolute error (MAE), and the goodness of fit (G) are shown. For example, the definition of these error indexes can be found in Afram et al. [29]. In general, the simulation performance is fair, as confirmed by values of G higher than 80%. As already mentioned, the identification using a step function is a relatively simple procedure. It means that there is a potential to calibrate better models using a better excitation signal, like PRBS. Instead, the main objective was to demonstrate that Eqs. 9 and 10 are an efficient way to characterize the energy flexibility potential of the building thermal mass. An essential aspect is that the simulation performance does not vary significantly between the different months of the space-heating season. The model has been calibrated using very specific operating conditions (i.e., the reference scenario 2 in Table 1). Nonetheless, the model performs equally well in all the months. This indicates that a linear time-invariant statespace model of second order is enough to model the relation between ΔPe(t) and ΔTi(t).

Discussion and Conclusions
A series of comment can be given regarding the hypothesis of linearity: x The difference between loading and unloading events is not specific to the present modeling framework. It is related to the linearity assumption. Therefore, any linear time-invariant state-space model would face the same issue. It does not matter if Ti or ΔTi is modeled.
x It is believed that assuming no active solar shading controlled on indoor temperature is a reasonable assumption. However, the assumption of constant effectiveness for heat recovery is more questionable. Firstly, the building is equipped with a constant air volume (CAV) ventilation. The airflow rate should not then lead to a change of effectiveness. This can be more challenging for variable air volume (VAV) ventilation systems, often implemented in office buildings. Secondly, for simplicity, the natural ventilation of the TEK87 building has been modeled by mechanical ventilation without heat recovery. In reality, the physics of natural ventilation is nonlinear. Furthermore, the internal doors between rooms have been assumed closed. With open internal doors, the bi-directional airflow between zones is also highly nonlinear.
The test case has been taken multi-zone. Firstly, in some countries, the temperature difference between zones can be considerable. For example, it is known that many Norwegian like colder bedrooms, even colder than 16°C [30]. The user preference regarding indoor temperature only influences the evaluation of the reference scenario. ΔTi is independent of the reference indoor temperature set-point defined by the occupants. As for internal gains, this makes the method less dependent on the user-behaviour.
In conclusion, the proposed data-driven model proved to be accurate enough to predict the change of indoor temperature resulting from a change of spaceheating power. The formulation of the model makes the aggregation to several buildings straightforward, and the analysis of the energy flexibility of neighborhoods possible. In future work, the method can be implemented in an energy planning model to optimize thermal mass activation as a part of a larger energy system.