Model Predictive Control Strategies to Activate the Energy Flexibility for Zones with Hydronic Radiant Systems

: This paper presents control strategies to activate energy ﬂexibility for zones with radiant heating systems in response to changes in electricity prices. The focus is on zones with radiant ﬂoor heating systems for which the hydronic pipes are located deep in the concrete and, therefore, there is a signiﬁcant thermal lag. A perimeter zone test-room equipped with a hydronic radiant ﬂoor system in an environmental chamber is used as a case study. A low order thermal network model for the perimeter zone, validated with experimental measurements, is utilized to study various control strategies in response to changes in the electrical grid price signal, including short term (nearly reactive) changes of the order of 10–15 min notice. An index is utilized to quantify the building energy ﬂexibility with the focus on peak demand reduction for speciﬁc periods of time when the electricity prices are higher than usual. It is shown that the developed control strategies can aid greatly in enhancing the zone energy ﬂexibility and minimizing the cost of electricity and up to 100% reduction in peak power demand and energy consumption is attained during the high-price and peak-demand periods, while maintaining acceptable comfort conditions.


Introduction
The world-wide demand for electricity is expected to double by 2050 [1] and the electricity demand in the building sector is projected to increase by 70% [2]. Demand response (DR) management of buildings has a significant role in the interaction between the supply side and demand side. Demand response involves modifying the client's electricity consumption through incentive-based dynamic pricing techniques [3][4][5][6]. Due to recent global efforts on reducing greenhouse gas emissions, world-wide adoption and integration of renewable energy sources (RES) has been increasing at an accelerated pace. For example, the total installed capacity for PV at the end of 2019 reached 627 GW with 115 GW of PV installed in 2019 only [7] and the total global installed capacity for wind reached 651 GW with 60.4 GW installed in 2019 [8]. The trend is anticipated to go on, since several countries have ambitious targets and are thriving to increase the RES share in energy system. The Renewable Energy Directive requires the European Union to satisfy minimum 32% of its total energy demands with renewables by 2030 (https://eur-lex.europa.eu/eli/dir/2018/2001/oj (accessed on 20 January 2021)), while Germany aiming to fulfil minimum 80% of its overall electricity demand via renewable energies by 2050 [9]. Some countries like Denmark, plan to cover 100% of their energy demand with renewables [10]. In Quebec, commercial and institutional buildings represent approximately 18% of the total energy consumption, of which 37% is from natural gas and 57% from hydro-electricity [11]. In order to reduce GHG emissions and thanks to being surrounded by about 500,000 lakes and 4500 rivers ( https://www.hydroquebec.com/about/our-energy.html (accessed on 20 January 2021)), a shift from natural gas to hydroelectricity as the major source of energy has been observed in recent years. Consequently, the electricity demand has been increasing significantly, notably during the heating period in winter and diurnal demand peaks are spotted in the morning and evening which are putting much pressure on the electrical grid to satisfy the peak power demand due to space heating on extreme cold days. Thus, more and more RES including solar panels and wind turbines are being considered for integration into the grid.
However, incorporating the highly variable RES into existing grids is a significant challenge. A consequence of electricity production by renewables on electricity prices is that in areas with a relatively high expected output from renewables, there is a much higher tendency to have more fluctuating prices than areas with lower renewables penetration [12]. In addition to the centralized plants, RES are also expanding in a distributed manner, as well as distributed storage systems (such as thermal storage systems, battery storage, etc.). The flexibility to cope with the discrepancy between consumption and generation can come from either the supply side (through the use of conventional power plants or storage) [13][14][15][16] or from the demand side [17][18][19][20]. The energy flexible buildings (EFB) concept was introduced by IEA EBC Annex 67 [18]. The building energy flexibility as defined by annex 67 is "the ability to manage building demand and generation according to local climate conditions, user needs, and energy network requirements" [21]. Characterizing buildings energy flexibility is challenging since it involves various aspects, and it is strongly dependent on the boundary conditions and assessment methods [19]. Hence, energy flexibility quantification in buildings is one of the high-interest topics in current research [22][23][24][25][26].
Model predictive control is one of the key approaches in activating the energy flexibility that can be offered by the building especially when considering the storage in building elements and thermal lag to meet smart grid needs. MPC utilizes a mathematical model of the building to predict its future behavior. Considering these predictions, MPC determines the optimal control actions based on a defined objective while considering thermal comfort and other considered constraints, and weather forecasts in a systematic and flexible way. Although there is an abundance of research work and pilot projects for MPC [27][28][29][30], the transfer of this technology to the building market is still in its early stages. One challenge in the building sector is that MPC is still relatively new and the sector generally considers satisfactory comfort conditions in buildings a priority [31]. Moreover, models are backbone of MPC and uniqueness of every building that requires customized modeling which is another major challenge that increases the engineering time and cost when applying model-based control strategies [32]. Therefore, there is a preference towards simple and low-order models when developing models for MPC. Low-order control-oriented models have several advantages compared to detailed models [33]. They are faster to create with limited available information about the building [34]. They are easier to calibrate and model parameters can be adjusted in real-time. The model resolution can be chosen depending on the objectives, control variables, etc. A methodology to develop a low-order control-oriented model for zones with radiant floor heating system is presented in the following sections.
Radiant floor heating systems have been receiving considerable attention due to the multiple advantages they offer such as improved thermal comfort in buildings and suitability for other related applications in cold climates. Hydronic radiant heating can utilize low temperature renewable energy heat sources including geothermal or solarsource heat pumps (water temperature as low as 35 • C) and provide significant flexibility to smart grids by storing energy and shifting heating demand [26]. In both cases heat pumps operate with electricity. The operation of these systems can be optimized by applying predictive control and further the energy costs can be reduced by optimizing their interaction with smart grids by utilizing the flexibility in their demand profiles. Due to the large thermal mass of the concrete floor, there is time lag between the heat supply of the radiant floor and the indoor temperature response, which could be up to several hours, depending on the pipe depth and thermal capacitance of the concrete floor. This inherent delay makes radiant floor heating system a suitable candidate for MPC. The focus of this paper is on zones equipped with high mass radiant systems (pipes buried deep in concrete) and such case study is used which is introduced in the following sections.

