Probabilistic FlexOffers in Residential Heat Pumps Considering Uncertain Weather Forecast

,

(1) Decrease dependency on fossil fuels supplying space heating and domestic hot water consumption. (2) Facilitate the integration of renewables into residential heating systems.
To achieve the aims, heat pumps can be operated to provide flexibility for the upstream energy networks. The flexibility potentials of residential heat pumps stem from the thermal inertia of the buildings. Also, thermal storage, e.g., water tanks, plays a pivotal role in increasing the flexibility potential of heat pumps. Occupants' thermal comfort is normally defined as a temperature setpoint with upper and lower thresholds. Consequently, the heat pumps can be operated near lower-/upper-temperature thresholds when the upstream network encounters a shortage/excess of renewable power.
In the last decade, many studies have been conducted to address the flexibility potential of residential heating systems. In 2010, the ground source heat pumps were studied with solar collectors in Sweden to optimize the energy consumption of the heating systems (Kjellsson et al. 2010). The simulations were coded in TRNSYS software. The study concluded that boreholes can influence the energy consumption pattern of households. In 2011, a novel photovoltaic model was designed in the UK to optimize the operation of heat pumps (Zhao et al. 2011). Based on the simulation results, the suggested approach achieved 55% and 19% of thermal and electrical efficiency, respectively. In 2012, a research study compared the operation of boilers and heat pumps for intermittency-friendly energy systems in Denmark (Blarke 2012). Simulation results showed that well-designed heat pumps are more cost-effective than electric boilers especially when the energy price increases. In 2013, the application of domestic heat pumps was developed for space heating and cooling patterns (Gupta and Irving 2013). The study used the regression methods addressing the UK standard model for household energy consumption. Also, it evaluated the emission reduction from the optimum operation of heat pumps. In 2014, a tariff-based energy price approach was discussed to shift the operation of heat pumps to off-peak periods (Kelly et al. 2014). The approach was examined in UK's detached dwellings. In 2015, the flexibility potentials of heat pumps were integrated into upstream networks to optimize the operation of smart grids (Kreuder and Spataru 2015). As the simulation results revealed, the demand response capacity of heat pumps has the potential to prevent the new peak demand on the grid level in winter. In 2016, the heat waste recovery of municipal wastewater treatment plants was used to provide stable and flexible operation for heat pumps (Chae and Ren 2016). Based on the results, the Coefficient of Performance (COP) of the heat pump reached optimum points of 4.06 and 3.64 for heating and cooling purposes, respectively. In 2017, a methodology was demonstrated to explore the flexibility of heat pump pools (Fischer et al. 2017). The study concluded that heat pump flexibility mainly depends on outdoor temperature and time of day. Although the electric backup heater provided power flexibility for shiftable energy, it decreased the energy efficiency of the whole system. In 2018, a twostep optimization model was presented to quantify the flexibility potential of residential heat pumps considering dimensions of time, energy, and cost (Oluleye et al. 2018). The study was implemented on 445 prosumers. The results indicated that the heat pumps can absorb a surplus of renewable energy, especially when integrated with thermal storage. In 2019, a Model Predictive Control (MPC) approach was implemented on variable speed heat pumps to unlock the flexibility potentials of buildings (Péan et al. 2019). The results showed that flexible control of heat pumps can provide load-shifting, cost and emission reduction, and energy flexibility.
With increasing attention towards heat pump operation, more recent studies have been focused on the flexibility opportunities of heat pumps. In 2020, the role of electricity price was investigated in the flexibility of residential heat pumps integrated with thermal storage (Fitzpatrick et al. 2020). The study suggested that real-time pricing is the most favorable scheme to unlock the highest flexibility potentials with the lowest cost. In (Bruninx et al. 2020), it was stated that the thermal capacity of residential heating systems is insufficient to provide enough flexibility in electricity markets; therefore, aggregators are proposed as mediators to monetize the heat flexibility of residential buildings through demand response programs. In 2021, a stochastic optimization approach was carried out to integrate the flexibility potentials of heat pumps into three uncertain electricity markets, including day-ahead, intraday, and balancing markets (Golmohamadi 2021). In research study (Feldhofer and Healy 2021), key approaches were characterized to increase the flexibility of detached houses under different weather conditions. In this study, the cumulative absolute residual load was introduced as the flexibility metric. The results showed that 36% flexibility is obtained annually in comparison to the base case. In 2022, an MPC was suggested for residential buildings supplied by low-temperature district heating. As the simulation results revealed, the heat controller optimizes the heat consumption of buildings in response to variable energy prices. Therefore, it shifted some part of energy consumption to off-peak periods, i.e., low price hours (Golmohamadi and Larsen 2021). An extensive review was conducted to classify the flexibility potentials of heat pumps, buildings, and district heating systems (Golmohamadi et al. 2022). Also, different generations of heat controllers were described in detail. The research study (Rasmussen et al. 2022) suggested that heat pumps can provide fast response flexibility for primary frequency support in power systems.
What is missing in the literature is to generate energy flexibility for heat pumps supplying residential buildings with different temperature zones under uncertain weather conditions. To fill the gaps, this paper generates FlexOffers for individual heat pumps in residential dwellings. To unleash the heat flexibility, the thermal dynamics of residential buildings are estimated using a data-driven approach. The estimated thermal dynamics are used to calculate the minimum and maximum flexibility potentials of the building in pessimistic and optimistic thermal energy consumption patterns. This way, the suggested approach is comprised of two main parts as follows: • Building model: in this section, the thermal dynamic model of the residential building is extracted using high-fidelity input data. The thermal dynamics are extracted using grey-box models in CTSM Software (Continuous-Time Stochastic Model). The output of the software is the estimation of the thermal dynamic and mathematical model of the buildings. The software is coded in the R language. • FlexOffers: in this section, the building model is integrated with an objective function to generate FlexOffers. The approach includes a success function to minimize and maximize the energy consumption of heat pumps. The FlexOffers are generated considering the existing gap between the minimum and maximum energy consumption.
The FlexOffers are formulated as queries on model strategies in UPPAAL-STRAT-EGO (UPPAAL STRATEGO Software 2022). UPPAAL is a free and not-for-profit software tool for research and teaching purposes that can be downloaded from the homepage (UPPAAL STRATEGO Software 2022). According to UPPAAL STRAT-EGO Software (2022), "UPPAAL-STRATEGO facilitates the generation, optimization, comparison as well as consequence and performance exploration of strategies for stochastic priced timed games in a user-friendly manner. The tool allows for efficient and flexible strategy-space exploration before adaptation in a final implementation by maintaining strategies as first-class objects in the model-checking query language. " • Based on the abovementioned facts, the main contributions of the study can be stated as follows: • Integration of the thermal dynamic model of residential buildings with multi-temperature zones from R language into UPPAAL software. The building's input data are used to train the data-driven approach, i.e., the CTSM. • Generating FlexOffers through calculation of minimum and maximum flexibility potentials of heat pumps in pessimistic and optimistic energy consumption patterns in UPPAAL-STRATEGO. • Investigating the impacts of the uncertain weather forecast on the flexibility potential of heat pumps.

