A Markov Process Approach to Ensemble Control of Smart Buildings

TL;DR

This paper introduces a Markov Process method for modeling and controlling energy consumption in smart buildings, simplifying complex physics models and enabling real-time data integration for optimal probabilistic management.

Contribution

It presents a novel procedure to convert physical building models into Markov Processes, reducing complexity and parameter dependence while supporting real-time data integration and control.

Findings

Reduced model complexity compared to traditional thermo-physics models

Enabled real-time data integration for system-level analysis

Validated approach using Belgian building data

Abstract

This paper describes a step-by-step procedure that converts a physical model of a building into a Markov Process that characterizes energy consumption of this and other similar buildings. Relative to existing thermo-physics-based building models, the proposed procedure reduces model complexity and depends on fewer parameters, while also maintaining accuracy and feasibility sufficient for system-level analyses. Furthermore, the proposed Markov Process approach makes it possible to leverage real-time data streams available from intelligent data acquisition systems, which are readily available in smart buildings, and merge it with physics-based and statistical models. Construction of the Markov Process naturally leads to a Markov Decision Process formulation, which describes optimal probabilistic control of a collection of similar buildings. The approach is illustrated using validated building data from Belgium.

Equations

$$\displaystyle\underset{\begin{subarray}{c}\rho,\mathcal{P},p,q\end{subarray}}{\text{min}}\sum_{t\in\mathcal{T}}\mathbb{E}_{\rho}\sum_{t\in\mathcal{T}}\!\sum_{\alpha\in\mathcal{A}}\!\left(-U_{t+1}^{\alpha}+\!{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum_{\beta\in\mathcal{A}}}\gamma_{t}^{\alpha\beta}\log\!\frac{\mathcal{P}_{t}^{\alpha\beta}}{\overline{\mathcal{P}}_{t}^{\alpha\beta}}\right)$$

$${\rho }_{t+1}^{\alpha }=\sum _{\beta \in \mathcal{A}}{\mathcal{P}}_t^{\alpha \beta }{\rho }_t^{\beta },\forall \alpha \in \mathcal{A},t\in \mathcal{T}$$

$$p_t=\sum _{\alpha \in \mathcal{A}}{p}^{\alpha }{\rho }_t^{\alpha },\forall t\in \mathcal{T},$$

$$q_t=\sum _{\alpha \in \mathcal{A}}{q}^{\alpha }{\rho }_t^{\alpha },\forall t\in \mathcal{T},$$

$$\sum _{\alpha \in \mathcal{A}}{\mathcal{P}}_t^{\alpha \beta }=1,\forall t\in \mathcal{T},\beta \in \mathcal{A},$$

$$\min \sum _{t\in \mathcal{T}}{d}_{j,t}^{\mathrm{H}}$$

$$\forall t\in \mathcal{T}:d_{j,t}^{\mathrm{H}}=p_{j,t}^{\mathrm{HP}}+p_{j,t}^{\mathrm{A}}$$

$$\forall t\in \mathcal{T}:p_{j,t}^{\mathrm{HP}}=p_{j,t}^{\mathrm{HP},\mathrm{SH}}+p_{j,t}^{\mathrm{HP},\mathrm{HW}}\le \overline{P_j^{\mathrm{HP}}}$$

$$\forall t\in \mathcal{T}:p_{j,t}^{\mathrm{A}}=p_{j,t}^{\mathrm{A},\mathrm{SH}}+p_{j,t}^{\mathrm{A},\mathrm{HW}}\le \overline{P_j^{\mathrm{A}}}$$

$$\forall t\in \mathcal{T}:\sum _{z\in \mathcal{Z}}{\dot{q}}_{j,t,z}^{\mathrm{SH}}=CO{P}_h^{\mathrm{SH}}\cdot p_{hj}^{\mathrm{HP},\mathrm{SH}}+p_{hj}^{\mathrm{A},\mathrm{SH}}$$

$$\forall p\in \mathcal{P},\forall t\in \mathcal{T}:t_{j,p,t}^{\mathrm{SH}}=A_{j,p}\cdot t_{j,p,t-1}^{\mathrm{SH}}$$

$$+\sum _{z\in \mathcal{Z}}{B}_{j,p,z}\cdot {\dot{q}}_{\mathrm{j},\mathrm{t},\mathrm{z}}^{\mathrm{SH}}+E_{j,p,t}^{\mathrm{SH}}$$

