Abstract
This paper presents a methodology to develop archetype gray-box models and use them in an economic model-based predictive control algorithm to simulate optimal heating load management in response to a newly-introduced static time-of-use tariff for Québec’s residential sector, rate Flex-D. The methodology is evaluated through a case study, wherein in situ measurements from a two-storey unoccupied research house of Hydro-Québec are used to develop an 11R6C network with a heuristic zoning-by-floor approach and compute the sequence of optimal electric heating input for the next control horizon. Properly-tuned economic model-based predictive control under rate Flex-D shows potential for an approximately 30% reduction in daily heating cost compared to the reference operation, with a minimal average deviation of indoor air temperature from the reference setpoint. Also, the analysis of the response’s sensitivity to weather forecast uncertainties indicates that the most influential uncontrolled input directing the performance of economic model-based predictive control is the structure price signal, rendering the impact of uncertainty in the weather forecast negligible.
Keywords
Introduction
Energy flexibility enables demand-side management and helps minimize the stress on the grid as demand peaks. Several national and international initiatives, including the International Energy Agency’s Energy in Buildings and Communities (IEA-EBC) Annex 67 (Jensen et al., 2017), Annex 81 (IEA EBC, n.d.-a), Annex 82 (IEA EBC, n.d.-b), Annex 84 (IEA EBC, n.d.-c), and the Grid-Interactive Efficient Buildings (GEB) program by the U.S. Department of Energy (Satchwell et al., 2021), focused on advancing the theoretical and practical aspects of energy flexibility in buildings.
The IEA-EBC Annex 67 has defined building energy flexibility as “the ability to manage short-term demand and generation according to local climate conditions, the user needs, and grid requirements without jeopardizing the technical capabilities of the operating systems and the comfort of occupants.” (Jensen et al., 2017). Li et al. (2023) have comprehensively reviewed data-driven energy flexibility Key Performance Indicators (KPIs) for buildings in operation, categorizing and analyzing them according to type, complexity, scope, key stakeholders, data requirement, baseline requirement, resolution, and popularity.
Advanced control strategies of heating and cooling systems, specifically Model-based Predictive Control (MPC), are proven to reduce energy consumption and, therefore, operational costs while ensuring occupant thermal comfort (Freund and Schmitz, 2021; Seal et al., 2020; Vallianos et al., 2019, 2023). MPC uses a mathematical model of the building and forecasts of weather and other uncontrolled inputs (disturbances) to predict its future behavior and derive optimal control actions (Drgoňa et al., 2020).
MPC has also shown promise in Demand-side Response (DR) applications, as it exploits the building’s inertia, that is, its thermal and electrical energy storage capacity, to shift consumption and thereby increase energy flexibility (Johra and Heiselberg, 2017; Johra et al., 2019; Murray et al., 2021; Oliveira Panão et al., 2019) or manage transactive interaction with the grid (Denholm and Hand, 2011). Challenges posed to the grids due to the electrification of heating and transport, for example, the growing popularity of heat pumps and electric vehicles, are pushing for the increased application of DR measures (Padovani et al., 2021; Ravi and Aziz, 2022).
The main approaches to building energy modeling are (1) physics-driven (white-box), based on the principles of energy/mass conservation and heat transfer, with a comprehensive description of building geometry, systems, and materials; (2) data-driven (black-box), purely based on statistics and machine learning methods, with no/little information about the building physics and systems; and (3) hybrid (gray-box) approach, based on a simplified description of building physics and systems, with physically-meaningful parameters identified from measurements (Li et al., 2021).
Most simulation tools, like EnergyPlus and TRNSYS, use detailed white-box models and therefore require expertise to create, are time-consuming to calibrate (Abtahi et al., 2022a; Candanedo et al., 2022) and require a great deal of knowledge about the building, which is often not available but assumed (Abtahi et al., 2023a, 2023b). For these reasons, they are not ideal candidates for control applications. On the other hand, black-box models are trained with measured data and do not require the expertise of a building engineer (Maddalena et al., 2020). Trained on enough data, black-box models are highly accurate; however, they cannot often extrapolate to situations not “seen” in the training set or from building to building (Afram et al., 2017). This poses a challenge when using black-box models for optimal operations as they are trained on data from suboptimal control strategies or for specific seasons or operating modes (Afram and Janabi-Sharifi, 2014).
Gray-box models require much less data for calibration (Arendt et al., 2018) and are more likely to stay reliable outside the training range (Afroz et al., 2018) compared to black-box models. Compared to the white-box models, gray-box models are faster to set up and calibrate and require much less prior knowledge/information. Therefore, gray-box models have been repeatedly recognized as the best candidates for short-term optimal decision-making and energy flexibility analysis (Abtahi et al., 2022b; Harish and Kumar, 2016; Wang and Chen, 2019).
A gray-box approach is not an oversimplification of the topology but a selection of relevant information. It requires knowledge of the application and its physics and a coherent choice of significant inputs and outputs. Common simplifications in gray-box models to facilitate the representation of buildings’ thermal dynamics are (1) linearization of heat transfer coefficients, (2) temporal discretization, and (3) reduction of spatial dimensions, also known as zoning (Athienitis and O’Brien, 2015). Essentially, an ideal gray-box model is the simplest one that describes all the dominant patterns and information embedded in the measurements (Bacher and Madsen, 2011). The most common practice in gray-box modeling is the use of thermal Resistance-Capacitance (RC) networks (Athienitis et al., 2015).
The heat balance in thermal RC networks is determined through heat exchange between capacitances and various sources that are either controlled or uncontrolled. Defining the network’s structure, that is, the order (number of capacitances) and the heat exchange paths (number and placement of resistances), is not straightforward. Low/high order affects the model’s structural identifiability (Agbi et al., 2012). An oversimplified model may be unable to adequately capture the thermal dynamics in different zones, whereas a high-order model may be over-parameterized and suffer from parameter redundancy, meaning that there are multiple sets of parameter values that yield the same model (Brastein et al., 2018). The structure of thermal RC networks is defined with either (1) a heuristic or (2) an automated approach.
A heuristic approach entails configuring the structure of thermal RC networks to encapsulate merely one or a limited number of zones, where the determination of the order and dominant heat exchange paths is based on either (1) established zoning approaches in the literature, for example, 1R1C or 3R2C (Fonti et al., 2017; Harb et al., 2016; John et al., 2018; Žáčeková et al., 2014) or (2) observations and control expertise for case-novel structures. As an example of the latter, Abtahi et al. (2022a) compared defining the structure based on the number of floors, that is, zoning by floor (6R3C), and based on the orientation of each room, that is, zoning by orientation (7R3C), as two zoning approaches in residential buildings with smart thermostats; in zoning by floor, each floor is a separate thermal zone represented by a lumped thermal capacitance, whereas, in zoning by orientation, the southern and northern rooms are separate thermal zones represented by a lumped thermal capacitance, which particularly helps minimize overheating in southern rooms during clear days with significant expected solar gains.
In contrast, an automated approach attempts to construct the most definitive structure through an iterative process; Bacher and Madsen (2011) defined the most accurate structure with a forward procedure that iteratively compares more complex structures to a base structure and tested it on an experimental setup; Prívara et al. (2012) proposed and validated (with surrogate data from TRNSYS) a two-stage procedure that first selects the most relevant inputs and then identifies the order which yields desirable accuracy; Wang et al. (2019) proposed and validated (with data from a real low-energy house) a procedure that starts with a very complex structure and iteratively deletes parameters based on the asymptotic confidence intervals of their estimates; Leprince et al. (2022) presented a methodology that automatically selects the best model structure from a number of pre-set structures; and, Vallianos et al. (2022) proposed and validated (with data from an experimental house) a methodology for multi-zone model generation that derives the most accurate structure by iteratively adding or removing model parameters.
Candanedo et al. (2022) have defined a Control-oriented Archetype Model (CAM) as a gray-box model that provides a generic representation of common zone/building types and, thus, enough information to evaluate the effect of various control sequences in short-term, helping to make near-optimal decisions. Vallianos et al. (2023) used the low-order CAMs, with a heuristic approach, from Abtahi et al. (2022a) and the multi-zone CAM, with an automated approach, from Vallianos et al. (2022) in an online Economic MPC (EMPC) framework and concluded that although low-order CAMs are faster to create and calibrate, achieving satisfactory accuracy requires more information for data processing and aggregation. On the other hand, the multi-zone CAM is fully automatic and is calibrated directly from data from smart thermostats, but it requires significantly more computational resources, which increase exponentially as the timestep of the model decreases.
Grid operators and control firms are increasingly promoting near-optimal practices to foster effective participation in DR programs, creating a “win-win” scenario for the grid and the customer. A well-designed DR program encourages grid-friendly behaviors without jeopardizing occupants’ comfort and causing technical or security issues; conversely, a poorly-designed program credits customers only at the expense of their comfort, leading to reduced participation. To ensure energy equity and encourage participation, it is essential to weigh the magnitude of potential response against the customers’ marginal benefits in a DR program. Therefore, this paper presents a methodology to develop CAMs and use them in an EMPC algorithm to simulate optimal heating load management in response to a recently introduced DR program for Québec’s residential sector; the methodology is then evaluated through a case study wherein the sensitivity of the response to weather forecast uncertainties is analyzed.
Methodology
Equation (1) is the differential equation governing the heat balance at any node n in a thermal RC network, where C n is the thermal capacitance (J/°C), Rn,m is the thermal resistances (°C/W) in heat exchange with any adjacent node m, T n is the node’s temperature (°C), and Q n is the heat flow (W) to/from the node. Equation (2) represents an explicit finite difference temporal discretization of equation (1), where k and k + 1 denote the current and next time steps, respectively, and δt is the interval between them (s). The heat flow is disaggregated into controlled electric heat input (Qe,n), uncontrolled solar heat gains (Qs,n = αnG), where G is the node’s effective irradiance (W/m2) and α n is the equivalent solar aperture (m2), and other uncontrolled internal gains (Qo,n); this equation serves as an explicit model to predict the node’s temperature at the next time-step, given its current temperature and valid parameter values (θ*). Algorithm 1 explains the procedure for CAM parametrization at every prediction horizon (Npre).
Next, the parametrized (recalibrated) CAM is used in a robust EMPC algorithm to minimize space-heating costs and deviation of indoor air temperature from the reference setpoint (D). The latter is incorporated into the objective function as a soft constraint with a weighting factor (P t ) that balances the cost of electricity against its equivalent cost. Algorithm 2 explains the procedure for EMPC at every control horizon (Ncon).
CAMs may incorporate nodes for capturing the long-term inertia associated with the effective-envelop mass in indoor thermal dynamics, the temperature at which is a hidden state as it is not directly measurable but influences the observed states, that is, the temperature at nodes designed to capture the short-term inertia associated with the effective-air mass in indoor thermal dynamics; a Moving Horizon Estimator (MHE) (Alessandri et al., 2008) leverages the calibrated model, electric heating input records, and weather data from the last κNpre steps to initialize the hidden states within the algorithm.
Case study
Québec’s context
Québec features three distinct climate regions according to the updated Köppen-Geiger climate classification world map (Peel et al., 2007), which are listed in Table 1. Projected to surpass 8.8 million by the end of 2023 (Statistique Québec, 2022), more than 80% of the province’s population reside in an approximate corridor of 300 km long and 100 km wide extending from Greater Montréal (the province’s most populous metropolitan area) to Québec City (the provincial capital) with a humid continental climate (Statistics Canada, 2021). The remaining population lives in scattered communities under subarctic and arctic climates, engaging in specific activities such as forestry, fishing, and mining (Britannica, n.d.).
Québec’s climate regions and characteristics based on the updated Köppen-Geiger climate classification (Peel et al., 2007).
Unsupervised classification of historical weather data yields typical daily patterns across various seasons. K-means clustering algorithms are widely used to partition one dataset into K (a pre-defined number) unlabeled datasets, aiming to maximize intra-similarity within each dataset while minimizing inter-similarity across different datasets, with each dataset represented by its centroid (mean of all the assigned data points).
The Standard Hartigan-Wong K-means clustering algorithm (Hartigan and Wong, 1979) is applied to classify historical weather data of Trois-Rivières from November 2015 to November 2023. Located roughly at the center of Québec’s populous corridor, it serves as an ideal location for identifying daily profiles of ambient air temperature and solar radiation during typical cold and very cold days in southern and south-western Québec with a humid continental climate. Historical weather is retrieved from the SIMEB (in French, SIMulation Energétique des Bâtiments) website (SIMEB, 2023), which provides (among others) ambient air temperature and Global Horizontal Irradiance (GHI) in hourly intervals from several weather stations in Québec; hourly data is interpolated to 15-min intervals using linear interpolation for temperature and backfilling for irradiance.
Figure 1 illustrates the results of historical weather data clustering. The top left plot classifies daily ambient air temperature profiles into three clusters, while the subsequent plots display the annual distribution for each cluster, and the bottom row provides corresponding information for daily GHI profiles. This clustering offers insights into typical daily profiles throughout the year, particularly focusing on the cold period (blue plot). Ambient air temperature and GHI have naturally different annual distributions due to their distinct seasonal and diurnal patterns.