Research Gap and Contribution
There have been several recent studies on different flexibility services that buildings can provide such as reducing peak loads and shifting peak demand through utilization of thermal mass [35][36][37][38], storage in batteries [39], charging of electric vehicles [40], and scheduling of the HVAC system utilization [21]. However, the current building energy flexibility literature is limited to few studies on control strategies that considers the electricity price prediction of the following days [41,42] and specifically, there is a limited number of such studies for buildings equipped with radiant floor heating systems. Additionally, no previous study was found on the short-notice changes in price signal and the strategies to deal with such situations in buildings. This is an important scenario that occurs when power grid faces unexpected load balances and buildings as the largest electricity consumers play a significant role in those situations. Hence, utilities such as Hydro-Quebec are looking for solutions and strategies to deal with such critical situations. This article is contributing towards increasing knowledge and presenting energy flexibility strategies for buildings with radiant floor heating systems.

Building Energy Flexibility Index (BEFI)
A quantitative assessment of the energy flexibility provided by structural thermal energy storage is essential to large scale deployment of thermal mass as in an active demand response (ADR) context. The available storage capacity expresses the amount of energy that can be added to the structure's thermal energy storage (STES) during a specific event. Therefore, the heat that can be stored within a dwelling not only depends upon the thermal properties of the building fabric, but also on the properties and actual use of the heating and ventilation systems [43]. Four performance indicators or characteristics for ADR are defined and quantification methods for the ADR potential of structural thermal storage are presented by Reynders et al. [22] that are mainly focused on the energy that is reduced/increased over a certain period. Here, an index is applied which is explicitly based on power (Watts) reduction/increase. A building energy flexibility index (BEFI) is defined by Athienitis et al. [44] as the average amount of power (kW) that can be increased or decreased relative to a reference power profile (Q ref , the "as-usual operation") by the building during a certain period of time (Dt): The schematic below in Figure 1 illustrates the BEFI concept. The grey area during the period Dt is the numerator of Equation (1) and shows the difference between the flexible and the reference heating loads.
The BEFI could be predicted continuously with a certain accuracy/uncertainty depending on the configuration and duration length. Generally, the shorter the duration, the less the uncertainty is expected. Another factor to consider is the time of activation. For some grid needs, activation time must be within seconds and must be totally autonomous. For others, it can be different and may rely on communication protocols, but once the signal is received by the building, BAS response must be automated. Then, for some others, the need may be planned half a day in advance and the building can be designed and operated to maximize its flexibility during the period of need with potential significant financial benefits in terms of reduced operational costs. In this paper, the BEFI is calculated with dynamic tariffs and a floor heating system.  The BEFI could be predicted continuously with a certain accuracy/uncertainty depending on the configuration and duration length. Generally, the shorter the duration, the less the uncertainty is expected. Another factor to consider is the time of activation. For some grid needs, activation time must be within seconds and must be totally autonomous. For others, it can be different and may rely on communication protocols, but once the signal is received by the building, BAS response must be automated. Then, for some others, the need may be planned half a day in advance and the building can be designed and operated to maximize its flexibility during the period of need with potential significant financial benefits in terms of reduced operational costs. In this paper, the BEFI is calculated with dynamic tariffs and a floor heating system.

Modeling of Thermal Zones with Radiant Floor Systems
The modeling approach used for this study is based on the low-order lumped-parameter thermal network models which have been shown to be practical for control studies and especially model predictive control [33,45]. This approach is generally valid when Biot number is less than 0.1. In this approach the thermal mass under consideration (radiant floor concrete slab) is discretized into a number of control volumes. Each of the discretized control volumes is represented by a node and considered to be isothermal. Each of the nodes, has a lumped thermal capacitance connected to it and thermal conductances connecting it to adjacent nodes. Considering the time interval p and time step Δt, the general form of the explicit finite-difference model for the nodes with a lumped thermal capacitance can be stated as where Ti,p represents the temperature of node "i" at time step "p", Tj,p represents temperature of node "j" at time step "p", Ci is the thermal capacitance of node "i", Ri,j is the thermal resistance between nodes i and j and Qi is the heat source at node i. A model with fewer parameters facilitates setting up the initial conditions which is a key parameter for control studies [33]. When the details of the construction are not known, a low order grey box modelling approach is practical and can be developed and calibrated