Problem formulation
In this section, the suggested approach is formulated mathematically. First of all, the thermal dynamics models are described. Afterward, the success function and FlexOffer generation are stated.

Thermal dynamics of buildings
Thermal dynamics of buildings describe how the temperature of indoor air, envelope, and heaters are changed in consecutive time slots in response to the heat consumption of the heating system and weather variables, e.g., ambient temperature and solar irradiation. In this section, the thermal dynamics are described using the following three differential equations ): (1) where the parameters and variables of the mathematical model are introduced in Table 1. Also, the notations are described graphically in Fig. 1.
The three differential Eqs. 1-3 describe consecutive changes in the temperature of indoor air, envelope, and heater, respectively. The heater may denote radiator or floor heating. In Eq. (1), the first term describes the heat transfer between heater and indoor air; the second term is the heat flux between indoor air and envelope; the third term is the heat energy captured by solar irradiation. In Eq. (2), the first term states the heat exchange between indoor air and envelope; the second term is the heat flux between the envelope of the target room with adjacent rooms; the third term denotes the heat exchange between the envelope of the room and envelopes surrounded by the outdoor environment. In Eq. (3), the first term is the heat transfer from the heater to the indoor air; the second term denotes the heat extraction from the heaters. Heat capacity of floor pipes for room r (kWh/ o C)