$$\forall p\in \mathcal{P},\forall t\in \mathcal{T}:\underline{T_{j,p,t}^{\mathrm{SH}}}\le t_{j,p,t}^{\mathrm{SH}}\le \overline{T_{j,p,t}^{\mathrm{SH}}}$$

$$\forall t\in \mathcal{T}:t_{j,t}^{\mathrm{HW}}=t_{j,t-1}^{\mathrm{HW}}-\frac{TP\cdot G_j}{C_j}\cdot (t_{j,t}^{\mathrm{HW}}-T^{\mathrm{E}})$$

$$+\frac{TP}{C_j}\cdot \left(CO{P}_j^{\mathrm{HW}}\cdot p_{j,t}^{\mathrm{HP},\mathrm{HW}}+p_{j,t}^{\mathrm{A},\mathrm{HW}}-{\dot{Q}}_{j,t}^{\mathrm{D}}\right)$$

$$\forall t\in \mathcal{T}:\underline{T_{j,t}^{\mathrm{HW}}}\le t_{j,t}^{\mathrm{HW}}\le \overline{T_{j,t}^{\mathrm{HP}}}$$

Full text

I Introduction

Controllable electrical and heat appliances in residential and commercial buildings can provide significant flexibility to electrical and heat distribution systems by adjusting their consumption to help meeting operational system limits and alleviate overloads. According to the US Federal Energy Regulatory Commission, the potential peak load reduction from all U.S. retail demand response programs was roughly 31,000 MW/year [1]. In the California Independent System Operator market, the total amount of proxy demand response capacity registered in 2016 was about 160 MW [2]. Currently, the ability of buildings to change their consumption is determined using field trips, remote on/off inspections, and equipment specifications [3]. These practices are labor-intensive and provide coarse, static, often weather-insensitive estimates. As a result, the true flexibility of buildings is often misestimated.

Assessing flexibility that each building can provide requires accounting for electric power and heat dynamics, which are driven by comfort and behavioral preferences of occupants. Currently, there are two large groups of methods to model and forecast building electricity and heat consumption: (i) modeling relevant physical processes (e.g. heat transport, electro-mechanical considerations, Kirchoff’s laws, evaporation, etc) and data-driven (e.g. statistical analyses and inference). Physical models not only use available measurements and static building parameters, e.g. location, floor area, number of stories, detailed information on the heating, ventilation, and air conditioning (HVAC) system, lights, coils, doors and windows, but also operate with specific models that govern dynamics of relevant characteristics [4]. EnergyPlus, for example, is a popular simulation tool for modeling energy needs of buildings using detailed thermo- and mass- modeling of energy flows inside the building [5]. EnergyPlus can also be used for an off-line and off-site analyses to determine set point adjustments of the energy consumption [6]. The advantage of using the physics-based models is in their ability to describe buildings without prior observations. However, the performance of these models is highly sensitive to the number and accuracy of the underlying modeling choices ans assumptions, as well as to input parameters. Physics-based models often require more inputs than existing data acquisition systems can provide [4], and therefore incur significant uncertainties in both model parameters and dynamic processes. Using such models for controlling an ensemble of buildings may lead to computational issues that would prevent their scalability and implementation for real-life decision-making. Due to these shortcomings, it is common to sacrifice modeling accuracy of the physics-based models, which may lead to a loss of their predictive power.

On the other hand, in lieu of the physics-based models, one can use machine learning and statistical modeling to perform data-driven studies of buildings using a vast amount of historical data available at the buildings equipped with smart meters. These models are trained using the historical energy consumption data and other parameters (e.g. weather conditions, daily operational schedules, and control functionality) [3, 6]. Then, the models can be used continuously to learn and predict energy usage from previously observed conditions. Availability of data is crucial for such approaches, especially when attempting to predict consumption with a minimum set of required inputs [3]. Notably, such data is publicly available at an urban scale. For example, New York City’s Local Law 84 (LL84) requires that all commercial (including multi-family) buildings of 50,000 square feet or more must report energy and water consumption on an annual basis. Although this data is very coarse, it has been used in a combination with other building information (e.g. year built, floor area, property-use type, occupancy) to develop more accurate data-driven building models [7]. On the other hand, the data-driven models are data-intensive and building-specific and require large amounts of data for re-training or re-calibration, even when minor changes are made to the buildings. This hinders scalability of the data-driven models and their ability to represent an ensemble of buildings with varying characteristics. Furthermore, numerous studies have revealed that data-driven models may yield discrepancies (up to 100%) between the models outputs and the observed data [7, 8]. To reduce the gap between the prediction and the actual performance, researchers carried out calibration studies to tune the various inputs to match the observations [9, 10, 11]. Nonetheless, calibration is still an over-specified and under-determined problem due to a relatively large number of inputs and a few measurable outputs [3].