Classification of daily profiles of ambient air temperature and global horizontal irradiance using K-means clustering.
The Québec residential sector contributes 18% of the province’s total energy consumption, with electric space heating expectedly accounting for approximately two-thirds of this share (Whitmore and Pineau, 2023), designating it as the primary end-use to shape the province’s demand profile and put the grid under stress on very cold weekday mornings and evenings (before and after a typical work schedule). In such a context, the high penetration of decentralized electric baseboard heaters in the province’s residential sector can be leveraged to activate a remarkable energy flexibility potential due to their capacity for individual control according to purpose and occupants’ preferences.
To incentivize off-peak consumption, Hydro-Québec, the public utility managing the generation, transmission, and distribution of electricity in Québec, has recently launched a DR program consisting of a static ToU tariff (rate Flex-D) (Hydro-Québec, 2022) and a smart home service with tailored instructions for reducing consumption through their subsidiary, Hilo. This program allows homeowners to actively reduce their cost of electricity consumption by shifting their load to off-peak periods and, as a result, relieve the pressure on the grid. Table 2 explains, in detail, the structure of rate Flex-D.
Structure of rate Flex-D: peak demand events may transpire on weekdays from December 1 to March 31, 6:00–9:00 and/or 16:00–20:00, with a projected range of 25–35 events per year (Hydro-Québec, 2022).
Simulation inputs
Perfect (no uncertainty) historical weather of two weekdays in winter 2023 is used as uncontrolled input to simulate optimal heating load management in response to this DR program. The first day, Wednesday, February 1, was among the 5% coldest days observed between 2015 and 2023 and featured both morning and evening peak demand events, while the second day, Wednesday, March 1, was a typical cold day with only a morning peak demand event. Figure 2 shows the uncontrolled weather input for both days.