Modeling of Thermal Zones with Radiant Floor Systems
The modeling approach used for this study is based on the low-order lumpedparameter thermal network models which have been shown to be practical for control studies and especially model predictive control [33,45]. This approach is generally valid when Biot number is less than 0.1. In this approach the thermal mass under consideration (radiant floor concrete slab) is discretized into a number of control volumes. Each of the discretized control volumes is represented by a node and considered to be isothermal. Each of the nodes, has a lumped thermal capacitance connected to it and thermal conductances connecting it to adjacent nodes. Considering the time interval p and time step ∆t, the general form of the explicit finite-difference model for the nodes with a lumped thermal capacitance can be stated as where T i,p represents the temperature of node "i" at time step "p", T j,p represents temperature of node "j" at time step "p", C i is the thermal capacitance of node "i", R i,j is the thermal resistance between nodes i and j and Q i is the heat source at node i. A model with fewer parameters facilitates setting up the initial conditions which is a key parameter for control studies [33]. When the details of the construction are not known, a low order grey box modelling approach is practical and can be developed and calibrated by means of real time data from the building. The models must be accurate enough to provide reliable information but also flexible enough for quick and computationally efficient decision-making [46,47], especially in reaction to electric grid's short notices on the change of the price signals.
An important part of a model for zones with radiant floor heating system is the radiative and convective heat transfer which are naturally nonlinear processes. However, the respective heat transfer coefficients are typically linearized so that the system energy balance equations can be solved with linear algebra techniques and represented with a linear thermal network. In the case of radiant floor heating, this linearization generally introduces less error for the for long-wave radiation heat transfer (h r ) than the convection heat exchange between the radiant floor surface and room air (h cf ) [48]. For example, in the case of radiant floor heating where usually the floor is hotter than the air and the heat flow is upward h cf is in the order of 3 W/(m 2 K) while for a cold floor and warmer air it is in the order of 1 W/(m 2 K). For heat flow upwards, h cf can be calculated as [49]: where T floor is the floor surface temperature and T air is the zone air temperature. Therefore, usually certain amount of calibration for the convection heat exchange between the radiant floor and room air is required for a model especially when considering the low order models.
It is demonstrated in the following section, that a well calibrated second-order model, for which the thermal network is shown in Figure 2, can accurately capture the most important thermal dynamics of a zone with a radiant floor heating system. Figure 2 also shows the boundary conditions for the model. Tg is the ground temperature, h c and h r are the convective and radiative heat transfer coefficients, respectively, and R inf is the infiltration thermal resistance.
cient decision-making [46,47], especially in reaction to electric grid's short notices on the change of the price signals.
An important part of a model for zones with radiant floor heating system is the radiative and convective heat transfer which are naturally nonlinear processes. However, the respective heat transfer coefficients are typically linearized so that the system energy balance equations can be solved with linear algebra techniques and represented with a linear thermal network. In the case of radiant floor heating, this linearization generally introduces less error for the for long-wave radiation heat transfer (hr) than the convection heat exchange between the radiant floor surface and room air (hcf) [48]. For example, in the case of radiant floor heating where usually the floor is hotter than the air and the heat flow is upward hcf is in the order of 3 W/(m 2 K) while for a cold floor and warmer air it is in the order of 1 W/(m 2 K). For heat flow upwards, hcf can be calculated as [49]: where Tfloor is the floor surface temperature and Tair is the zone air temperature. Therefore, usually certain amount of calibration for the convection heat exchange between the radiant floor and room air is required for a model especially when considering the low order models. It is demonstrated in the following section, that a well calibrated second-order model, for which the thermal network is shown in Figure 2, can accurately capture the most important thermal dynamics of a zone with a radiant floor heating system. Figure 2 also shows the boundary conditions for the model. Tg is the ground temperature, hc and hr are the convective and radiative heat transfer coefficients, respectively, and Rinf is the infiltration thermal resistance. In the above thermal network for the case study (Perimeter Zone Test Cell (PZTC)), the hydronic pipes are located at the bottom of the slab and there is a large thermal mass above the pipes that causes significant thermal lag. There can be cases where pipes are placed in the middle of the slab in which case the thermal dynamics will be different. In the above thermal network for the case study (Perimeter Zone Test Cell (PZTC)), the hydronic pipes are located at the bottom of the slab and there is a large thermal mass above the pipes that causes significant thermal lag. There can be cases where pipes are placed in the middle of the slab in which case the thermal dynamics will be different. There can also be a case where the location of the pipes is not exactly known, and data from the building automation system can be utilized to calibrate the model.