Alternatively, machine learning and data-driven techniques can be leveraged to inform physics-based building models. In this manuscript, we construct a Markov Process (MP) to represent the energy consumption of building appliances to assess the building flexibility for various applications in electric power and heat distribution systems. First, the physical building model is used to characterize the MP using the probability transition matrix. Based on this MP, we formulate the Markov Decision Process (MDP) that can in turn be used to optimally control electric and heat appliances within buildings either by the local utility, or third-party aggregators, or building managers. The ultimate objective of this paper is to combine the scalability of the MDP with the accuracy of the physics-based models.

II MDP for Building Modeling

This section describes an approach to reducing multi-parametric and computationally demanding physical building models to a MP model. First, Section II-A reviews the MP approach to model an ensemble of homogeneous electrical or heat appliances and describes how this MP representation can be used for non-disruptive ensemble control. Second, Section II-B describes a physically accurate, state-space model of a building that serves as a benchmark for the proposed approach. Finally, Section II-C details reduction of the physically accurate, state-space model of Section II-B to the MP described in Section II-A.

II-A MP Representation for Ensemble Control

As explored in our previous work [12, 13, 14], an ensemble of homogeneous electrical and heat appliances can be modeled as an MP. The MP modeling is advantageous due to a sufficiently high accuracy and fidelity provided without rendering prohibitive computational complexity. This paper, leveraging the MP representation from [12, 13, 14], assumes that all electrical and heat appliances within the ensemble are geographically co-located, i.e. placed near one another, and the ensemble includes a sufficiently large number of appliances. These two assumptions are hardly restrictive for modeling real-life buildings, especially in urban areas with relatively high population densities, where normally hundreds of appliances are installed. Under these two assumptions, [12, 13, 14] represent the ensemble of electrical or heat appliances as a discrete-time and discrete-space MP that can be extended to an MDP to control the aggregated state of the ensemble.

MP is considered over a finite time horizon with discrete time periods $0\leq t\leq T$ that constitute set $\mathcal{T}$. Accordingly, each ensemble is represented as a MP with a given number of discrete states, as illustrated in Fig. 1. These states are obtained, given operating range and limits of each appliance, by discretizing the range of electric power or heat consumptions within the ensemble. The set of possible states is defined as $\mathcal{A}$ and $\alpha\neq\beta\in\mathcal{A}$ are indices of the individual states. Using the MP representation in Fig. 1, the probability that the ensemble is observed in state $\alpha$ at time $t$ is described as $\rho_{t}^{\alpha}$. The temporal evolution of the ensemble is accounted for by $\mathcal{\overline{P}}_{t,b}^{\alpha\beta}$, representing transition probability from the state $\beta$ at time $t$ to the state $\alpha$ at time $t+1$. Note that $\mathcal{\overline{P}}_{t,b}^{\alpha\beta}$ stands for default evolution of the ensemble, i.e. it is driven by consumption preferences and choices of consumers that use electrical and heat appliances, and does not include the effect of control actions of the utility or third-party aggregators.

Similarly to [12, 13, 14], the MP approach introduced above leads to the following MDP formulation describing control of an aggregated state of the ensemble:

[TABLE]

subject to:

[TABLE]