Perfect historical weather of two weekdays as uncontrolled input for EMPC simulation: February 1, 2023 (left) and March 1, 2023 (right).
In situ measurements and metadata from an unoccupied research house 35 km north-west of Trois-Rivières are used to develop a CAM for this case study. It is a detached two-storey house with an excavated, uninsulated basement and an attached garage, measuring 7.6 m × 7.9 m (60 m2 footprint) with a floor height of 3 m. The first floor consists of the living room, dining room, kitchen, and a small powder room, while the second floor has three bedrooms and a full bathroom. The house is a representative Québec residence, with R-20 insulation for the walls, R-30 insulation for the roof, 1 and a total of 19 m2 of fenestration consisting of vinyl-framed double-glazing windows with an air gap. It is oriented 35° west of south and is heated with an electric baseboard (η e = 1) in each room controlled by an individual smart thermostat.
Measurements consist of perfect (no uncertainty) smart thermostat records of indoor air temperature, heating setpoint, and electric heating input, with 15-min temporal resolution. A Raspberry Pi 4 serves as an access point (node) within the house Internet-of-Things (IoT) network and directly communicates with the smart thermostats, collecting and storing their data in a local database before transferring to an offline server. The layout of the house is shown in Figure 3, and Table 3 shows the location and installed capacity for each baseboard heater.