Case Study: Perimeter Zone Test Cell (PZTC)
The Solar Simulator-Environmental Chamber (SSEC) laboratory is an experimental facility located at Concordia University in downtown Montreal, Canada. This facility allows accurate and repeatable testing of solar systems and advanced building envelopes under standard test conditions with simulated solar radiation and indoor plus outdoor conditions. The temperature test range of the environmental chamber is −40 • C to +50 • C. Two solar simulators (large-scale and mobile) are designed to emulate solar radiation to test solar systems such as PV and PV/T modules, solar air collectors, solar water collectors and building-integrated solar systems.
The perimeter zone test cell (PZTC) (case study) is a 3 m × 3 m × 3 m room placed inside the solar simulator/environmental chamber laboratory and it contains a radiant floor heating/cooling system. The side and back walls and ceiling consists of 10 cm of insulation with thermal resistance of 5.64 m 2 K/W. The front wall is a photovoltaic/thermal  Table 1. The pipes of the radiant floor system are made of conventional cross-linked polyethylene (PEX) and have an external diameter of 1.75 cm. The pipes are installed in a foam matrix of insulating material that also facilitates keeping them in place. The pipes are spaced at 15 cm. Figure 3 shows the piping configuration of the radiant floor before the concrete was poured (left), during pouring of concrete and final look of the radiant floor slab (right).
facility located at Concordia University in downtown Montreal, Canada. This facility allows accurate and repeatable testing of solar systems and advanced building envelopes under standard test conditions with simulated solar radiation and indoor plus outdoor conditions. The temperature test range of the environmental chamber is −40 °C to +50 °C. Two solar simulators (large-scale and mobile) are designed to emulate solar radiation to test solar systems such as PV and PV/T modules, solar air collectors, solar water collectors and building-integrated solar systems.
The perimeter zone test cell (PZTC) (case study) is a 3 m × 3 m × 3 m room placed inside the solar simulator/environmental chamber laboratory and it contains a radiant floor heating/cooling system. The side and back walls and ceiling consists of 10 cm of insulation with thermal resistance of 5.64 m 2 K/W. The front wall is a photovoltaic/thermal solar test façade with 0.88 m 2 K/W thermal resistance which was not in operation during the experiment. The floor is made of an approximately 8 cm thick concrete slab, with the piping at the bottom of the concrete and insulation of 7.4 m 2 K/W under the slab. Thermal properties of the concrete are shown in Table 1. The pipes of the radiant floor system are made of conventional cross-linked polyethylene (PEX) and have an external diameter of 1.75 cm. The pipes are installed in a foam matrix of insulating material that also facilitates keeping them in place. The pipes are spaced at 15 cm. Figure 3 shows the piping configuration of the radiant floor before the concrete was poured (left), during pouring of concrete and final look of the radiant floor slab (right).  A second-order model, for which the thermal network was shown in Figure 2, was developed for the PZTC (Bi ≈ 0.065 < 0.1). The floor convective heat transfer can vary significantly depending on the temperature difference between the floor surface temperature  A second-order model, for which the thermal network was shown in Figure 2, was developed for the PZTC (Bi ≈ 0.065 < 0.1). The floor convective heat transfer can vary significantly depending on the temperature difference between the floor surface temperature and the air in the zone. Therefore, an objective function is defined to find the effective values of h cf so that the floor surface temperature calculated from the simulations (T floor, simulated ) matches with experiment (T floor, measured ) as accurately as possible. Therefore, the objective function is defined as Optimization was done in MATLAB using fmincon function which uses the simplex algorithm. The optimization result is shown in Figure 5. As observed, a well calibrated second-order model, for which the thermal network shown in Figure 2, can capture thermal dynamics of a zone that contains radiant floor heating system with a good accuracy and the maximum error between the model and the measurements is about 0.45 • C. However, for the majority of the time the model error is less than 0.1 • C. This validated model is used in the following sections to study optimal control strategies while dealing with different price signals from the grid. The strategies to deal with the various penalty signals presented in the following sections of this paper are divided into two main categories which are: 1. Predictive strategies with short notice in the order of 10 to 15 mins; 2. Predictive strategies with a prediction horizon of about 1-2 days.
It should be noted that in Quebec two peak demand periods are observed in a day during the heating season which are typically in the morning from 6 to 9 a.m. and in the evening from 4 to 8 p.m. Therefore, the price of the electricity can be different and higher during these periods compared to other times of the day as part of a potential dynamic tariff structure intended to motivate reduced demand for electricity during these periods. The price signals that are assumed and utilized in the following sections are created based on those typical peak demand periods in Quebec.

Control Strategies with a Short Notice
In this section thermal response of the radiant floor in reaction to a sudden change in the grid price signal is presented. The control of the radiant floor heating is based on the floor surface temperature since it provides a much smoother heating load profile due to the thermal mass of the floor compared to controlling the zone air temperature which has a very low thermal capacity and hence it causes lots of fluctuations in the heating load. The heating is controlled through proportional control as where Tsp = floor surface temperature setpoint, Tfloor = floor surface temperature, and Kp is Model error (°C) This validated model is used in the following sections to study optimal control strategies while dealing with different price signals from the grid. The strategies to deal with the various penalty signals presented in the following sections of this paper are divided into two main categories which are: 1. Predictive strategies with short notice in the order of 10 to 15 mins; 2. Predictive strategies with a prediction horizon of about 1-2 days.
It should be noted that in Quebec two peak demand periods are observed in a day during the heating season which are typically in the morning from 6 to 9 a.m. and in the evening from 4 to 8 p.m. Therefore, the price of the electricity can be different and higher during these periods compared to other times of the day as part of a potential dynamic tariff structure intended to motivate reduced demand for electricity during these periods. The price signals that are assumed and utilized in the following sections are created based on those typical peak demand periods in Quebec.