where $\mathcal{P}_{t}^{\alpha\beta}$, $\rho_{t}^{\alpha}$, $\rho_{t}^{\beta}$, $p_{t}$ and $q_{t}$ are decision variables, $\overline{\mathcal{P}}_{t}^{\alpha\beta}$, $U_{t+1}^{\alpha}$, $p_{t}^{\alpha}$ and $q_{t}^{\beta}$ are given parameters, and operator $\mathbb{E}_{\rho}$ stands for evaluating expectations over $\rho$. Eq. (1)-(5) describes optimal evolution of the ensemble under control from state $\beta$ at time $t$ to state $\alpha$ at time $t+1$. Eq. (1) is the objective function that includes two terms. The first term maximizes the total utility of the ensemble observed in state $\alpha$ at time $t+1$, denoted as $U_{t+1}^{\alpha}$, while the second term introduces cost of the discomfort due to control decisions ${\mathcal{P}}_{t}^{\alpha\beta}$. Decisions ${\mathcal{P}}_{t}^{\alpha\beta}$ represent the state transitions that are requested by the utility or third-party aggregator and, therefore, are likely to deviate from the default evolution described by $\overline{\mathcal{P}}_{t}^{\alpha\beta}$. Naturally, the difference in transitions between the requested and default states, i.e. between ${\mathcal{P}}_{t}^{\alpha\beta}$ and $\overline{\mathcal{P}}_{t}^{\alpha\beta}$, measures the cost of discomfort computed using the Kullback-Leibler (KL) divergence111The choice of the Kullback-Leibler divergence for the discomfort cost is motivated by its wide use for modeling randomness of discrete and continuous time series. Other penalty functions can also be used instead. in the second term of the objective function using penalty $\gamma^{\alpha\beta}$. Respectively, if ${\mathcal{P}}_{t}^{\alpha\beta}=\overline{\mathcal{P}}_{t}^{\alpha\beta}$, there is no discomfort cost and the second term in Eq. (1) is equal to zero. Eq. (2) describes temporal evolution of the ensemble from $t$ to $t+1$, based on the states and transitions illustrated in Fig. 1. Eqs. (3)-(4) compute the expected active and reactive power injections of the TCL ensemble into the distribution system, where parameters $p^{\alpha}$ and $q^{\alpha}$ describe active and reactive power consumptions in the state $\alpha$. Eq. (5) imposes the integrality constraint on the control decisions ${\mathcal{P}}_{t}^{\alpha\beta}$ such that their total probability is equal to one (normalized).

For the sake of simplicity, in the following we set $\gamma^{\alpha\beta}_{t}=1,\forall\alpha\neq\beta\in\mathcal{A},t\in\mathcal{T}$, which converts the MDP in Eq. (1)-(5) to a Linearly Solvable MDP that can be solved analytically (see [15]).

II-B Physical Building Model

Default transition probabilities $\overline{\mathcal{P}}_{t}^{\alpha\beta}$ can be viewed as a reduced-form representation of the actual changes in electrical and heat measurements observed in a real-life building or as an aggregated consumption metric of a detailed, physical building model. Therefore, if one represents a physical model of a building as $\overline{\mathcal{P}}_{t}^{\alpha\beta}$, the MDP described by Eq. (1)-(5) can be used to control the ensemble of appliances within a given building. This section explains the physical building model used below to construct $\overline{\mathcal{P}}_{t}^{\alpha\beta}$ for a given building.

The physical building model used in this paper is a dynamic, linear, state-space model of a building and its heating system constructed based on [16]. This model captures thermal dynamics of a given building envelope and its heating system using a nine-state representation over the two thermal zones (e.g. differentiating between day and night). Accuracy of the model is validated via detailed physical simulations using the IDEAS library [17] in Modelica, as described in [18]. Furthermore, as shown in [16] and [18], this linear state-space model approximates nonlinear dynamics of residential buildings with heat pumps with a sufficient level of accuracy for system-level studies, i.e. from the perspective of electrical or heat distribution systems. For example, these models have been used to access controllability of an ensemble of thermostatically controlled loads in [19].

The building heating system consists of an air-coupled heat pump (HP) and a back-up electric resistance heater (A), which supply heat to the floor heating system in the day and night zones, i.e. space heating (SH), and to the storage tank for domestic hot water (HW) [18, 20]. Under the assumption of perfectly known external and internal states and gains (e.g. outside temperature, solar radiance, temperature inside the building) and user-defined settings (i.e., temperature bounds $\underline{T_{pj}^{SH}}$ and $\overline{T_{pj}^{SH}}$), the consumption of the heat pump ($P^{HP}_{t}$) and auxiliary heater ($P^{A}_{t}$) for space heating and hot water production can be minimized for each building $j\in\mathcal{J}$, where $\mathcal{J}$ is a set of buildings, using the optimization problem formulated in [19] and stated below:

[TABLE]

subject to

[TABLE]

Eq. (7) aggregates the electrical power required by the heat pump ($p^{\mathrm{HP}}_{j,t}$) and auxiliary heater ($p^{\mathrm{A}}_{j,t}$), which are limited to their nameplate capacity $\overline{P^{\mathrm{HP}}_{j}}$ and $\overline{P^{\mathrm{A}}_{j}}$ in (8) and (9). We distinguish between the electrical power required for space heating ($p^{\mathrm{HP,SH}}_{\mathrm{{j,t}}}\geq 0$ and $p^{\mathrm{A,SH}}_{j,t}\geq 0$) and hot water production ($p^{\mathrm{HP,HW}}_{j,t}\geq 0$ and $p^{\mathrm{A,HW}}_{j,t}\geq 0$) supplied by the heat pumps and auxiliary heaters. The electrical power consumed by the heat pump and the auxiliary heater is related to the required thermal power for space heating $\dot{q}^{\mathrm{SH}}_{j,t,z}\geq 0$ in each zone $z$ (set $\mathcal{Z}$) using the coefficient of performance $COP^{\mathrm{SH}}_{j}$ of each building $j$ as given by Eq.  (10). Eq. (11) describes the temperature evolution ($t^{\mathrm{SH}}_{j,p,t}\geq 0$) in each temperature state $p$ (set $P$) based on the thermal power supplied to each building zone ($\dot{q}^{\mathrm{SH}}_{j,t,z}$), thermal losses, internal and external gains $E^{\mathrm{SH}}_{j,p,t}$. Equation (12) enforces the user-defined comfort constraints on the indoor air temperature ($\underline{T^{\mathrm{SH}}_{j,p,t}}$, $\overline{T^{\mathrm{SH}}_{j,p,t}}$). Eq. (13) describes evolution of temperature $t^{\mathrm{HW}}_{j,t}$ in the hot water storage tank, which depends on the thermal conductance of the storage tank, $G_{j}$, the temperature of its surroundings $T^{\mathrm{E}}$, the thermal power supplied to the hot water storage tank, the thermal capacity of the storage tank and its contents $C_{j}$. $\dot{Q}^{\mathrm{D}}_{j,t}$ represents the withdrawal of thermal power from the hot water storage tank due to the hot water consumption, which follows a pre-defined, consumer-specific profile. $TP$ is the temporal resolution of the model. Equation (14) ensures the user-defined availability of hot water at temperatures between $\underline{T^{\mathrm{HW}}_{j,t}}$ and $\overline{T^{\mathrm{HP}}_{j,t}}$.

II-C From building data to MP

Given the physically accurate building model discussed in Section II-B, the goal of this section is to describe reduction of the optimization in Eq. (6)-(14) to MDP in Eq. (1)-(5). The reduction is illustrated in Fig. 2. It includes the following three steps: (i) building data generation and aggregation, (ii) state-space definition and reduction, and (iii) construction and validation of the resulting MP. These steps are further detailed below.

II-C1 Step 1: Building Energy Models

This step applies the physical building model in Eq. (6)-(14) to a given set of buildings to assess their performance under different use scenarios and external and internal gains. To simplify the physical modeling process, we consider an ensemble of physically homogeneous buildings with different occupancy levels and comfort constraints. Each building in the ensemble is influenced by stochastic, external gains, representing outside temperature, solar radiance, wind etc. Similar in their physical properties and exposed to similar climatic geographical conditions, buildings may be described by a slightly different set of parameters and they may also be subject to some external (e.g. weather-driven) and internal (e.g. windows opened/closed, thermostat set-point changes due to user preferences, fluctuations in occupancy etc.) fluctuations. To account for the internal fluctuations, the physical building model was simulated for different occupancy profiles, provided via binary sequence of “presence/absence” states of occupants. These simulations provide a state-space representation, which includes stochasticity in the internal gains.

II-C2 Step 2: State-space definition and reduction