Hydro-Québec’s Experimental House for Building Energetics (EHBE) and its floor plans: the basement (top right), the first floor (bottom left), and the second floor (bottom right).
Location and installed capacity of each electric baseboard heater.
The authors have acquired insights from a preceding study (Vallianos et al., 2023), suggesting that the zoning-by-floor approach in Abtahi et al. (2022a) (6R3C network), wherein a single lumped capacitance represents an entire floor in a house, may not accurately capture the building’s thermal dynamics, particularly in free-floating conditions, during which the indoor air discharges considerably faster than the building’s effective thermal mass, depicting it unrealistic to combine the two into one lumped capacitance; therefore, for this case study, an 11R6C network is developed with a heuristic zoning-by-floor approach, where each floor is a separate thermal zone represented by two lumped thermal capacitances: one for capturing the short-term inertia associated with the effective-air mass and another for capturing the long-term inertia associated with the effective-envelop mass. Since the building is unoccupied, the only heat sources considered are the controlled electric heating input and the uncontrolled solar heat gains.
This archetype network is shown in Figure 4, and its information is listed in Table 4.

Case study CAM structure: 11R6C network with the heuristic zoning-by-floor approach.
Case study CAM information: Q e is the electric heat input, and Qs is the solar heat gain.
Algorithm 1 recalibrates the CAM at the beginning of each day (Npre = 24 h), using the records of the past 3 days for training (Ntrain* = 72 h), and initializes with the calibrated values from the previous execution (θ0 = θ*,−1). The training data for each day includes both upward and downward temperature variations within the comfort range of 18–24°C, ensuring a robust representation of indoor thermal dynamics under varying conditions. Calibration results for both days are summarized in Tables 5 and 6.
Case study CAM parameters: valid parameter values for simulation on February 1; R (°C/W), C (MJ/°C), and α (m2).
Case study CAM parameters: valid parameter values for simulation on March 1; R (°C/W), C (MJ/°C), and α (m2).
The capacitance values (MJ/°C) reflect the building’s overall ability to store thermal energy. The values obtained from calibration are consistent with the expected physical properties of building materials and the volume of indoor air, with indoor air nodes corresponding to the air volume in each zone and effective thermal mass nodes representing the heat storage capacity of building materials, showing similar values for the first and second floors and higher values for the basement, consistent with its heavier construction materials. The variations in capacitance values between different days indicate changes in environmental conditions and operational factors, providing an accurate and adaptable model.
The resistance values (°C/W), represent the thermal resistance between different nodes, characterizing the heat transfer within the building. The variability in resistance values between February 1 and March 1 reflects distinct responses to various weather conditions; notably, the resistances toward the ambient air temperature show significant variations between the two cases, with lower values on the colder day (February 1) and higher values on the milder day (March 1).
The solar aperture values (m2), represent the effective area through which solar radiation contributes to the building’s heat gain. These values are distributed between the air and envelope nodes to reflect the building’s design and material properties. Higher solar aperture values for air nodes indicate significant direct solar gains affecting indoor air temperature, while envelope nodes capture the thermal inertia associated with solar heat absorbed by the building materials. The presence of α values on both air and envelope nodes is justified by the need to account for both immediate and delayed effects of solar radiation. The differences in solar aperture values between the 2 days underscore the model’s sensitivity to varying solar conditions, ensuring accurate predictions of the building’s thermal response to solar gains.
Figure 5 compares the actual indoor air temperature (ground truth) in a thick dashed black line against the CAM’s 24 h-ahead predictions in thin colored lines; at each timestep, a thin colored line represents the prediction of indoor air temperature starting from that moment and projecting forward for 24 h based on previous predictions. Nearer-term predictions tend to be more accurate and closely align with the ground truth compared to predictions further into the future, reflecting the accumulation of prediction errors and uncertainties over time, which leads to more significant divergence from the actual observed behavior as the prediction horizon extends. To address this, a decay function is introduced to underscore the significance of nearer-term predictions over those toward the end of the prediction horizon. Therefore, the discrepancy pertains to the variance between further-term predictions and the ground truth. However, it’s worth noting that the effect of this discrepancy is mitigated in the MPC.