R ie
Heat resistance between indoor air and envelope ( o C/ kW)

R r ie
Heat resistance between indoor air and envelope for room r ( o C/ kW) Heat resistance between envelopes of rooms r and r' ( o C/ kW) Heat resistance between envelopes of rooms r and ambient ( o C/ kW)

R r ih
Heat resistance between indoor air and floor pipes for room r ( o C/ kW)

R ea
Heat resistance between envelope and ambient ( o C/ kW)

Estimation of thermal dynamics
In the thermal dynamics, Eqs. 1-3, constant coefficients should be estimated for the target building. The thermal coefficients include heat resistance, heat capacity, and the fraction of solar power absorption. The set of constant coefficients is stated as follows: The thermal coefficients (4) depend on the building characteristics, e.g., dimension, windows, insulation quality, envelope material, etc. The crux of the matter is that the estimation of thermal coefficients requires sophisticated machine-learning approaches with adequate sensor data.
To estimate the coefficients, a grey box approach is used based on CTSM. The CTSM is coded in R language by DTU (Technical University of Denmark) and is publicly available (DTU Compute 2021). The software estimates the constant coefficient of the thermal dynamics using measurement data. The input (measurement) data are stated as follows (Bacher and Madsen 2011): where k is the point of data measurement time; Y k is the measurement data; e k is the measurement error.
The input sensor data is imported to the grey box of CTSM to obtain maximum likelihood estimation as follows (Bacher and Madsen 2011): where p(.) is a conditional density to describe the probability of observing the sensor data Y k given the previous observations and the target parameters θ; and p 0 (.) denotes the initial conditions. Finally, the maximum likelihood for the target parameter θ s is found by (7).

Probabilistic flexOffers
Probabilistic FlexOffers are an extension to the FlexOffers concept, in which electric consumption flexibility of consumer and/or prosumer is extracted, and then traded in the energy market with other peers participating in the electric community. What this new concept of Probabilistic FlexOffers adds to the original one, is the quantification of the uncertainty of the FlexOffers bounds; this uncertainty on the energy consumption/production flexibility can for instance be due to deviations from forecasted weather conditions or building consumption patterns. Considering the unpredictable conditions during the process of FlexOffers generation and trading, allows us to provide valuable information on the likelihood of prediction errors, and to increase or reduce the flexibility interval based on desired confidence. Figures 2 and 3 show an example of the representation of a slice of a probabilistic FlexOffers. The minimum and maximum bounds of the slice are expressed by probability distributions, Normal distributions in this example. Let the minimum/maximum consumption distributions be referred to as min and max, respectively. Then, for each energy input x, the schedule success function succ is given as follows: where CDF refers to the associated Cumulative Distribution Function. The function describes the probability that the system is able to follow a given schedule.
This success probability distribution function comes into the problem when FlexOffers is assigned to the consumer/prosumer offering this flexibility amount. Indeed, if the buyer assigns a schedule with low probability, then a lower penalty is incurred for not following the schedule. This can be used by the buyer to evaluate if they are willing to take a risk if a grid overload is severe enough.

Building thermal dynamic model mapping
As stated before, the building thermal dynamic is fully described by the three Ordinary Differential Equations (ODEs) which express Indoor Temperature (T i ), Envelope Temperature (T e ), and Heater Temperature (T h ) behaviors in each room presented in the building. Once the thermal dynamic coefficients have been estimated, these ODEs are imported into UPPAAL and used to predict the temperature evolution in response to the model inputs. A few of the most common model inputs which are usually kept into consideration are heating system output thermal power (P h ), solar irradiance power (P s ), and outdoor ambient temperature (T am ).