Control Strategies with a Short Notice
In this section thermal response of the radiant floor in reaction to a sudden change in the grid price signal is presented. The control of the radiant floor heating is based on the floor surface temperature since it provides a much smoother heating load profile due to the thermal mass of the floor compared to controlling the zone air temperature which has a very low thermal capacity and hence it causes lots of fluctuations in the heating load. The heating is controlled through proportional control as where T sp = floor surface temperature setpoint, T floor = floor surface temperature, and Kp is the proportional constant equal to 900 W/ • C. The maximum heating output is limited to the size of the system which is 2.7 kW. The ambient temperature profile considered is shown in Figure 6. The ambient temperature changes between −7 • C and −18 • C with the mean value of about −12 • C-a relatively cold day in Quebec when high demand for electric space heating is expected.  Figure 7 shows the heating load, floor surface temperature and zone air temperature for maintaining a constant floor surface temperature of 24 °C which is considered the reference surface temperature here. The results shown in Figure 5 is for a day when the electricity price remains constant for the whole day at the reference value which is assumed to be $10 per kW for power.   3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25 10 $/kW Figure 6. Ambient temperature profile. Figure 7 shows the heating load, floor surface temperature and zone air temperature for maintaining a constant floor surface temperature of 24 • C which is considered the reference surface temperature here. The results shown in Figure 5 is for a day when the electricity price remains constant for the whole day at the reference value which is assumed to be $10 per kW for power.
As can be observed from Figure 7, when the floor surface temperature is kept around 24 • C, the zone air temperature stays between 18 • C and 20 • C which is within the typical indoor air temperature comfort range (18-24 • C) for occupants [50] and the operative temperature also stays within the same range. Figure 8 shows the short-notice predictive (almost reactive) behavior to sudden increase, from 10 $/kW to 20 $/kW, in the price signal for a period of 3.5 h (hours between 30 and 33.5) which is the typical morning peak demand period in Quebec (6-9 a.m.). It is assumed that the price of energy (kWh) is constant and only the power (kW) price is changing. The simplest strategy to deal with this sudden increase in the price signal is to turn off the heating system. It can be observed that due to shutdown of the heating system during that the period, the floor surface temperature (red line) falls to 22.2 • C and air temperature (purple line) to 18 • C. Considering the typical comfort boundaries of 18-24 • C for the zone air temperature, the air temperature is at the lower threshold of the comfort boundary. It can be calculated that the BEFI would be roughly about 400 W for 3.5 h. The floor surface temperature is also shown for the case with non-linear convective heat transfer coefficient between floor surface and air temperature (h cf ). As observed, there is no significant difference between the two curves, thus confirming that assuming a linear convective coefficient is an acceptable approximation.  Figure 7 shows the heating load, floor surface temperature and zone air temperature for maintaining a constant floor surface temperature of 24 °C which is considered the reference surface temperature here. The results shown in Figure 5 is for a day when the electricity price remains constant for the whole day at the reference value which is assumed to be $10 per kW for power.   indoor air temperature comfort range (18-24 °C) for occupants [50] and the operative temperature also stays within the same range. Figure 8 shows the short-notice predictive (almost reactive) behavior to sudden increase, from 10 $/kW to 20 $/kW, in the price signal for a period of 3.5 h (hours between 30 and 33.5) which is the typical morning peak demand period in Quebec (6-9 a.m.). It is assumed that the price of energy (kWh) is constant and only the power (kW) price is changing. The simplest strategy to deal with this sudden increase in the price signal is to turn off the heating system. It can be observed that due to shutdown of the heating system during that the period, the floor surface temperature (red line) falls to 22.2 °C and air temperature (purple line) to 18 °C. Considering the typical comfort boundaries of 18-24 °C for the zone air temperature, the air temperature is at the lower threshold of the comfort boundary. It can be calculated that the BEFI would be roughly about 400 W for 3.5 h. The floor surface temperature is also shown for the case with non-linear convective heat transfer coefficient between floor surface and air temperature (hcf). As observed, there is no significant difference between the two curves, thus confirming that assuming a linear convective coefficient is an acceptable approximation.  Figure 9 shows another increase in the price signal during the evening on the same day. This second increase is during the evening peak demand hours (4-8 p.m.) in Quebec during winter. As observed, since the second increase event happens when the zone air temperature was higher compared to the first event and generally the system was at a higher state of charge, the air temperature drops to 18.5 °C at lowest and the building can offer the average BEFI of 400 W for 4 h during this period. The operative temperature stays between 19.5 °C and 21 °C.  Figure 9 shows another increase in the price signal during the evening on the same day. This second increase is during the evening peak demand hours (4-8 p.m.) in Quebec during winter. As observed, since the second increase event happens when the zone air temperature was higher compared to the first event and generally the system was at a higher state of charge, the air temperature drops to 18.5 • C at lowest and the building can offer the average BEFI of 400 W for 4 h during this period. The operative temperature stays between 19.5 • C and 21 • C. The intermittent nature of renewable energy sources increases the necessity for backup power to deal with the unexpected swings in power production. Therefore, change in the price signal can be more frequent and the price signal from the grid can vary every hour during the day depending on the availability of electricity from renewable energy sources, weather conditions, emergency events, etc. For example, the dynamic price signal for the following 48 h may look like the profile shown in Figure 10. As can be observed from Figure 10, the price signal is changing frequently during the day and if we consider the $10 per kW as the normal, reference operating price signal, there are some periods with price signals higher, as well as lower than the reference. The intermittent nature of renewable energy sources increases the necessity for backup power to deal with the unexpected swings in power production. Therefore, change in the price signal can be more frequent and the price signal from the grid can vary every hour during the day depending on the availability of electricity from renewable energy sources, weather conditions, emergency events, etc. For example, the dynamic price signal for the following 48 h may look like the profile shown in Figure 10. As can be observed from Figure 10, the price signal is changing frequently during the day and if we consider the $10 per kW as the normal, reference operating price signal, there are some periods with price signals higher, as well as lower than the reference. The intermittent nature of renewable energy sources increases the necessity for backup power to deal with the unexpected swings in power production. Therefore, change in the price signal can be more frequent and the price signal from the grid can vary every hour during the day depending on the availability of electricity from renewable energy sources, weather conditions, emergency events, etc. For example, the dynamic price signal for the following 48 h may look like the profile shown in Figure 10. As can be observed from Figure 10, the price signal is changing frequently during the day and if we consider the $10 per kW as the normal, reference operating price signal, there are some periods with price signals higher, as well as lower than the reference. It should be noted that a short notice (10-15 min or shorter) of a price change from the grid can also be received. The following control strategy is presented for these situations. First let us categorize the upcoming hour price signal into two categories:

1.
Price signal lower than the normal reference operating price signal that is: P(t) < P reference In this case since the price signal is lower than reference then we can just keep the heating control as usual and operating in the reference case (Q aux,ref ) for which the details were explained in equation (5). Therefore, the heating during the low-price period, Q aux,low-price , is inserted as 2. Price signals higher than the normal reference operating price signal: P(t) > P reference In this case the strategy will be to keep the floor surface temperature at 24 • C; however, the heating load will be a fraction of the reference heating and this fraction gets smaller as the price increases. Therefore, the strategy is defined as where P reference is the reference price and P(t) is the price at time t. Figure 11 shows the application of this strategy to the thermal zone: It should be noted that a short notice (10-15 min or shorter) of a price change from the grid can also be received. The following control strategy is presented for these situations. First let us categorize the upcoming hour price signal into two categories: 1. Price signal lower than the normal reference operating price signal that is: In this case since the price signal is lower than reference then we can just keep the heating control as usual and operating in the reference case (Qaux,ref) for which the details were explained in equation (5). Therefore, the heating during the low-price period, Qaux,low-price, is inserted as , , 2. Price signals higher than the normal reference operating price signal: In this case the strategy will be to keep the floor surface temperature at 24 °C; however, the heating load will be a fraction of the reference heating and this fraction gets smaller as the price increases. Therefore, the strategy is defined as where Preference is the reference price and P(t) is the price at time t. Figure 11 shows the application of this strategy to the thermal zone: Figure 11. Application of the nearly reactive strategy to deal with the changing price signal.
As can be observed, during the high price signal periods, whenever necessary, the heating load was turned on as a fraction of the reference heating load to keep the floor surface temperature near the setpoint hence reducing the amount of heating power and reducing the cost. Figure 12 shows the comparison between the reference heating load and the flexible heating load from the strategy above. It can be observed how the above strategy decreases the heating load during the high price periods. As mentioned earlier, BEFI is the amount of power (W) that the building can reduce (or increase) during certain periods of time. Therefore, BEFI during high price signal is the difference between the two curves which is shown by the area covered with the grey arrows.   Figure 11. Application of the nearly reactive strategy to deal with the changing price signal.
As can be observed, during the high price signal periods, whenever necessary, the heating load was turned on as a fraction of the reference heating load to keep the floor surface temperature near the setpoint hence reducing the amount of heating power and reducing the cost. Figure 12 shows the comparison between the reference heating load and the flexible heating load from the strategy above. It can be observed how the above strategy decreases the heating load during the high price periods. As mentioned earlier, BEFI is the amount of power (W) that the building can reduce (or increase) during certain periods of time. Therefore, BEFI during high price signal is the difference between the two curves which is shown by the area covered with the grey arrows.

Predictive control strategies with longer prediction horizons
When the prediction of the upcoming days for weather as well the grid price signals are available, the system can be operated in order to minimize the power and/or energy consumption during the high price signal period. In this section the predictive strategies for the price signals given in the previous sections are presented.
Here, the aim is to find a setpoint for the floor surface temperature proportional controller so that the surface temperature stays within the comfort range with 29 °C as maximum based on the ASHRAE standard [50]. The objective function is defined as minimize so that : 23 C 29 aux floor J Psig Q T C = × →°≤ ≤° (8) where Psig is the value of the price signal at every time step. Considering the above range for the floor surface temperature, the operative temperature (Top) stays between 19 °C and 25 °C. Optimization was done in MATLAB using fmincon function which uses the simplex algorithm. Figure 13 shows the results for the first price signal. As observed, the optimized floor surface temperature setpoint modulates the heating load of the floor in a way to keep the floor surface temperature very close to the minimum of 23 °C during the periods with reference price signal but as we get closer to the high price signal period, it starts preheating the floor slab by increasing the setpoint as high as 26 °C. Therefore, as observed, the heating load is zero during the high price period. The cost for whole period is $18.60. Now if no optimization was performed and the reference strategy (Figure 7) was used without considering the change in the price, then the total cost would be $20.70. Therefore, a decrease of more than 10% in the cost is observed.