Case study CAM performance: indoor air temperature records (ground truth) versus the CAM’s 24 h—ahead predictions—February 1, 2023 (left) and March 1, 2023 (right).
Subsequently, Algorithm 2 uses the recalibrated CAM every 6 h to compute the sequence of optimal electric heating input for the next 6 h (Ncon = 6 h). The daily reference setpoint profile is derived from a recent study (Henao et al., 2018), wherein the analysis of data from 300 smart thermostats situated in 30 houses in Trois-Rivières indicated that the most probable times for heating setpoint adjustments are between 5:00 and 6:00 for morning set-ups, at 8:00 for morning set-backs, between 16:00 and 17:00 for evening set-ups, and at 22:00 for evening set-backs; this study also reported that the most likely heating setpoint is 21°C. Therefore, the daily reference setpoint profile in the first and second floors is established with a high of 21°C, a low of 19°C, a morning set-up at 5:00, a set-back at 8:00, an evening set-up at 17:00, and a set-back at 22:00; in the basement, the setpoint profile is a constant 17°C.
Indoor air temperatures deviating below the reference setpoint or exceeding it by more than 3° (D = 3°C) are penalized with a default equivalent cost of P t = 275 ȼ/°Ch, which implies that 1 °Ch of deviation from the reference setpoint costs approximately as much as 5 kWh of electricity consumption during the peak-demand events and 55 kWh during off-peak periods. All the zones are equally weighted (ωn = μn = 1), while nearer-term predictions are emphasized over predictions toward the end of the prediction horizon with a gentle decay factor (ϒ1 = ϒ2 = 0.01). The reference BaU operation is assumed to reactively maintain the indoor air temperature at the reference setpoint without considering the price of electricity; therefore, it is highly expected that the latter closely tracks the former.
Figure 6 (top plot) depicts the flexible cost-optimal heating demand profile in comparison to the reference heating demand profile, along with their difference (latter−former); this figure (bottom plot) also presents the daily reference setpoint profile versus the predicted indoor air temperature, along with their difference (latter−former). The horizontal green band indicates the range in which deviation from the reference setpoint incurs no penalties, while deviations outside this range are linearly penalized; also, the vertical violet bands highlight peak demand events, where electricity is more expensive than in the white areas. EMPC leverages the lower electricity price to charge the building’s thermal mass during off-peak hours and later relaxes the heating system unless the indoor air temperature approaches the lower limit of the non-penalized range; this significantly reduces heating demand during peak demand events, leading to substantial cost savings.