Weather forecast stochasticity
An important step forward is given by the way stochasticity is added to weather forecast data. Indeed, in the previous approach (Agesen et al. 2017), the weather forecast sample (e.g., ambient temperature) at time t was extracted from a Gaussian Distribution with a mean equal to the value retrieved from weather historical data, and fixed variance set by the user. The issue is that with that approach, samples at times t and t + 1 are uncorrelated, and this is not a realistic hypothesis.
In the proposed approach, two consecutive samples are correlated with each other, and the overall sample trend follows the weather forecasted behavior.
Starting from a weather forecasted trend retrieved from some meteorological institute, the data with added stochasticity (denoted with the apex "s") are computed using Eqs. (9) and (10). The variable err(n) is defined as the current deviation from the forecasted data; it is extracted from a Gaussian Distribution centered on the previous err value, and with a variance that is proportional to how much time in the future we are predicting the data sample (i.e., the more in the future we are trying to predict, the less reliable can be the forecast data). Once the current err has been calculated, it is summed to the weather forecast data sample, but given that we do not want to deviate too much from the forecasted trend, we will not add an error greater than err max . In the case taken into consideration in Eqs. (9) and (10), the stochasticity is added to the forecasted ambient temperature T a . Then, T a s represents the generated (stochastic) data. Figure 4 depicts three different outdoor temperature behaviors obtained through the beforehand explained algorithm starting from the same weather forecast data. The three subfigures represent a case study where the difference between forecasted and uncertain data is marginal (top-left plot), and two cases where the generated expected temperature is particularly lower or higher than the forecasted one (top-right and bottom plots, respectively). The same method can be applied to any other weather variable (e.g., wind speed); indeed, in the following, we use this method to add stochasticity also to the solar (9) err(n) = N(err(n − 1), σ ) (10) T s a (n) = T a (n) + err(n) if |err(n)| < err max T a (n) + sign(err(n)).err max otherwise irradiance power forecast. Note that the values of σ and err max to generate the stochastic data for T am and P s are stated in Table 2.

Scheduling strategy and probabilistic flexoffers generation
In this paper, it is assumed that the model includes two continuous variables kWh and time representing the overall energy consumption and the global time, respectively. The optimization capabilities of UPPAAL-STRATEGO are then used to generate two strategies, including σ H min and σ H max , that minimizes (resp. maximize) the expected value of kWh for the horizon H.
These two strategies can be extracted by running the two following queries: Then, under these two strategies, the expected value of the minimum and maximum energy balance for a given number of runs N are obtained using the following queries: • E[< = H; N] (min:kWh) under minkWh.  The resulting probability distribution constitutes the bounds of a probabilistic FlexOffers slice of duration H.
Finally, UPPAAL simulation output is parsed by the means of a regular expression, and the following insights are extracted and plotted for both the maximum and minimum consumption strategies:

Case study
The case study consists of a 150 m 2 4-room residential building with floor heating. Water is heated by a 2.5 kW heat pump. The floor plan of the building is depicted in Fig. 5. To extract the parameters of building thermal dynamics, a dataset containing 50 days of   (Vinther et al. 2017).
To run the simulation, the following assumptions are taken into consideration: • In a real implementation scenario, the first assumption is not necessary because we would have access to the actual state of the system at the beginning of the simulation. On the contrary, in this case, we need to hypothesize the initial state of the building to predict the temperature evolutions and generate a heat pump power schedule that satisfies the imposed temperature and consumption constraints.
The heat power of the four rooms is within the intervals of [0, 3500] W, [0, 3000] W, [0, 1000] W, and [0, 1200] W for rooms 1 to 4 respectively. To determine the COP, more input data about the heat exchanger are required which is out of the scope of this study. In this specific scenario, its value can be hypothesized to be around 3.7.
First, the historical input data are imported to the CTSM-R. The CTSM estimates the thermal dynamics of the simulated buildings. This way, 5 days' worth of sensor data is used to obtain the maximum likelihood estimation of thermal dynamic parameters (see Eq. (6), Fig. 6 Building historical data; indoor and outdoor temperatures over two days (7)). The data is equal to 7200 min (on a minute basis). Figures 6 describes the indoor temperature of rooms and ambient temperature over a two-days horizon.
Once fed with the provided dataset, CTSM estimates the thermal coefficients and converges to a reliable solution in approximately 5-10 min. The computation time varies based on the size and characteristics of rooms. For example, room 1 has the highest computation time due to large windows and envelopes. In contrast, room 3 has the lowest computation time due to small envelopes and window dimensions.
Given that this process needs to be performed only once per season and not in real-time, the computation time and resources required are more than sustainable both from computational power and energy consumption points of view. Without loss of generality, the specific ODEs characterizing the understudy building are formulated in Eqs. (11)-(13). The estimated thermal coefficients for each room, extracted from CTSM-R, are listed in Table 4.
It is worth mentioning that to increase the estimation accuracy of thermal dynamic estimation, the parameters Ω and Ψ are retrieved from the simple multiplication/division of heat resistance and capacity, i.e., R and C parameters of Eqs. (1-3). Therefore, the parameters Ω and Ψ are thermal coefficients and have a different meaning from the R and C values. This simplification is obtained based on "try and error" to increase the estimation accuracy and decrease the computation burden of the thermal dynamic estimation approach in R language.
The variable Q, i.e., heat consumption of the heating system, is the control variable. This is the decision variable to optimize the operation of the heating system generating FlexOffers. The heat consumption is a function of mass flow and inlet/outlet mass temperature as follows: where c ρ is the specific heat of water; m is mass flow; T in and T out are inlet and outlet temperature of the water, respectively. This variable can be controlled by the means of the four valves installed in each of the four rooms.
The objective is to compute the available consumption flexibility, i.e., generate Probabilistic FlexOffers. To achieve the aim, the power consumption of the heat pump is optimized while meeting the indoor temperature within residents' comfort bound. The residents' comfort bound is comprised of lower and upper thresholds as follows: The building thermal model described beforehand can be represented in UPPAAL-STRATEGO through the network of timed automata depicted in Fig. 7. In the first transition, from location Choose_Power to Wait, the position for each of the four valves is selected, i.e., open or closed. If it is set to open, the maximum amount of thermal power is applied to the r-th room for the whole duration of the next Flex-Offers slice. In the second transition, from location Wait to Choose_Power, time is allowed to pass, and function apply_flow() is used to: Update environment information, i.e., ambient temperature and solar irradiance, resp. T am and P s .
Update temperatures for each room, i.e., indoor, heater, and envelope temperatures, resp. T i r , T h r , and T e r , according to input variable using ODEs formulated in Eqs. (1-3). Update overall energy consumption, i.e., variable kWh. This cycle is repeated H times, where H represents the FlexOffers time horizon, i.e., the number of FlexOffer timeslots in the future for which we want to generate flexibility offers and heat pump schedules. By simulating the system evolution using the UPPAAL queries introduced beforehand, it is possible to produce the results which are discussed in the following section.