Predictive control strategies with longer prediction horizons
When the prediction of the upcoming days for weather as well the grid price signals are available, the system can be operated in order to minimize the power and/or energy consumption during the high price signal period. In this section the predictive strategies for the price signals given in the previous sections are presented.
Here, the aim is to find a setpoint for the floor surface temperature proportional controller so that the surface temperature stays within the comfort range with 29 • C as maximum based on the ASHRAE standard [50]. The objective function is defined as where Psig is the value of the price signal at every time step. Considering the above range for the floor surface temperature, the operative temperature (T op ) stays between 19 • C and 25 • C. Optimization was done in MATLAB using fmincon function which uses the simplex algorithm. Figure 13 shows the results for the first price signal. As observed, the optimized floor surface temperature setpoint modulates the heating load of the floor in a way to keep the floor surface temperature very close to the minimum of 23 • C during the periods with reference price signal but as we get closer to the high price signal period, it starts preheating the floor slab by increasing the setpoint as high as 26 • C. Therefore, as observed, the heating load is zero during the high price period. The cost for whole period is $18.60. Now if no optimization was performed and the reference strategy (Figure 7) was used without considering the change in the price, then the total cost would be $20.70. Therefore, a decrease of more than 10% in the cost is observed.
Next, if another price signal was in the forecast in which two periods of high price are expected during the next 48 h, then a new setpoint found from the optimization is applied as shown in Figure 14. Next, if another price signal was in the forecast in which two periods of high price are expected during the next 48 h, then a new setpoint found from the optimization is applied as shown in Figure 14. As can be observed from both Figure 13 and Figure 14, the setpoint performs in such a way that the zone is heated up until just before the high price periods and then the setpoint is relatively low during that period so that due the thermal mass of the radiant floor concrete slab, the surface temperature does not drop below that setpoint and therefore, no heating is needed. The cost for the whole period is about $19, and compared to using reference strategy (cost = $22), a decrease of 13.63% in the cost is observed. Next, if another price signal was in the forecast in which two periods of high price are expected during the next 48 h, then a new setpoint found from the optimization is applied as shown in Figure 14. As can be observed from both Figure 13 and Figure 14, the setpoint performs in such a way that the zone is heated up until just before the high price periods and then the setpoint is relatively low during that period so that due the thermal mass of the radiant floor concrete slab, the surface temperature does not drop below that setpoint and therefore, no heating is needed. The cost for the whole period is about $19, and compared to using reference strategy (cost = $22), a decrease of 13.63% in the cost is observed. As can be observed from both Figures 13 and 14, the setpoint performs in such a way that the zone is heated up until just before the high price periods and then the setpoint is relatively low during that period so that due the thermal mass of the radiant floor concrete slab, the surface temperature does not drop below that setpoint and therefore, no heating is needed. The cost for the whole period is about $19, and compared to using reference strategy (cost = $22), a decrease of 13.63% in the cost is observed.
There can be more complicated price signals in the forecast as well. Figure 15 shows the predictive strategy for the price signal that was introduced in the previous section which changes almost every hour. Here, this price signal is known the day before from the forecast. There can be more complicated price signals in the forecast as well. Figure 15 shows the predictive strategy for the price signal that was introduced in the previous section which changes almost every hour. Here, this price signal is known the day before from the forecast. As can be observed from Figure 15, the building is heated up during the low-price period (until around 6 a.m.) and the setpoint is increased to the maximum upper threshold for floor surface temperature which is 29 °C. Then, after 6 a.m., there is no heating during the peak demand period (6-9 a.m.) and the warm slab thermal mass lets the heating be turned off until hour 13. Then, the setpoint is again increased to 29 °C to preheat the building and in order to minimize the heating during the afternoon peak demand hours (4-8 p.m.). A similar trend is observed during the next day as well. The total cost with the optimized setpoint is $18.80 compared to $24.12 with the reference strategy, resulting in a decrease of 22% in cost.

Emergency and Non-Predicted Change in the Price Signal
Contingency or operating reserve is the amount of power that utility may call when needed to face the loss of a generation unit or other unexpected load balance [51]. This reserve needs to be available typically within a ten-minute notification time, for a duration of one hour. Therefore, the building automation system or smart thermostat should be able to quantify the amount of power that it can decrease from its current consumption profile for the following hour. Figures 11-13 demonstrate the strategy when prediction of the weather as well as the price signal is available. Figure 16 shows the BEFI during a contingency/emergency event when the original predictive strategy before this event was based on the price signal #2 ( Figure 14).
As can be observed from Figure 15, the building is heated up during the low-price period (until around 6 a.m.) and the setpoint is increased to the maximum upper threshold for floor surface temperature which is 29 • C. Then, after 6 a.m., there is no heating during the peak demand period (6-9 a.m.) and the warm slab thermal mass lets the heating be turned off until hour 13. Then, the setpoint is again increased to 29 • C to preheat the building and in order to minimize the heating during the afternoon peak demand hours (4-8 p.m.). A similar trend is observed during the next day as well. The total cost with the optimized setpoint is $18.80 compared to $24.12 with the reference strategy, resulting in a decrease of 22% in cost.

Emergency and Non-Predicted Change in the Price Signal
Contingency or operating reserve is the amount of power that utility may call when needed to face the loss of a generation unit or other unexpected load balance [51]. This reserve needs to be available typically within a ten-minute notification time, for a duration of one hour. Therefore, the building automation system or smart thermostat should be able to quantify the amount of power that it can decrease from its current consumption profile for the following hour. Figures 11-13 demonstrate the strategy when prediction of the weather as well as the price signal is available. Figure 16 shows the BEFI during a contingency/emergency event when the original predictive strategy before this event was based on the price signal #2 ( Figure 14).
As observed, there is a sudden unexpected increase in the price signal between hours 33.5 and 35.5 (9:30-11:30 a.m.), and therefore the heating was turned off immediately. As shown in Figure 16 turning off the heating during that period has little impact on the floor surface and air temperatures. The BEFI which is the amount of power that is reduced compared to the reference case ( Figure 14) and is calculated on average as, BEFI = 270 W.  As observed, there is a sudden unexpected increase in the price signal between hours 33.5 and 35.5 (9:30-11:30 a.m.), and therefore the heating was turned off immediately. As shown in Figure 16 turning off the heating during that period has little impact on the floor surface and air temperatures. The BEFI which is the amount of power that is reduced compared to the reference case ( Figure 14) and is calculated on average as, BEFI = 270 W.
The timing of the unexpected change in the price signal is also particularly important. For example, if the unexpected change happens in a different hour that is closer to the period with a high price, then it could be more challenging to deal with this sudden change. As shown in Figure 17, if the sudden increase in price signal happens between hours 38.5 and 39.5, then the controller which is acting based on the predictions would not have enough time to follow back the optimized setpoint profile and as the result, some heating is observed during the predicted high price signal period as shown in Figure 17 with black dashed circles. However, as can be observed from Figure 17, since the floor surface temperature is already above the lower threshold of 23 °C, there is no need for heating. Therefore, in this case the setpoint temperature can be modified for the period between 40 and 44 h and put at the lower threshold of 23 °C. This modification will enable zero power consumption during that period as shown in Figure 18. The timing of the unexpected change in the price signal is also particularly important. For example, if the unexpected change happens in a different hour that is closer to the period with a high price, then it could be more challenging to deal with this sudden change. As shown in Figure 17, if the sudden increase in price signal happens between hours 38.5 and 39.5, then the controller which is acting based on the predictions would not have enough time to follow back the optimized setpoint profile and as the result, some heating is observed during the predicted high price signal period as shown in Figure 17 with black dashed circles. However, as can be observed from Figure 17, since the floor surface temperature is already above the lower threshold of 23 • C, there is no need for heating. Therefore, in this case the setpoint temperature can be modified for the period between 40 and 44 h and put at the lower threshold of 23 • C. This modification will enable zero power consumption during that period as shown in Figure 18.