Given a full set of data obtained using the simulations described in Step 1, as well as available building measurements over a given finite time horizon, one constructs a high-dimensional state-space model to control an ensemble of buildings. However, in practice it is not feasible to monitor and capture all state parameters and, therefore, we derive a $d$-dimensional state-space $S\in\mathbbmtt{R}$, where $d$ is a user-defined number which is smaller than the degree of the original building model. Note that the degree of this reduction depends on the original physical model (state variables and their number) and also on the tolerance to modeling inaccuracies. Therefore, in every case, this decision should be made cautiously based on the intended model application. In the following, we will use a two-state model ($d=2$) to simplify visualization of our results. From the physical perspective, the reduced states should constitute a meaningful physical process. Thus, changing a one-state variable should result in changes in another state variable. Therefore, given the optimization described in Eq. (6)-(14), it is reasonable to select the inside temperature ($T_{in}$) and the electric power for heating ($P^{HP}_{j,t}$ or $P^{A}_{j,t}$) because these state variables are mutually dependent (i.e. an increase in electric power should result in a temperature increase). Furthermore, these variables are of importance to consumers since $T_{in}$ defines comfort preferences and $P^{HP}_{j,t}$ and $P^{A}_{j,t}$ affect the electricity bill.

Given these selected state variables, we construct the reduced state-space building model by uniformly discretizing each state variable into numerically ordered ranges to represent their dispatch range, e.g. as shown in Fig. 1. For example, the range for $T_{in}$ is represented by $n$ states and the range for $P^{HP}_{j,t}$ or $P^{A}_{j,t}$ is represented by $m$ states. Then, the list of possible states will be represented by all possible cross-combinations of these states for $T_{in}$, $P^{HP}_{j,t}$ or $P^{A}_{j,t}$, i.e. $n\times m$. Under these assumptions, $\overline{\mathcal{P}}_{t}^{\alpha\beta}$ and ${\mathcal{P}}_{t}^{\alpha\beta}$ will also have dimensions $n\times m$.

II-C3 Step 3: Construct the MP

Finally, we construct the transition probability matrix $\overline{\mathcal{P}}_{t}^{\alpha\beta}$, which describes the steady-state evolution of the system. Given a reduced finite state-space representation obtained in Step 2, the transition probabilities $\overline{P}_{\alpha\beta}$ are computed for all states by tracing the number of transitions between states $\beta$ and $\alpha$ for all available simulations regardless of their duration and time interval. The resulting $d$-dimensional array should then be normalized to meet the integrality criteria of the MP setting.

Note that this approach for constructing the MP may lead to some non-physical states, which cannot be implemented in practice. To avoid such infeasibilities, respective elements of the matrix $\overline{P}_{\alpha\beta}$ are set to 0.

III Case study

This case study aims to illustrate how the proposed MP approach represents an ensemble of buildings and to demonstrate its utility for non-intrusive control of appliances in the building.

III-A Data for Constructing Physical Building Models

We construct a physical building model, described in Section II-B, using parameters for individual buildings, their heating systems, occupancy and behavioral choices of occupants reported in [16, 21, 22, 18, 23]. The key parameters are also summarized below.

We consider an ‘average’ low-energy building, in which the day and night zones have a surface area of $132$ $m^{2}$ and $138$ $m^{2}$, respectively, [24]. All these buildings are assumed to have undergone a renovation of windows, air tightness, walls, floor and roof resulting in an average U-value of $0.3$ $W/K$ and a ventilation rate of $0.4$ $ACH$ (air changes per hour). The heat pump is sized to meet $80\%$ of the peak heat demand and the remainder is assumed to be supplied by the back-up electric resistance heater. As in [25, 26], the COP of the heat pump is calculated assuming a nominal supply water temperature of the floor heating of $35\,^{\circ}C$. The hot water storage tanks have the volume of either $200$ $l$ or $300$ $l$ (vary based on the maximum daily hot water demand at $50\,^{\circ}C$). The maximum supply temperature of the heat pump is $60\,^{\circ}C$.