Implementation results
In this section, the simulation results of the case study are presented and discussed. Figures 8 and 9 depict the indoor temperature evolutions and the optimized thermal power schedule in the two different scenarios including the minimum consumption and the maximum consumption, respectively. In both cases, the plotted results are an average of 1000 simulations where different weather conditions have been taken into consideration. The whole simulation and scheduling process took about 160 s, on a standard Linux virtual machine equipped with an 8-core CPU @3.0 GHz and 16 GB of RAM. As a result, this approach allows us to schedule the heat pump operation for the next 24 h with an acceptable execution time.
Regarding the thermal power schedule, the bar graph represents the cumulative average heating power provided by the floor heating system in the associated Flex-Offer slice of 15 min. These bars also show how power is provided by each of the four-room heating systems. Moreover, the overall energy consumption evolution is depicted for the day. The final energy consumptions are 72.30 kWh and 112.31 kWh, for the minimum and maximum consumption scenarios, respectively. Therefore, it can provide 40.01 kWh of thermal energy flexibility for the upstream energy network. Fig. 8 The scenario of minimum energy consumption in terms of rooms' temperature, power scheduling, and cumulative energy consumption In the scenario of minimum energy consumption, it is shown that the indoor temperatures are kept as close as possible to the lower comfort threshold; consequently, the heating power is only used to maintain the indoor temperatures close to T i min threshold.
On the contrary, in the maximum consumption scenario, temperatures tend to reach and maintain a value around the T i max threshold.
We need to remark that the power and energy values under discussion are referred to as the thermal power (P h ) provided by the underfloor heating system, and it is different from the instant electrical power (P e ) consumed by the heat pump. The latter can be retrieved by dividing P h by the heat pump COP (i.e., Coefficient Of Performance). Considering a COP value of 3.7, we can convert the previously obtained thermal energy flexibility into electrical energy flexibility. This way, the building can offer 10.81 kWh of demand flexibility to the power grid. Table 5 elaborates on the electrical energy flexibility in different scenarios.
In this table, optimal consumption refers to the scenario where the heat pump maintains the rooms' temperatures as close as possible to the residents' setpoint. If we consider the final energy consumption value as the average of the maximum and minimum consumption scenarios, then the heat pump provides flexibility corresponding to 43.3% of the ideal building energy consumption. Fig. 9 The scenario of maximum energy consumption in terms of rooms' temperature, power scheduling, and cumulative energy consumption Finally, in Fig. 10, the generated Probabilistic FlexOffers are described. Blue and red lines represent the PDF (Probability Distribution Function) of energy consumption in the case of minimum and maximum consumption scenarios, respectively. Whereas the green line shows, for each energy consumption level, which is the likelihood that the assigned schedule can be followed. The two FlexOffers bounds are described by two Gaussian distributions including N (72.30, 0.66) for the lower bound and N (112.31, 0.63) for the upper bound. It can be observed that the variance of both normal distributions is relatively small despite the variance imposed on the weather information. This relies on the fact that, in this specific case, the building energy consumption is not highly affected by the weather conditions; that may be because of the high-quality thermal insulation which reduces the heat exchange between indoor and outdoor environments. Also, the unavailability of on-site renewable energy generation, e.g., roof-top photovoltaic sites, decreases the dependency on weather conditions. On the contrary, if the building is equipped, for instance, with a photovoltaic (PV) system, then the interdependency between the household energy profile and the solar irradiation increases. This way, an increase in PV power generation will decrease the energy demand from the power grid. Also, in some cases, it is possible to inject surplus PV generation into the power grid. Consequently, the flexibility potentials of the building increase significantly not only to operate the flexible heat pumps but also to integrate power flexibility of self-generation facilities into the upstream network. Adversely, a cloudier day than the forecasted one will reduce the PV energy generation and increase the energy supply from the power grid. Therefore, it will be interesting in future studies to investigate how Probabilistic FlexOffers variance is affected by weather conditions in buildings that are equipped with renewable self-generation units. As stated before, the 95% confidence interval on the generated FlexOffers is very valuable information for potential buyers, because it allows them to quantify the risk-bearing capacity. If the energy consumption schedule assigned to the consumer is included in the interval [73.39, 111.27] kWh, it will be successfully followed with a likelihood greater than 0.95; therefore, it is almost guaranteed that the assigned consumption pattern will be respected.
In Table 6, other confidence intervals with associated success probability are listed. In particular, an energy schedule from the 0.1% probability success interval would be chosen only if the power grid encounters a severe power deficit or excess. Therefore, the maximum amount of demand flexibility is required to provide up-/down-regulation in the opposite direction of power system imbalance.

Conclusion and future work
This paper discusses how the residential heat pumps can be operated to provide demand flexibility for the power grid. The approach addresses a probabilistic success function to minimize and maximize the energy consumption of the heat pump in response to the deficit and surplus of renewable power generation on the supply side. The residents' comfort bound is defined as a setpoint with lower and upper thresholds.
To examine the suggested approach, a high-fidelity building model is chosen under the Danish building regulations from 2010. The simulation results show that the FlexOffers approach minimizes energy consumption when a power shortage occurs on the power grid. This way, the indoor temperature is scheduled close to the lower threshold. On the contrary, when the power system encounters an excess of renewable power, the FlexOffers approach maximizes the energy consumption of the heat pump while maintaining the indoor temperature around the upper threshold.
Also, different scenarios are presented to address the uncertain weather forecast. Therefore, the approach generates the FlexOffers considering uncertainties associated with ambient temperature and solar power. Despite the fact that the weather conditions affect the heat pump operation, self-generation facilities, e.g., roof-top photovoltaic sites, can increase the flexibility potential of the residential buildings. To sum up, the results confirm that the suggested FlexOffers approach can provide flexibility on the opposite side of power system imbalance when a deficit/excess of renewable power occurs.
We plan to develop more coarse-grained thermal building models because often only a single indoor temperature is available. As part of this, we will look for alternatives for estimating the model parameters.