Dynamic Pricing of Energy
There are certain situations in which the dynamic pricing is explicitly mentioned for energy consumption (kWh) instead of power (kW). In this case the optimization must be done for energy and the cost for energy (kWh) is minimized instead of power (kW). The unit of the price signal will be cents per kWh (ȼ/kWh). Figure 19 shows the results.

Dynamic Pricing of Energy
There are certain situations in which the dynamic pricing is explicitly mentioned for energy consumption (kWh) instead of power (kW). In this case the optimization must be done for energy and the cost for energy (kWh) is minimized instead of power (kW). The unit of the price signal will be cents per kWh (¢/kWh). Figure 19 shows the results.  As observed from Figure 19, the heating works in a similar way to when it was used to minimize the power consumption. The moving mean of the setpoint profile obtained from the optimization is used since the optimized setpoint contains lots of spikes in the temperature which consequently causes spikes in the heating loads and much on/off cycling which is not practical. Therefore, as can be observed, when minimizing the cost of power, the cost of energy is also minimized.

Conclusions
The major contributions of this article can be listed as 1. a methodology to develop low-order models for zones with high-mass radiant floor heating systems; 2. a method for short notice predictive control in response to change in the grid price signal for zones with As observed from Figure 19, the heating works in a similar way to when it was used to minimize the power consumption. The moving mean of the setpoint profile obtained from the optimization is used since the optimized setpoint contains lots of spikes in the temperature which consequently causes spikes in the heating loads and much on/off cycling which is not practical. Therefore, as can be observed, when minimizing the cost of power, the cost of energy is also minimized.

Conclusions
The major contributions of this article can be listed as 1. a methodology to develop loworder models for zones with high-mass radiant floor heating systems; 2. a method for short notice predictive control in response to change in the grid price signal for zones with hydronic floor heating systems; and 3. a method for day-ahead predictive control in response to the change in the grid price signal for the zones with hydronic floor heating systems.
Optimal Control strategies in response to different price signals from the electricity grid with short and long-term predictions were presented to enhance the energy flexibility for a thermal zone equipped with a hydronic floor heating system and significant amount of thermal mass. The flexibility achieved through the floor thermal mass is particularly valuable in interaction with the smart grid to reduce the peak power demand significantly during the critical periods. A methodology to create low order thermal network models was presented and utilized to develop the thermal model of the zone. The model was validated with experimental measurements and subsequently used to study the control strategies. In general, in a building without much knowledge about the details of the construction and zones, a grey-box model similar to the presented second order model can be calibrated in real time by the building automation system to optimize the building performance. Therefore, this modeling methodology can be applied in any location to any zone equipped with a radiant floor heating system.
A method was presented to deal with a short notice price signal that changes every hour during the day and for which the prices can be either higher or lower than the reference operating condition. It was observed that this method can reduce the cost of power consumption while still ensuring the occupant comfort is not violated. Then, an optimization methodology was introduced to minimize the cost function defined based on the predicted price signal for the following days. The optimized floor surface temperature setpoint was found for different price signals to minimize the cost of the heating power (kW) as well as the energy (kWh). It was shown that the strategy can be modified if the price signal prediction would not be as expected and sudden increase in the price happens especially during the critical periods. The full implementation of the developed MPC strategies to a real building will be the next step in the research.  Acknowledgments: This research was supported through NSERC/Hydro-Quebec Industrial Research Chair in "Optimized Operation and Energy Efficiency: Towards High Performance Buildings" at Concordia University.

Conflicts of Interest:
The authors declare no conflict of interest.

A i
Area of surface represented by node i (m 2 ) BEFI Building energy flexibility index C i Thermal capacitance of node i(J/K) k Thermal conductivity of materials (W/(m.k)) Kp Proportional control constant (W/ • C) l i Thickness of surface i (m) P sig Price signal P reference Reference price signal

Q aux
Auxiliary heating source (W) Q i Source entering node i (W) Q max Maximum heating capacity (system size) Q ref Reference heating power profile (W) R i,j Thermal resistance between nodes i and j (m 2 K/W) R inf Infiltration thermal resistance (m 2 W/K) R W Thermal resistance of walls (m 2 W/K) T i Temperature of node i ( • C) T sp Air setpoint temperature ( • C) T o Outdoor temperature ( • C) T floor Floor surface temperature ( • C) t Time