When occupants are in the building, the lower bounds for the indoor temperature set points are $20\,^{\circ}C$ and $18\,^{\circ}C$ for the day and night zones, respectively, while the upper bounds are $22\,^{\circ}C$ and $20\,^{\circ}C$, respectively [27]. When unoccupied, the lower bound is set to $16\,^{\circ}C$ and the upper bound is unchanged. To account for the stochasticity of internal gains, we generate 52 user behavior profiles using STroBE (https://github.com/open-ideas/StROBe), a statistical tool described in [22], which allows to capture different levels of consumers’ presence (e.g. occupied or not). External gains are modeled using weather data for Uccle, Belgium metered in 2013. The weather measurements include solar gains computed using the model in [28]. The heating load was simulated during four weeks of calendar year 2013 with a 15 minute-time step.

III-B Construction of MP

As explained above, we construct the MP using two state-space variables. In this representation, the power variables, either electric power of the heat pump ($p^{HP}_{j,t}$) or auxiliary heater ($p^{A}_{j,t}$), are considered as independent variables, while the inside building temperature ($T_{in}^{p}$ ) is considered as a dependent variable. This selection of dependent and independent variables allows to use the model of Section II-B that establishes a connection between the temperature and power variables.

Given this modeling assumptions we have experimented with the following cases:

• Case 1: Construct MP for the power consumption of the heat pump ($p^{HP}_{j,t}$), while ignoring the auxiliary heater ($p^{A}_{j,t}$).
• Case 2: Construct MP for the power consumption of the auxiliary heater ($p^{A}_{j,t}$), while ignoring the heat pump ($p^{HP}_{j,t}$).

This split in two cases is motivated by the different operating principles and low-level controls available for the heat pumps and auxiliary heaters that can affect their operations and, thus, thermal dynamics in the buildings. Note that in both cases, the inside building temperature ($T_{in}^{p}$) is a dependent variable.

III-B1 Case 1: MP for $p^{HP}_{j,t}$

Based on our experiments, we established that the most physically accurate parameters of the MP are obtained for the cases when the range of feasible values on the inside temperature is uniformly discretized in either 10 or 17 states, i.e. $\text{card}(\mathcal{P})=10$ or $\text{card}(\mathcal{P})=17$. In this case, the range of possible states for independent variable $p^{HP}_{j,t}$ is also uniformly discretized in 10 or 17 states. Accordingly, the state space can be represented by $10\times 10$ and $17\times 17$ states, respectively, and matrix $\overline{\mathcal{P}}^{\alpha\beta}$ has dimensions of $10\times 10$ and $17\times 17$, respectively.

Using the physical model and occupancy data of buildings described above, we obtain transition probability matrices $\overline{\mathcal{P}}^{\alpha\beta}$ with dimensions $10\times 10$ and $17\times 17$. These transition probability matrices are shown in Fig. 3 and Fig. 4, respectively. Note that as the number of states increases, the relative weight of elements on the main diagonal increases, while the relative weight of the off-diagonal elements reduces. Also note that increasing the number of states leads to a sparser matrix. This sparsity indicates a fairly small number and magnitude of exogenous disturbances under normal building operations. If such disturbances are more frequent and larger due to internal or external gains (e.g. natural disasters), the transition matrix would become less sparse.

III-B2 Case 2: MP for $p^{A}_{j,t}$

We repeat the process of Section III-B1 and obtain the MP for the case with $p^{A}_{j,t}$ modeled as an independent variable. In this case, our experiments suggest that the most physically accurate parameters of the MP are obtained for the cases when the range of feasible states is discretized in 14 and 32 states, respectively. The resulting transition probability matrices are illustrated in Fig. 5 and Fig. 6. Comparing to results in Fig. 3 and Fig. 4, we note that the individual elements of the matrices obtained with different independent variables vary significantly.

III-C MDP for Ensemble Control

The MP constructed in Section III-B can now be used to implement the control of building ensembles as described in Section II-A. Using the solution procedure described in [14], we solve the optimization problem in Eq. (1)-(5) with the values of $\overline{\mathcal{P}}^{\alpha\beta}$ described in Fig. 4 and Fig. 5 and the optimization horizon set to 96 hourly time intervals. Optimal decisions ($\mathcal{P}_{t}^{\alpha\beta}$) for each case lead to the active power consumption displayed in Fig. 7 and Fig. 8, respectively. Similarly to the observations [14], these optimal decisions change over time following changing prices during the optimization horizon. Note that the results in Fig. 7 and Fig. 8 demonstrate that the optimal decisions converge to one or may be a very few states due to more revenues opportunities in these states.

Acknowledgements

The work at New York University was in part supported by NSF Award # CMMI-1825212. The work at LANL was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396. The work was partially supported by DOE/OE/GMLC and LANL/LDRD/CNLS projects.