Case study EMPC performance: reference heating demand profile versus flexible heating demand profile (top); daily reference setpoint profile versus the predicted indoor air temperature (bottom)—February 1, 2023 (left) and March 1, 2023 (right).
Finally, the average demand deviation from t to t + Δt is calculated using equation (5) (Athienitis et al., 2020), where Pe,ref and Pe,flex are the reference and the flexible profiles, respectively; equation (6) calculates the average deviation of indoor air temperature from the reference setpoint for the same period.
The simulation is repeated with two different values for the equivalent price of deviation from the reference setpoint, 1375 and 55 ȼ/°Ch, and the results are summarized in Tables 7 and 8.
Case study EMPC performance: average deviation from the reference BaU operation; ΔPavg (kW) and Davg (°C).
Case study EMPC performance: daily average deviation of indoor air temperature from the reference setpoint versus percentage of daily heating cost reduction.
A negative cost reduction indicates an increase in heating costs compared to the reference scenario.
Table 7 shows a trade-off between P t and ΔPavg; the more deviation of indoor air temperature from the reference setpoint is penalized (higher P t ), the less EMPC uses the building’s thermal mass for pre-heating during off-peak hours and free-floating during peak demand events (lower ΔPavg). As per Table 8, EMPC achieves an approximately 30% reduction in daily heating cost with a minimal average deviation of indoor air temperature from the reference setpoint.
However, introducing a large value for P t may result in an increase in heating cost; using P t = 1375 ȼ/°Ch for the simulation on February 1, the daily heating cost reaches 11.0 $, representing 124% increase compared to the reference daily heating cost of 4.9 $, due to heating during peak demand events, where electricity is 12 times more expensive.
The simulation is repeated, incorporating weather forecast uncertainty for the next control horizon as a noise added to perfect historical weather, with a maximum of 1.5% for the ambient air temperature and a maximum of 15% for the global horizontal irradiance (Voyant et al., 2017); the results indicate that the most influential uncontrolled input directing the performance of EMPC is the structure of rate Flex-D, featuring substantial step changes and rendering the impact of uncertainty in the weather forecast negligible, with less than a 1% difference in predicted cost reduction.
Conclusion
Grid operators and control firms are increasingly promoting near-optimal practices to foster effective participation in DR programs, creating a “win-win” scenario for both the grid and the customer. To ensure energy equity and encourage participation, it is essential to weigh the magnitude of potential response against the customers’ marginal benefits in a DR program. To address this requirement, this paper presented a methodology to develop CAMs and use them in an EMPC algorithm to simulate optimal heating load management in response to a recently-introduced static ToU tariff for Québec’s residential sector, rate Flex-D; the methodology was then evaluated through a case study.
For the case study, perfect in situ measurements from a two-storey unoccupied research house of Hydro-Québec were used to develop an 11R6C network with a heuristic zoning-by-floor approach, where each floor is a separate thermal zone represented by two lumped thermal capacitances: one for capturing the short-term inertia associated with the effective-air mass and another for capturing the long-term inertia associated with the effective-envelop mass.
Given the building’s unoccupancy, heat sources were assumed to be exclusively the controlled electric heating input and the uncontrolled solar heat gains. Perfect historical weather of two weekdays in winter 2023 was used as uncontrolled input for the simulation; the first day, Wednesday, February 1, was among the 5% coldest days observed between 2015 and 2023 and featured both morning and evening peak demand events; and the second day, Wednesday, March 1, was a typical cold day with only a morning peak demand event only. The model was recalibrated at the beginning of each day, using the records of the past 3 days for training and was initialized with the calibrated values from the previous execution.
Subsequently, it was used in EMPC under rate Flex-D every 6 h to compute the sequence of optimal electric heating input for the next 6 h. The daily reference setpoint profile in the first and second floors was established with a high of 21°C, a low of 19°C, a morning set-up at 5:00, a set-back at 8:00, an evening set-up at 17:00, and a set-back at 22:00; in the basement, the setpoint profile was a constant 17°C. Indoor air temperatures deviating below the reference setpoint or exceeding it by more than 3°C were penalized, and the reference operation was assumed to reactively maintain the indoor air temperature at the reference setpoint, without considering the price of electricity.
The daily average deviation of indoor air temperature from the reference setpoint and the percentage of daily heating cost reduction were compared, given different penalization factors (P t ); generally, EMPC achieved an approximately 30% reduction in daily heating cost compared to the reference operation, with a minimal average deviation of indoor air temperature from the reference setpoint. Finally, the analysis of the response’s sensitivity to weather forecast uncertainties indicated that the most influential uncontrolled input directing the performance of EMPC is the structure of rate Flex-D, featuring substantial step changes and rendering the impact of uncertainty in the weather forecast negligible, with less than a 1% difference in predicted cost reduction.
Limitations and future work
The Case Study section primarily serves as a demonstration of using the proposed methodology for recursive day-ahead optimal space-heating load management in the context of the Québec residential sector; the prediction horizon in the case study is 24 h, and at the beginning of every day, the model parameters change, effectively representing a new model iteration (same structure, new parameter values). Therefore, extending the analysis beyond 1 day would not evaluate a unique model but rather multiple iterations of the same structure with varying parameter values. For example, if the analysis were performed over a consecutive 1-week period, it would essentially entail evaluating seven distinct models (one for each day). While this demonstration is limited to 2 days (February 1 and March 1) for the sake of brevity, the methodology is scalable and adaptable to different scenarios. However, incorporating more days into the future analysis offers a more accurate distribution of potential flexibility metrics.
The present study uses fixed spatial and temporal weights in Algorithms 1 and 2, evenly distributing importance among all zones (ωn = μn = 1) and emphasizing near-term predictions with gentle decay factors (ϒ1 = ϒ2 = 0.01). These heuristic values have been effective for this specific case study but may not be optimal for all scenarios, potentially limiting the accuracy and adaptability of the predictive model across different building configurations and seasonal variations. Future work involves refining (tuning) these spatial and temporal weights through sensitivity analysis and optimization techniques to better account for the relative importance of each node within the thermal RC network and to adjust the emphasis on predictions over different horizons. This refinement will consider seasonal variations and differing building configurations to enhance model robustness and performance. Fine-tuning these weights precisely improves the accuracy and applicability of the proposed methodology across a broader range of scenarios.
The current study assumes all measurements to be perfect, without considering any measurement uncertainties. This assumption could affect the robustness and accuracy of the predictive model, as real-world data often contain some degree of error. Future research will address this limitation by incorporating realistic measurement uncertainties, providing a more accurate and resilient model that better reflects practical scenarios.
Another limitation is the lack of comparison between the validation metrics and existing standards. This omission limits the contextual understanding of the model’s performance relative to established benchmarks. To enhance the comprehensiveness of future research, validation metrics will be compared against relevant standards, ensuring a clearer assessment of the model’s accuracy and reliability.
Moreover, the study’s findings are based on simulations, which may not fully represent real-world conditions. As a future step, the methodology presented will be implemented in the research house modeled in this paper. This real-world implementation will validate the simulation results and provide practical insights, ensuring the methodology’s applicability and effectiveness in real-world scenarios.
Footnotes
Appendix
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors would like to gratefully acknowledge Hydro-Québec financial support through an industrial research chair with the Natural Sciences and Engineering Research Council of Canada (NSERC) and a Business Strategy Internship (BSI) program with Mitacs. The authors would also like to thank the research team at Hydro-Québec Research Centre Laboratoire des Technologies de l’Energie (LTE) for their collaborative efforts throughout the project. M.A. would also like to thank the financial support of Fonds de Recherche du Québec-Nature et Technologies (FRQNT) through a Doctoral Research Scholarship.
