Demand Flexibility Management for Buildings-to-Grid Integration with Uncertain Generation

Abstract

We present a Buildings-to-Grid (BtG) integration framework with intermittent wind-power generation and demand flexibility management provided by buildings. First, we extend the existing BtG models by introducing uncertain wind-power generation and reformulating the interactions between the Transmission System Operator (TSO), Distribution System Operators (DSO), and buildings. We then develop a unified BtG control framework to deal with forecast errors in the wind power, by considering ancillary services from both reserves and demand-side flexibility. The resulting framework is formulated as a finite-horizon stochastic model predictive control (MPC) problem, which is generally hard to solve due to the unknown distribution of the wind-power generation. To overcome this limitation, we present a tractable robust reformulation, together with probabilistic feasibility guarantees. We demonstrate that the proposed demand flexibility management can substitute the traditional reserve scheduling services in power systems with high levels of uncertain generation. Moreover, we show that this change does not jeopardize the stability of the grid or violate thermal comfort constraints of buildings. We finally provide a large-scale Monte Carlo simulation study to confirm the impact of achievements.

1. Introduction

An essential requirement of power grids with high wind-power penetration is the ancillary reserve power service, which can reduce unwanted power curtailment and enable higher integration of renewable generation. The reserve scheduling task of the Transmission System Operator (TSO) in power grids deals with day-ahead scheduling of the generator reserve power, in order to compensate for mismatches between the forecast and actual wind power [1]. Due to the growing penetration levels of wind-power generation, TSOs need to deal with increasing levels of uncertainty, thus imposing novel challenges and responsibilities for TSOs to avoid blackouts and other contingencies. This trend highlights the necessity for TSOs to introduce new types of ancillary services by enabling end-users (buildings) demand-side flexibility.

Approximately 40% of the global energy is consumed by buildings, with half of this being directly related to Heating, Ventilation, and Air Conditioning (HVAC) [2,3]. Buildings control decisions are typically only optimized locally to minimize energy consumption, but not in the wider scope of optimal control of the electricity grid. Demand-side flexibility of buildings represents the capability to shift production and/or consumption of electricity in time, while still satisfying consumer comfort requirements, and without changing the total production or consumption [4]. Building demand-side flexibility has been identified as an important asset for providing ancillary services being managed by the Distribution System Operators (DSOs) [5,6].

However, TSOs still have a much more prominent role in scheduling ancillary services, and DSOs often miss a majority of the scheduling information [7]. For a resilient, sustainable power system, DSOs need monitoring, simulation and control strategies that allow them to act as an active system operator of the future smart grid. Due to their uncertain and intermittent nature, renewable energy sources, such as wind power, contribute to power imbalance in the grid, making the task of frequency control more difficult [8]. Consequently, DSOs have the inevitable task to engage in providing ancillary services to the grid, in order to secure grid frequency, voltage levels, and the power balance [9].

It is worth noting that flexibility-based services are already available in some countries. For example, in the United Kingdom, National Grid ESO saw significant reductions in the electricity demand as a result of the COVID-19 pandemic [10]. As a response, the system operator introduced a novel Optional Downward Flexibility Management Service, which can curtail energy production from renewables when demand is lower than forecast [11]. Moreover, the European Network Code on Requirements for Generators (RfG) was issued in 2016, to harmonize standards that all power generating units (both conventional and renewable) must respect [12]. To show compliance, all new generators above 50 MW must undergo extensive simulation and testing, before they are allowed to connect to the grid [13]. However, despite these illustrative examples, more fundamental research is required to show the potential and limitations of such flexibility-based services.

The use of consumer appliances to provide dynamical demand response has been studied extensively in recent work, such as [5,6,14]. The potential of using building HVAC control as a source of ancillary services for the grid is shown by [15,16]. In [17,18], the authors develop the so-called Buildings-to-Grid (BtG) integration framework, which couples buildings and grid dynamics explicitly. Similarly, refs. [19,20] study how building demand response programs can be enabled in BtG systems via model predictive control (MPC). However, none of these studies explicitly coupled appliances, e.g., storage systems or building HVAC power load, to the grid dynamics in the presence of uncertain generation, which highlights the potential of developing a unified control framework to handle uncertain renewable energy sources in the BtG setting using building demand-side flexibility.

Research Contributions

In this paper, we propose a unified dynamical BtG framework by modeling the building demand-side flexibility, while also explicitly formulating the hierarchical interactions between the TSO, DSOs, and buildings. Such a demand-side flexibility asset can be considered to be a short-term operating reserve, which may yield an economic benefit in a robust fashion. The main novelty of this paper is that the proposed demand-side flexibility model in the BtG framework with uncertain generation can substitute the traditional reserve scheduling service in power systems with wind farms, without losing stability properties of the power grid and violating the buildings desired thermal comfort. This paper extends our previous work in [21,22], by presenting explicit expressions to quantify the available building demand-side flexibility from specific assets, and by providing more detailed insights in both the theoretical and experimental results.

The remaining part of this paper is organized as follows. In Section 2, we present the BtG model dynamics, and ancillary service deployment via traditional reserves and demand-side flexibility is discussed in Section 3. Using the proposed BtG model dynamics, we provide a stochastic MPC formulation in Section 4, and then a robust tractable reformulation in Section 5 to achieve a desired level of constraints fulfillment with a-priori probabilistic certificates. We simulate an extended version of an IEEE benchmark case study in Section 6, to demonstrate the functionality of our proposed hierarchical integrated energy system model and control technique. Finally, our main conclusions are summarized in Section 7.

2. Buildings-to-Grid System Modeling

Consider a hierarchical electrical energy network consisting of one TSO network together with a wind farm production unit connected to multiple DSO networks, to which individual buildings are connected.

2.1. Wind Integrated TSO Model Dynamics

Consider $\mathcal{T}=\{1,\dots ,n_t\}$ to be the set of TSO buses (nodes), $\mathcal{G}=\{1,\dots ,n_g\}$ to be the set of TSO generators, and $\mathcal{D}=\{1,\dots ,n_d\}$ to be the set of DSOs connected to the TSO network. Define also $\Gamma \in {\mathbb{R}}^{n_t\times n_g}$ to be an incidence matrix such that the entries, ${\Gamma }_{k,m}$, relate generators to the TSO buses, respectively. The set of neighboring nodes of TSO bus k is defined by ${\mathcal{T}}_k^n$.

The TSO network is modeled using the existing swing equation model, e.g., [21,23], where we extend the active power swing equation for bus k by integrating wind-power production in the following form:

$$\begin{array}{cc}\hfill m_k{\ddot{\delta }}_k\left(t\right)+d_k{\dot{\delta }}_k\left(t\right)& ={\sum }_{m=1}^{n_g}{\Gamma }_{k,m}\left(P_{{\mathrm{GR}}_m}\left(t\right)+R_m\left(t\right)\right)+P_{w_k}\left(t\right)-P_{{\mathrm{LD}}_k}\left(t\right)-S_{{\mathrm{LD}}_k}\left(t\right)\hfill \\ \hfill & -{\sum }_{j\in {\mathcal{T}}_k^n}{b}_{kj}sin({\delta }_k\left(t\right)-{\delta }_j\left(t\right)),\hfill \end{array}$$

(1)

where ${\delta }_k\left(t\right)$, ${\dot{\delta }}_k\left(t\right)$ and ${\ddot{\delta }}_k\left(t\right)$ are the voltage angle, angular velocity (frequency), and angular acceleration of bus k, and $m_k$ and $d_k$ are inertia and damping coefficients, respectively. The power flow between buses is considered to be purely reactive, and is characterized by the line susceptance, $b_{kj}=b_{jk}$, and the difference in voltage angle, ${\delta }_k\left(t\right)-{\delta }_j\left(t\right)$. The susceptance is the imaginary part of admittance, which is a measure of how easily a circuit will allow a current to flow through a line [24]. $P_{{\mathrm{GR}}_m}$ denotes the power dispatch of generator m, $P_{w_k}$ is the wind farm power injection at bus k, and $P_{{\mathrm{LD}}_k}$ is the load of any DSO network connected to bus k. Furthermore, $R_m$ and $S_{{\mathrm{LD}}_k}$ are generator reserve and demand-side flexibility power contributions, respectively, which are defined formally in Section 3, and are used to compensate for errors in the wind-power forecast. Without loss of generality, it is assumed that $P_{w_k}$ is a realization of an unknown stochastic process defined on some probability space $(\mathcal{W},\mathfrak{B}(\mathcal{W}),\mathbb{P})$. It is important to note that we do not require the sample space $\mathcal{W}$ and the probability measure $\mathbb{P}$ to be known explicitly, as will be explained later. We only need a finite number of realizations of the uncertain variable $P_{w_k}$, and it is sufficient to consider that they are independent and identically distributed (i.i.d).

The second order differential equation for TSO bus k in Equation (1) can be rewritten as two first-order equations, by defining the angular frequency deviation ${\omega }_k={\omega }_k^{true}-{\omega }_0$, where ${\omega }_k^{true}$ is the absolute frequency of bus k, and ${\omega }_0$ is the synchronous frequency. Using this notation, the complete TSO network dynamics are captured in the following state-space model:

$$\begin{array}{cc}\hfill {\dot{x}}_t\left(t\right)& =A_t{x}_t\left(t\right)+B_{\mathrm{GR}}(P_{\mathrm{GR}}\left(t\right)+R_m\left(t\right))+P_w\left(t\right)-B_{\mathrm{LD}}(P_{\mathrm{LD}}\left(t\right)+S_{{\mathrm{LD}}_k}\left(t\right))+\mathrm{\Psi }\left(\delta \left(t\right)\right),\hfill \end{array}$$

with state variable $x_t\left(t\right)={[{\delta }_1,\dots ,{\delta }_{n_t},{\omega }_1,\dots ,{\omega }_{n_t}]}^{\top }$, and appropriate matrices $A_t$, $B_{\mathrm{GR}}$, and $B_{\mathrm{LD}}$. For the full derivation of the system parameters for a similar model without wind-power integration, the reader is referred to [17,21,25]. The element of vector $\mathrm{\Psi }\left(\delta \left(t\right)\right)\in {\mathbb{R}}^{n_t}$ for bus k is defined as ${\psi }_k={\sum }_{j\in {\mathcal{T}}_k^n}{b}_{kj}sin({\delta }_k\left(t\right)-{\delta }_j\left(t\right))$. $P_{\mathrm{LD}}\left(t\right)\in {\mathbb{R}}^{n_t}$ and $S_{\mathrm{LD}}\left(t\right)\in {\mathbb{R}}^{n_t}$ represent the power demand and flexibility from all DSO networks connected to the TSO, respectively, and are obtained as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{P}_{\mathrm{LD}}\left(t\right)\hfill & =\sum _{i\in \mathcal{D}}^{}{P}_{\mathrm{LD}}^i\left(t\right)=\sum _{i\in \mathcal{D}}^{}\left(P_{\mathrm{IMP}}^i\left(t\right){\mathbf{1}}_{n_d^i}\right)\hfill \\ S_{\mathrm{LD}}\left(t\right)\hfill & =\sum _{i\in \mathcal{D}}^{}{S}_{\mathrm{LD}}^i\left(t\right)=\sum _{i\in \mathcal{D}}^{}\left(S_{\mathrm{IMP}}^i\left(t\right){\mathbf{1}}_{n_d^i}\right)\hfill \end{array}\right.,\end{array}$$

where $P_{\mathrm{LD}}^i$, $S_{\mathrm{LD}}^i\in {\mathbb{R}}^{n_t}$ represent the load and flexibility of DSO i on each TSO bus, respectively. $P_{\mathrm{IMP}}^i\in {\mathbb{R}}^{n_t\times n_d^i}$ with entries $P_{\mathrm{IMP},kl}^i$ denotes the active power flow from TSO bus k to bus l of DSO i, and is nonzero only for border buses (i.e., buses that are physically connected via a line). Similarly, $S_{\mathrm{IMP}}^i$ represents the flexibility power flow between the TSO and DSO i.

Remark 1.

To determine the power flow across the TSO and DSO networks, we follow a DC power approximation and follow the small-angle approximation to assume that $sin({\delta }_k\left(t\right)-{\delta }_j\left(t\right))\approx ({\delta }_k\left(t\right)-{\delta }_j\left(t\right))$. Using such an approximation and letting ${\mathcal{L}}_t$ be the set of lines in the TSO network, we can define the Laplacian matrix $L^{\mathrm{TSO}}\in {\mathbb{R}}^{|{\mathcal{L}}_t|\times n_t}$ for the TSO network with entries defined by

$$L_{k,j}^{\mathrm{TSO}}=\left\{\begin{array}{cc}\sum _{m\in {\mathcal{T}}_k^n}[-b_{km}]\hfill & \mathrm{if}k=j,\hfill \\ b_{kj}\hfill & \mathrm{if}k\ne j,j\in {\mathcal{T}}_k^n,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(2)

We can now construct the following augmented matrix which can be used to obtain the power flow in the TSO network lines:

$${\tilde{L}}^{\mathrm{TSO}}=\left[\begin{array}{cc}{L}^{\mathrm{TSO}}& {\mathbf{0}}_{n_t\times n_t}\\ {\mathbf{0}}_{n_t\times n_t}& {\mathbf{0}}_{n_t\times n_t}\end{array}\right].$$

(3)

A similar expression can be written for DSO network to obtain ${\tilde{L}}^{\mathrm{DSO},i}$ for DSO i.

2.2. DSO Model Dynamics

Consider ${\mathcal{D}}_i=\{1,\dots ,n_d^i\}$ to be the set of buses (nodes) of DSO i. Denote by ${\mathcal{D}}_l^{n,i}$ the neighborhood set of node l of DSO i. Each DSO network is connected to the parent TSO through one or multiple border buses. DSO network dynamics are the same as in Equation (1), but with the additional assumption that no generators are connected to DSO buses, leading to $m=0$ and $d=0$ in the DSO swing equations. Please note that this is not a restrictive assumption, as the model can easily be extended to include generation in the DSO network as well. The dynamics at bus l of DSO i are described as

$$\begin{array}{cc}\hfill {\widehat{d}}_l^i{\dot{\delta }}_l^i\left(t\right)& ={\left[{[P_{\mathrm{IMP}}^i\left(t\right)+S_{\mathrm{IMP}}^i\left(t\right)]}^{\top }{\mathbf{1}}_{n_t}\right]}_l-P_{{\mathrm{BD}}_l}^i\left(t\right)-S_{{\mathrm{BD}}_l}^i\left(t\right)-{\sum }_{j\in {\mathcal{D}}_l^{n,i}}{b}_{lj}^{i}sin({\delta }_l^i\left(t\right)-{\delta }_j^i\left(t\right)),\hfill \end{array}$$

(4)

where the subscript in ${\left[·\right]}_l$ denotes element l of the respective vector, ${\widehat{d}}_l^i{\dot{\delta }}_l^i\left(t\right)$ is the frequency-sensitive portion of the uncontrollable load at bus l, and $P_{{\mathrm{BD}}_l}^i\left(t\right)$, $S_{{\mathrm{BD}}_l}^i\left(t\right)$ denote the total building demand and flexibility at bus l of DSO i. The resulting state-space model for DSO i is given by:

$$\begin{array}{cc}\hfill {\dot{x}}_d^i\left(t\right)& =A_d^i{x}_d^i\left(t\right)+B_d^i\left({\left[P_{\mathrm{IMP}}^i\left(t\right)+S_{\mathrm{IMP}}^i\left(t\right)\right]}^{\top }{\mathbf{1}}_{n_t}-P_{\mathrm{BD}}^i\left(t\right)-S_{\mathrm{BD}}^i\left(t\right)-{\mathrm{\Psi }}^i\left({\delta }^i\left(t\right)\right)\right),\hfill \end{array}$$

where $x_d^i\left(t\right)={[{\delta }_1^i,\dots ,{\delta }_{n_d^i}^i,{\omega }_1^i,\dots ,{\omega }_{n_d^i}^i]}^{\top }$ is the state vector, and $A_d^i$ and $B_d^i$ are system matrices. The parameter ${\mathrm{\Psi }}^i\left({\delta }^i\left(t\right)\right)\in {\mathbb{R}}^{n_d^i}$ and its element for bus l are defined as ${\psi }_l^i={\sum }_{j\in {\mathcal{D}}_l^{n,i}}{b}_{lj}^{i}sin({\delta }_l^i\left(t\right)-{\delta }_j^i\left(t\right))$. The total building demand $P_{\mathrm{BD}}^i\left(t\right)$ and flexibility $S_{\mathrm{BD}}^i\left(t\right)$ at node l can be determined as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{P}_{\mathrm{BD}}^i\left(t\right)\hfill & ={[P_{{\mathrm{BD}}_1}^i,\dots ,P_{{\mathrm{BD}}_{n_d^i}}^i]}^{\top }\hfill \\ \hfill & ={\mathrm{\Pi }}^i(P_{\mathrm{stor}}\left(t\right)+P_{\mathrm{hvac}}\left(t\right)+P_{\mathrm{misc}})\hfill \\ S_{\mathrm{BD}}^i\left(t\right)\hfill & ={[S_{{\mathrm{BD}}_1}^i,\dots ,S_{{\mathrm{BD}}_{n_d^i}}^i]}^{\top }\hfill \\ \hfill & ={\mathrm{\Pi }}^i(S_{\mathrm{stor}}\left(t\right)+S_{\mathrm{hvac}}\left(t\right))\hfill \end{array}\right.,\end{array}$$

where incidence matrix ${\mathrm{\Pi }}^i\in {\mathbb{R}}^{n_d^i\times n_b}$ relates the buses of DSO i to the $n_b$ buildings in the network. Each building is assumed to be connected to exactly one DSO bus, and thus, ${\sum }_{i\in \mathcal{D}}\left[{\left({\mathrm{\Pi }}^i\right)}^{\top }{\mathbf{1}}_{n_d^i}\right]={\mathbf{1}}_{n_b}$. $P_{\mathrm{stor}}\left(t\right)$, $P_{\mathrm{hvac}}\left(t\right)$, $P_{\mathrm{misc}}\left(t\right)$$\in {\mathbb{R}}^{n_b}$ are the power storage, HVAC power demand, and uncontrollable miscellaneous power consumption of all buildings, respectively. Finally, $S_{\mathrm{stor}}\left(t\right)$, $S_{\mathrm{hvac}}\left(t\right)\in {\mathbb{R}}^{n_b}$ represent the total storage and the total building HVAC flexibility, respectively.

2.3. Building and Storage Model Dynamics

2.3.1. Building Thermal Comfort Model

The framework developed in this paper explicitly couples building decision variables to the power grid. Traditional resistance and capacitor (RC) networks are widely used to model thermodynamics of building envelopes [15]. Both highly detailed models with variables for many zones, e.g., [26,27,28], and low order models with a single lumped zone, e.g., [17,18], have been proposed in the literature. As the envisioned purpose of the framework is to integrate large clusters of buildings into the grid, a low order thermal building model is most appropriate, considering the resulting computational complexity. Using the 3R-2C circuits model for each building adopted from [17], consider the thermal dynamics:

$$\begin{array}{cc}\hfill {\dot{T}}_{\mathrm{wall}}\left(t\right)& =\frac{T_{\mathrm{amb}}\left(t\right)-T_{\mathrm{wall}}\left(t\right)}{C_{\mathrm{wall}}{R}_2}+\frac{T_{\mathrm{zone}}\left(t\right)-T_{\mathrm{wall}}\left(t\right)}{C_{\mathrm{wall}}{R}_1}+\frac{{\dot{Q}}_{\mathrm{sol}}\left(t\right)}{C_{\mathrm{wall}}},\hfill \\ \hfill {\dot{T}}_{\mathrm{zone}}\left(t\right)& =\frac{T_{\mathrm{wall}}\left(t\right)-T_{\mathrm{zone}}\left(t\right)}{C_{\mathrm{zone}}{R}_1}+\frac{T_{\mathrm{amb}}\left(t\right)-T_{\mathrm{zone}}\left(t\right)}{C_{\mathrm{zone}}{R}_{\mathrm{win}}}+\frac{{\dot{Q}}_{\mathrm{int}}\left(t\right)-{\dot{Q}}_{\mathrm{hvac}}\left(t\right)}{C_{\mathrm{zone}}},\hfill \end{array}$$

(5)

where $T_{\mathrm{wall}}\left(t\right)$ and $T_{\mathrm{zone}}\left(t\right)$ are the wall and zone (interior room) temperatures, respectively, and $T_{\mathrm{amb}}\left(t\right)$ is the ambient temperature. The resistance of the external walls, internal walls, and windows are given by the resistance parameters $R_1$, $R_2$ and $R_{\mathrm{win}}$. The lumped thermal capacities of the building exterior walls and the zone are denoted by $C_{\mathrm{wall}}$ and $C_{\mathrm{zone}}$, respectively. Moreover, ${\dot{Q}}_{\mathrm{sol}}\left(t\right)$ represents the sum of the solar radiation absorbed on the external walls, and ${\dot{Q}}_{\mathrm{int}}\left(t\right)$ is the total heat gain from internal sources. The room temperature is controlled by the cooling load ${\dot{Q}}_{\mathrm{hvac}}\left(t\right)$, which is proportional to the HVAC power consumption via ${\dot{Q}}_{\mathrm{hvac}}\left(t\right)={\mu }_{\mathrm{hvac}}(P_{\mathrm{hvac}}\left(t\right)+S_{\mathrm{hvac}}\left(t\right))$. The thermal dynamics in Equation (5) for an individual building are written as the following state-space model:

$${\dot{x}}_b^l\left(t\right)=A_b^l{x}_b^l\left(t\right)+B^l(P_{\mathrm{hvac}}^l\left(t\right)+S_{\mathrm{hvac}}^l\left(t\right))+{\overline{B}}^l{\overline{w}}_b^l\left(t\right),$$

(6)

where $x_b^l\left(t\right)={\left[T_{\mathrm{wall}}^l{T}_{\mathrm{zone}}^l\right]}^{\top }$ is the state vector and ${\overline{w}}_b^l={\left[T_{\mathrm{amb}}^l{\dot{Q}}_{\mathrm{sol}}^l{\dot{Q}}_{\mathrm{int}}^l\right]}^{\top }$ is an uncertain vector. Please note that for the buildings model, we consider ${\overline{w}}_b^l$ to be equivalent to the forecast value of the uncertain parameters. For the explicit derivation of system matrices $A_b^l$, $B^l$, and ${\overline{B}}^l$, we refer the reader to [21]. Since we aim to describe the dynamics of clusters of buildings, let $n_b$ be the total number of buildings, and denote by $\mathcal{B}=\{1,\dots ,n_b\}$ the full set of buildings connected to the grid. In the absence of communication between buildings, the dynamics of all buildings combined are described by the following state-space model with block diagonal system matrices (full matrix definitions are omitted for brevity):

$${\dot{x}}_b\left(t\right)=A_b{x}_b\left(t\right)+B(P_{\mathrm{hvac}}\left(t\right)+S_{\mathrm{hvac}}\left(t\right))+\overline{B}{\overline{w}}_b\left(t\right).$$

(7)

2.3.2. Electrical Energy Storage Model

Define $x_s^j\left(t\right)\in \mathbb{R}$ to be the energy state variable (state of charge, or SoC) and $P_{\mathrm{stor}}^j\left(t\right),S_{\mathrm{stor}}^j\left(t\right)\in \mathbb{R}$ the normal and flexibility power input rate variable, respectively, for a dedicated electrical storage unit j of building $j\in \mathcal{B}$. Consider now the discrete-time dynamical model of storage unit j as:

$$\begin{array}{c}\hfill x_s^j(k+1)=\xi x_s^j\left(k\right)+h\eta (P_{\mathrm{stor}}^j\left(k\right)+S_{\mathrm{stor}}^j\left(k\right)),\end{array}$$

(8)

where $x_s^j\left(0\right)=x_s^{j,0}$ is the given initial SoC, h is the discretization step size, and $\xi$ and $\eta$ are efficiency coefficients of the storage unit. Although the current equation models energy loss proportionally, higher order equations for losses can be modeled as well. Denote by $\mathrm{\Xi }$ and $\mathrm{\Omega }$ the diagonal matrices composed of ${\xi }^l$ and ${\eta }^l$ for all buildings $l\in \mathcal{B}$, respectively.

Remark 2.

The proposed models for the TSO, DSO, and building thermal comfort dynamics in Equations (1), (4), and (6), respectively, are continuous-time. To formulate a discrete optimal control problem in the MPC paradigm, the dynamics are discretized using the first-order backward Euler implicit method proposed in [29]. The discretized dynamics in the three distinct model areas are then denoted by the following functions:

$$\begin{array}{c}{x}_t(k+1)=f_t\left(x_t\left(k\right),P_{\mathrm{GR}}\left(k\right),R\left(k\right),P_{\mathrm{LD}}\left(k\right),S_{\mathrm{LD}}\left(k\right),P_w\left(k\right)\right),\hfill \\ x_d^i(k+1)=f_d^i\left(x_d^i\left(k\right),P_{\mathrm{IMP}}^i\left(k\right),S_{\mathrm{IMP}}^i\left(k\right),P_{\mathrm{BD}}^i\left(k\right),S_{\mathrm{BD}}^i\left(k\right)\right),\hfill \\ x_b(k+1)=f_b\left(x_b\left(k\right),P_{\mathrm{hvac}}\left(k\right),S_{\mathrm{hvac}}\left(k\right)\right),\hfill \end{array}$$

(9)

for all time steps $k=1,2,\dots ,N_{\mathrm{sim}},$ such that $N_{\mathrm{sim}}$ is the finite-time step of the simulation study.

3. Ancillary Service Deployment

In this section, we first describe common practice in the TSO network to deal with highly fluctuating wind-power integration, the so-called reserve scheduling service, and then present a novel formulation for building demand flexibility to be deployed as an ancillary service to handle uncertain wind power in the energy network.

3.1. Reserve Scheduling Formulation

Wind-power generation suffers from uncertainty and limited predictability [30]. Consider $P_w^f\in {\mathbb{R}}^{n_w}$ to be the forecast value of the wind power for every wind farm in the set $\mathcal{F}=\{1,\dots ,n_w\}$, and denote the error between forecast and actual power by $\mathrm{\Delta }{P}_w=P_w-P_w^f$. When $\mathrm{\Delta }{P}_w\ne 0$, the power balance in the TSO network between production and demand is not satisfied anymore. To restore power balance, the common practice is to deploy so-called reserve power such that by adjusting the generator power output, forecast errors are compensated with altered generation (denoted by $R\in {\mathbb{R}}^{n_g}$). Increasing the generator output is called up-spinning reserve ($R^{\mathrm{us}}\ge 0$), and decreasing the output is down-spinning reserve ($R^{\mathrm{ds}}\ge 0$). To schedule the reserve power, one can define the reserve power as follows:

$$-\sum _{i\in \mathcal{G}}{R}_i\left(k\right)=\sum _{m\in \mathcal{F}}\mathrm{\Delta }{P}_{w_m}\left(k\right),$$

(10)

and consider $-R^{\mathrm{ds}}\left(k\right)\le R\left(k\right)\le R^{\mathrm{us}}\left(k\right),$ to determine $R^{\mathrm{ds}}\left(k\right)$ and $R^{\mathrm{us}}\left(k\right)$ at each time step k. Please note that the new scheduled power should also satisfy the generation limits, i.e., $P_{\mathrm{GR}}^{min}\le P_{\mathrm{GR}}\left(k\right)+R\left(k\right)\le P_{\mathrm{GR}}^{max}$, where $P_{\mathrm{GR}}^{min}$ and $P_{\mathrm{GR}}^{max}$ are the minimum and maximum limits of power generation units, respectively.

3.2. Building Flexibility Formulation

Define $S\in {\mathbb{R}}^{n_b}$ to be the building flexibility and consider the capability to increase energy consumption as increased-demand flexibility ($S^{\mathrm{id}}\ge 0)$ and the capability to decrease energy consumption as decreased-demand flexibility ($S^{\mathrm{dd}}\ge 0$). If instead of generator reserve power, demand flexibility is used to mitigate the wind-power forecast error, one can schedule the flexibility as follows:

$$\sum _{l\in \mathcal{B}}{S}_l\left(k\right)=\sum _{m\in \mathcal{F}}\mathrm{\Delta }{P}_{w_m}\left(k\right),$$

(11)

and consider $-S^{\mathrm{dd}}\left(k\right)\le S\left(k\right)\le S^{\mathrm{id}}\left(k\right),$ to determine $S^{\mathrm{dd}}\left(k\right)$ and $S^{\mathrm{id}}\left(k\right)$ at each time step k.

Consider now the flexibility contribution of individual buildings in the network via two sources: (1) building storage systems, and (2) building HVAC loads, which yields the following bounds on the building flexibility at time step k:

$$-S_{\mathrm{stor}}^{\mathrm{dd}}\left(k\right)-S_{\mathrm{hvac}}^{\mathrm{dd}}\left(k\right)\le S\left(k\right)\le S_{\mathrm{stor}}^{\mathrm{id}}\left(k\right)+S_{\mathrm{hvac}}^{\mathrm{id}}\left(k\right),$$

(12)

where $S_{\mathrm{stor}}^{\mathrm{id}}$, $S_{\mathrm{hvac}}^{\mathrm{id}}$, $S_{\mathrm{stor}}^{\mathrm{dd}}$, and $S_{\mathrm{hvac}}^{\mathrm{dd}}$ represent the increased- and decreased-demand flexibility using the storage unit and HVAC load, respectively. The available storage flexibility of building $l\in \mathcal{B}$ can be formulated as:

$$\begin{array}{cc}\hfill S_{\mathrm{stor}}^{\mathrm{id},l}\left(k\right)& =min\left\{P_{\mathrm{stor}}^{\max }-P_{\mathrm{stor}}^l\left(k\right),\frac{x_s^{\max }-\xi x_s^l\left(k\right)}{h\eta }-P_{\mathrm{stor}}^l\left(k\right)\right\},\hfill \\ \hfill S_{\mathrm{stor}}^{\mathrm{dd},l}\left(k\right)& =max\left\{P_{\mathrm{stor}}^{\min }-P_{\mathrm{stor}}^l\left(k\right),\frac{x_s^{\min }-\xi x_s^l\left(k\right)}{h\eta }-P_{\mathrm{stor}}^l\left(k\right)\right\},\hfill \end{array}$$

(13)

where h is the discretization step size, and $x_s^{\min }$, $x_s^{\max }$, $P_{\mathrm{stor}}^{\min }$, and $P_{\mathrm{stor}}^{\max }$ are the SoC limits and storage power limits, respectively. Storage flexibility is defined as the margin between the current storage power $P_{\mathrm{stor}}\left(k\right)$ and the upper and lower bounds, respectively, constrained by the SoC limits on $x_s\left(k\right)$. As indicated in Figure 1, increased-demand storage flexibility ($S_{\mathrm{stor}}^{\mathrm{id}}$) is the maximum amount of power with which $P_{\mathrm{stor}}\left(k\right)$ can be increased for one time step width h, such that both the maximum storage power injection constraint (Figure 1a) and the storage energy limit (Figure 1b) are still satisfied.

A similar expression is derived for the available HVAC flexibility of building $l\in \mathcal{B}$, such that both HVAC power limits and building comfort levels are still satisfied:

$$\begin{array}{cc}\hfill S_{\mathrm{hvac}}^{\mathrm{id},l}\left(k\right)& =min\left\{P_{\mathrm{hvac}}^{\max }-P_{\mathrm{hvac}}^l\left(k\right),P_{\mathrm{hvac}}^{\mathrm{id},l}\left(k\right)-P_{\mathrm{hvac}}^l\left(k\right)\right\},\hfill \\ \hfill S_{\mathrm{hvac}}^{\mathrm{dd},l}\left(k\right)& =max\left\{P_{\mathrm{hvac}}^{\min }-P_{\mathrm{hvac}}^l\left(k\right),P_{\mathrm{hvac}}^{\mathrm{dd},l}\left(k\right)-P_{\mathrm{hvac}}^l\left(k\right)\right\},\hfill \end{array}$$

(14)

where $P_{\mathrm{hvac}}^{min}$, $P_{\mathrm{hvac}}^{max}$ are the minimum and maximum limits of HVAC power usage, respectively. $P_{\mathrm{hvac}}^{\mathrm{dd},l}\left(k\right)$ and $P_{\mathrm{hvac}}^{\mathrm{id},l}\left(k\right)$ are the minimum and maximum HVAC power such that the comfort level constraints are still satisfied:

$$\begin{array}{cc}\hfill P_{\mathrm{hvac}}^{\mathrm{dd}}\left(k\right)& =\frac{T_{\mathrm{wall}}-T_{\mathrm{zone}}}{R_1{\mu }_{\mathrm{hvac}}}+\frac{T_{\mathrm{amb}}-T_{\mathrm{zone}}}{R_{\mathrm{win}}{\mu }_{\mathrm{hvac}}}+\frac{{\dot{Q}}_{\mathrm{int}}}{{\mu }_{\mathrm{hvac}}}-\frac{C_{\mathrm{zone}}}{h{\mu }_{\mathrm{hvac}}}\left[\min \left\{T_{\mathrm{zone}}^{\max }\left(k\right),T_{\mathrm{zone}}^{\max }(k+h)\right\}-T_{\mathrm{zone}}\left(k\right)\right],\hfill \\ \hfill P_{\mathrm{hvac}}^{\mathrm{id}}\left(k\right)& =\frac{T_{\mathrm{wall}}-T_{\mathrm{zone}}}{R_1{\mu }_{\mathrm{hvac}}}+\frac{T_{\mathrm{amb}}-T_{\mathrm{zone}}}{R_{\mathrm{win}}{\mu }_{\mathrm{hvac}}}+\frac{{\dot{Q}}_{\mathrm{int}}}{{\mu }_{\mathrm{hvac}}}-\frac{C_{\mathrm{zone}}}{h{\mu }_{\mathrm{hvac}}}\left[\max \left\{T_{\mathrm{zone}}^{\min }\left(k\right),T_{\mathrm{zone}}^{\min }(k+h)\right\}-T_{\mathrm{zone}}\left(k\right)\right],\hfill \end{array}$$

where $T_{\mathrm{zone}}^{\max }\left(k\right)$ and $T_{\mathrm{zone}}^{\min }\left(k\right)$ represent the thermal comfort level limits at time k. Positive HVAC flexibility is visualized in Figure 2, showing that the limiting term is either the HVAC power (Figure 2a) or the thermal comfort level (Figure 2b).

3.3. Reserve Scheduling Together with Building Flexibility

Including both reserve and flexibility in the grid power balance yields

$$\sum _{i\in \mathcal{G}}(P_{{\mathrm{GR}}_i}+R_i)+\sum _{m\in \mathcal{F}}(P_{w_m}^f+\mathrm{\Delta }{P}_{w_m})=\sum _{l\in \mathcal{B}}(P_{{\mathrm{BD}}_l}+S_l),$$

(15)

where the sum of reserve and flexibility mitigates the total wind-power forecast error:

$$\sum _{l\in \mathcal{B}}{S}_l-\sum _{i\in \mathcal{G}}{R}_i=\sum _{m\in \mathcal{F}}\mathrm{\Delta }{P}_{w_m}.$$

(16)

Based on Equation (12), the following constraint encodes that the scheduled reserve and flexibility is always sufficient to compensate the wind-power error:

$$-\sum _{i\in \mathcal{G}}{R}_i^{\mathrm{us}}-\sum _{l\in \mathcal{B}}{S}_l^{\mathrm{dd}}\le \sum _{m\in \mathcal{F}}\mathrm{\Delta }{P}_{w_m}\le \sum _{i\in \mathcal{G}}{R}_i^{\mathrm{ds}}+\sum _{l\in \mathcal{B}}{S}_l^{\mathrm{id}}.$$

(17)

4. Stochastic MPC Formulation

Model predictive control (MPC) is a flexible paradigm that defines receding-horizon-based optimization problems, enabling the specification of time-domain objectives together with the ability to explicitly enforce constraints on system dynamics. Two extensions to MPC exist when the system dynamics and/or constraints are subject to uncertainties, namely robust and stochastic MPC. Robust MPC is able to handle the uncertainties using the so-called worst-case approach while still ensuring that the state constraints are met. Alternatively, stochastic MPC has attractive features, due to its ability to handle uncertain systems in a less conservative way. Stochastic MPC considers the stochastic characteristics of the uncertainties and treats system constraints in a probabilistic manner, i.e., using chance constraints. In this section, we formulate a receding horizon stochastic optimization problem to compute an optimal decisions sequence that minimizes a given objective function, subject to the uncertain BtG model dynamics and chance constraints.

Given the set of prediction time steps ${\mathcal{N}}_h=\{0,1,...,N_h\}$ such that $N_h$ is the length of the prediction horizon, consider the concatenated vector of control decision variables to be

$$\begin{array}{c}\hfill {\mathsf{U}}_k={[P_{\mathrm{GR}}(\ell |k),P_{\mathrm{hvac}}(\ell |k),P_{\mathrm{stor}}(\ell |k),R^{\mathrm{us}}(\ell |k),R^{\mathrm{ds}}(\ell |k),S_{\mathrm{hvac}}^{\mathrm{id}}(\ell |k),S_{\mathrm{hvac}}^{\mathrm{dd}}(\ell |k),S_{\mathrm{stor}}^{\mathrm{id}}(\ell |k),S_{\mathrm{stor}}^{\mathrm{dd}}(\ell |k)]}_{\ell \in {\mathcal{N}}_h}.\end{array}$$

The objective function consists of two parts. The first part penalizes TSO and DSO grid frequency deviations, while part two accounts for the operating costs, e.g., power generation, etc. The initial grid state variables are given by $x_k:={[x_t\left(0\right|k),{\left[x_d^i\left(0\right|k)\right]}_{\forall i\in \mathcal{D}}]}^{\top }$, the vector of uncertainty is ${\mathsf{P}}_{w_k}:=\left[P_w{(\ell |k)}_{\ell \in {\mathcal{N}}_h}\right]$, and $Q_t$, ${\left\{Q_d^i\right\}}_{i\in \mathcal{D}}$, $Q_{\mathrm{GR}}$, and $Q_{\mathrm{hvac}}$ are diagonal cost matrices associated with the TSO and DSO states, generation, and HVAC power usage respectively. Finally, $q_{\mathrm{S}}^{\mathrm{id}}$, $q_{\mathrm{S}}^{\mathrm{dd}}$, $q_{\mathrm{R}}^{\mathrm{us}}$, and $q_{\mathrm{R}}^{\mathrm{ds}}$ are the cost vectors related to increased- and decreased-demand, up- and down-spinning reserves, respectively. Using these definitions, the cost function can be formulated as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \tilde{J}(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k})=\sum _{\ell \in {\mathcal{N}}_h}(& x_t^{\top }(\ell |k)Q_t{x}_t(\ell |k)+\sum _{i=1}^{n_d}\left(x_d^{i\top }(\ell |k)Q_d^i{x}_d^i(\ell |k)\right)+P_{\mathrm{GR}}^{\top }(\ell |k)Q_{\mathrm{GR}}{P}_{\mathrm{GR}}(\ell |k)\hfill \\ \hfill & +P_{\mathrm{hvac}}^{\top }(\ell |k)Q_{\mathrm{hvac}}{P}_{\mathrm{hvac}}(\ell |k)+q_{\mathrm{S}}^{\mathrm{id}\top }\left(S_{\mathrm{stor}}^{\mathrm{id}}(\ell |k)+S_{\mathrm{hvac}}^{\mathrm{id}}(\ell |k)\right)\hfill \\ \hfill & +q_{\mathrm{S}}^{\mathrm{dd}\top }\left(S_{\mathrm{stor}}^{\mathrm{dd}}(\ell |k)+S_{\mathrm{hvac}}^{\mathrm{dd}}(\ell |k)\right)+q_{\mathrm{R}}^{\mathrm{us}\top }{R}^{\mathrm{us}}(\ell |k)+q_{\mathrm{R}}^{\mathrm{ds}\top }{R}^{\mathrm{ds}}(\ell |k)).\hfill \end{array}\end{array}$$

(18)

The cost function in Equation (18) is a random variable, since it depends on the uncertain TSO and DSO state variables. Therefore, we consider $J(x_k,{\mathsf{U}}_k):=\mathbb{E}\left[\tilde{J}(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k})\right]$ to obtain a deterministic cost function, which can be approximated empirically following the approach in [31] by averaging the value of its argument for some number $S_0$ of different realizations of the uncertain variable, which plays a tuning parameter role, i.e., $J(x_k,{\mathsf{U}}_k):=\frac{1}{S_0}{\sum }_{s=1}^{S_0}\mathbb{E}\left[\tilde{J}(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k}^s)\right]$.

The discretized dynamics in Equation (9) are non-deterministic, due to the uncertainty in the wind, reserve, and flexibility power. To obtain a set of deterministic dynamics, consider now the special case when $P_w=P_w^f$, and denote the corresponding TSO, DSO, building, and storage dynamics by $x_t^f$, $x_d^{f,i}$, $x_b^f$, $x_s^f$. In this case, $\mathrm{\Delta }{P}_w=0$, which means that the optimal reserve and flexibility contributions as determined in Equation (16) are both zero. Using this notation, we are in a position to formulate a finite-horizon stochastic control problem for each sampling time k using the following optimization program:

$$\begin{array}{c}\hfill \underset{{\mathsf{U}}_k}{minimize}J(x_k,{\mathsf{U}}_k)\mathrm{subject}\mathrm{to},\forall \ell \in {\mathcal{N}}_h:\end{array}$$

• TSO deterministic frequency model dynamics:

  $$\begin{array}{c}\hfill x_t^f(\ell +1|k)=f_t\left(x_t^f(\ell |k),P_{\mathrm{GR}}(\ell |k),P_{\mathrm{LD}}(\ell |k),P_w^f(\ell |k)\right).\end{array}$$

  (19a)
• TSO generation, ramping, line, and balance constraints:

  $$\begin{array}{c}\begin{array}{cc}\hfill & P_{\mathrm{GR}}^{min}\le P_{\mathrm{GR}}(\ell |k)\le P_{\mathrm{GR}}^{max},\hfill \\ \hfill & P_{\mathrm{GR}}^{\mathrm{down}}\le P_{\mathrm{GR}}(\ell +1|k)-P_{\mathrm{GR}}(\ell |k)\le P_{\mathrm{GR}}^{\mathrm{up}},\hfill \\ \hfill & L^{\min }\le {\tilde{L}}^{\mathrm{TSO}}\left[x_t^f(\ell +1|k)\right]\le L^{\max },\hfill \\ \hfill & {\left[\Gamma P_{\mathrm{GR}}(\ell |k)+P_w^f(\ell |k)-P_{\mathrm{BL}}(\ell |k)\right]}^{\top }{\mathbf{1}}_{n_t}=0,\hfill \end{array}\hfill \end{array}$$

  (19b)

  where $P_{\mathrm{GR}}^{\mathrm{down}}$ and $P_{\mathrm{GR}}^{\mathrm{up}}$ are ramping limits of power generation units, ${\tilde{L}}^{\mathrm{TSO}}$ is defined in Remark 1, and $L^{\min }$ and $L^{\max }$ are the TSO power flow limits.
• DSO deterministic frequency dynamics, $\forall i\in \mathcal{D}$:

  $$\begin{array}{cc}\hfill & x_d^{f,i}(\ell +1|k)=f_d^i\left(x_d^{f,i}(\ell |k),P_{\mathrm{IMP}}^i(\ell |k),P_{\mathrm{BD}}^i(\ell |k)\right).\hfill \end{array}$$

  (19c)
• DSO power line and balance constraints, $\forall i\in \mathcal{D}$:

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill & L^{i,\min }\le {\tilde{L}}^{\mathrm{DSO},\mathrm{i}}\left[x_d^{f,i}(\ell +1|k)\right]\le L^{i,\max },\hfill \\ \hfill & {\left[P_{\mathrm{IMP}}^i{(\ell |k)}^{\top }{\mathbf{1}}_{n_t}-P_{\mathrm{BD}}^i(\ell |k)\right]}^{\top }{\mathbf{1}}_{n_d^i}=0,\hfill \end{array}\end{array}$$

  (19d)

  where ${\tilde{L}}^{\mathrm{DSO},\mathrm{i}}$ is defined in Remark 1, and $L^{i,\min }$ and $L^{i,\max }$ are the DSO power flow limits.
• Buildings deterministic thermal comfort dynamics:

  $$\begin{array}{c}\hfill x_b^f(\ell +1|k)=f_b(x_b^f(\ell |k),P_{\mathrm{hvac}}(\ell |k)).\end{array}$$

  (19e)
• Buildings thermal comfort, power balance, and HVAC usage constraints:

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill & x_b^{min}\le x_b^f(\ell +1|k)\le x_b^{max},\hfill \\ \hfill & {\sum }_{i=1}^{n_d}{P}_{\mathrm{BD}}^i(\ell |k)=P_{\mathrm{hvac}}(\ell |k)+P_{\mathrm{stor}}(\ell |k)+P_{\mathrm{misc}}(\ell |k),\hfill \\ \hfill & P_{\mathrm{hvac}}^{min}\le P_{\mathrm{hvac}}(\ell |k)\le P_{\mathrm{hvac}}^{max},\hfill \end{array}\end{array}$$

  (19f)

  where $x_b^{min}$ and $x_b^{max}$ are the minimum and maximum bounds for the desired thermal comfort of buildings.
• Buildings electrical storage unit dynamics:

  $$\begin{array}{c}\hfill x_s^f(\ell +1|k)=\mathrm{\Xi }{x}_s^f(\ell |k)+\mathrm{\Omega }\left(P_{\mathrm{stor}}(\ell |k)\right),\end{array}$$

  (19g)
• Buildings storage capacity and power usage constraints:

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill & x_s^{min}\le x_s^f(\ell +1|k)\le x_s^{max},\hfill \\ \hfill & P_{\mathrm{stor}}^{min}\le P_{\mathrm{stor}}(\ell |k)\le P_{\mathrm{stor}}^{max},\hfill \end{array}\end{array}$$

  (19h)
• Probabilistic constraint:

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbb{P}\{{\mathsf{P}}_{w_k}\in \mathcal{W}|& L^{\min }\le {\tilde{L}}^{\mathrm{TSO}}\left[x_t(\ell +1|k)\right]\le L^{\max },\hfill \\ \hfill & L^{i,\min }\le {\tilde{L}}^{\mathrm{DSO},\mathrm{i}}\left[x_d^i(\ell +1|k)\right]\le L^{i,\max },\hfill \\ \hfill & P_{\mathrm{GR}}^{min}\le P_{\mathrm{GR}}(\ell |k)+R(\ell |k)\le P_{\mathrm{GR}}^{max},\hfill \\ \hfill & {\left[S_{\mathrm{IMP}}^i{(\ell |k)}^{\top }{\mathbf{1}}_{n_t}-S_{\mathrm{BD}}^i(\ell |k)\right]}^{\top }{\mathbf{1}}_{n_d^i}=0,\hfill \\ \hfill & {\sum }_{i=1}^{n_d}{S}_{\mathrm{BD}}^i(\ell |k)=S_{\mathrm{hvac}}(\ell |k)+S_{\mathrm{stor}}(\ell |k),\hfill \\ \hfill & P_{\mathrm{hvac}}^{min}\le P_{\mathrm{hvac}}(\ell |k)+S_{\mathrm{hvac}}(\ell |k)\le P_{\mathrm{hvac}}^{max},\hfill \\ \hfill & P_{\mathrm{stor}}^{min}\le P_{\mathrm{stor}}(\ell |k)+S_{\mathrm{stor}}(\ell |k)\le P_{\mathrm{stor}}^{max},\hfill \\ \hfill & -a(\ell |k)\le b(\ell \left|k\right)\le c(\ell |k),\hfill \\ \hfill & \forall i\in \mathcal{D}\mathrm{and}\forall \ell \in {\mathcal{N}}_h\}\ge 1-\epsilon ,\hfill \end{array}\end{array}$$

  (19i)

  where $x_t\left(·\right)$ and $x_d^i\left(·\right)$ are given by Equation (9), $\epsilon \in (0,1)$ is the level of admissible constraint violation, $a={\sum }_{i\in \mathcal{G}}{R}^{\mathrm{us}}(\ell |k)+{\sum }_{l\in \mathcal{B}}(S_{\mathrm{hvac}}^{\mathrm{dd},\mathrm{l}}(\ell |k)+S_{\mathrm{stor}}^{\mathrm{dd},\mathrm{l}}(\ell |k))$, $b={\sum }_{m\in \mathcal{F}}(P_w(\ell |k)-P_w^f(\ell |k))$, and $c={\sum }_{i\in \mathcal{G}}{R}^{\mathrm{ds}}(\ell |k)+{\sum }_{l\in \mathcal{B}}(S_{\mathrm{hvac}}^{\mathrm{id},\mathrm{l}}(\ell |k)+S_{\mathrm{stor}}^{\mathrm{id},\mathrm{l}}(\ell |k))$, such that the formulation for $S_{\mathrm{hvac}}^{\mathrm{id}}(\ell |k)$, $S_{\mathrm{hvac}}^{\mathrm{dd}}(\ell |k)$, $S_{\mathrm{stor}}^{\mathrm{id}}(\ell |k)$, and $S_{\mathrm{stor}}^{\mathrm{dd}}(\ell |k)$ are given in Section 3.2.

The constraint Equation (19i) is a chance constraint, which ensures that the feasibility probability for the constraint is above the specified level. Using the proposed chance constraint Equation (19i), stochastic MPC offers an alternative approach to robust MPC. Since stochastic MPC directly incorporates the trade-off between constraint feasibility and control performance, the corresponding solutions are generally less conservative.

Please note that Equations (19a)–(19h) are all based on the deterministic BtG dynamics for the special case when no ancillary services are required, i.e., $\mathrm{\Delta }{P}_w$, R, $S_{\mathrm{LD}}$, $S_{\mathrm{IMP}},S_{\mathrm{BD}}$, $S_{\mathrm{hvac}}$, and $S_{\mathrm{stor}}$ are all zero. On the contrary, the probabilistic constraint Equation (19i) depends on the non-deterministic dynamics. Hence, the optimization problem in Equation (19) is a finite-horizon, chance-constrained quadratic program, whose stages are coupled by the dynamics of the TSO, DSO, buildings thermal comfort, and storage systems at each sampling time k. We note that the proposed optimization program in Equation (19) is in general a non-convex problem which is hard to solve. The feasible set of Equation (19) is non-convex and hard to determine explicitly in the presence of chance constraints Equation (19i). In the following section, we will develop a tractable framework to obtain a probabilistically feasible solution.

5. Tractable Robust MPC Reformulation

Using a more compact notation, the chance constraints Equation (19i) can be written as

$$\mathbb{P}\{{\mathsf{P}}_{w_k}\in \mathcal{W}|g(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k})\le 0\}\ge 1-\epsilon .$$

(20)

To this end, we approximate the proposed probabilistic constraints in Equation (20) using the following robust reformulation:

$$g(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k})\le 0,\forall {\mathsf{P}}_{w_k}\in \mathbb{W},$$

(21)

where $\mathbb{W}$ is a bounded uncertainty set compared to $\mathcal{W}$ which is an unbounded and unknown uncertainty set in Equation (20). In the perspective of power systems, the proposed tractable reformulation Equation (21) translates the unbounded set of wind-power scenarios $\mathcal{W}$ to a bounded set $\mathbb{W}$ which contains the more probable scenarios of wind power with a high level of confidence. Please note that the proposed constraint Equation (21) is a robust constraint, since it should be satisfied for all ${\mathsf{P}}_{w_k}\in \mathbb{W}$. Obviously, for any feasible solution for Equation (21), ${\mathsf{P}}_{w_k}\in \mathbb{W}$ implies $g(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k})\le 0,{\mathsf{P}}_{w_k}\in \mathcal{W}$. Therefore, by choosing $\mathbb{W}$ that covers at least a $(1-\epsilon )$ content of ${\mathsf{P}}_{w_k}$, i.e., $\mathbb{W}$ satisfies $\mathbb{P}\{{\mathsf{P}}_{w_k}\in \mathbb{W}\}\ge 1-\epsilon$, any feasible solution for Equation (21) must satisfy $\mathbb{P}\{{\mathsf{P}}_{w_k}\in \mathcal{W}|g(x_k,{\mathsf{U}}_k,{\mathsf{P}}_{w_k})\le 0\}\ge \mathbb{P}\{{\mathsf{P}}_{w_k}\in \mathbb{W}\}\ge 1-\epsilon$, implying that it is also feasible for Equation (20) (see [32] for more detailed descriptions). In other words,

Lemma 1.

Any feasible solution for Equation (21) using a $(1-\epsilon )$-content set $\mathbb{W}$ is a feasible solution for the probabilistic constraint in Equation (20) with probability $1-\epsilon$.

Proof. 

The proof is straightforward, and we omit it for the sake of brevity. □

Let us now introduce $\mathbb{W}\subseteq {\mathbb{R}}^{n_w=1}$ as a bounded set. We assume for simplicity that $\mathbb{W}$ is an axis-aligned hyper-rectangular set. Please note that this is not a restrictive assumption and any convex set, e.g., ellipsoids and polytopes, could have been chosen instead as described in [31]. We can define $\mathbb{W}:=[-\omega ,\omega ]$ as an interval, where the vector $\omega \in {\mathbb{R}}^{n_w}$ defines the hyper-rectangle bounds.

Consider now the following optimization problem that aims to determine the set $\mathbb{W}$ with minimal volume:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\underset{\omega \in {\mathbb{R}}^{n_w}}{min}\hfill & {\parallel \omega \parallel }_1\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathsf{P}}_{w_k}^l\in [-\omega ,\omega ],\forall l=1,\dots ,S\hfill \end{array}\right.,\end{array}$$

(22)

where S is a finite number of i.i.d samples ${\mathsf{P}}_{w_k}^l\in \mathcal{W}$ needed to be available from either a known distribution $\mathcal{W}$ or through historical observations. If we denote by $\mathbb{W}=[-\tilde{\omega },\tilde{\omega }]$ the optimal solution of Equation (22), then the following proposition provides an explicit relation between the number of required samples S and the $(1-\epsilon )$-content set $\mathbb{W}$.

Proposition 1.

Fix $\epsilon \in (0,1)$, $\beta \in (0,1)$, determine $S\ge [\frac{2}{\epsilon }(n_w+ln\frac{1}{\beta })]$, and solve Equation (22) to obtain its optimal solution $[-\tilde{\omega },\tilde{\omega }]=\mathbb{W}$. Then, with probability of at least $1-\beta$,

$$\mathbb{P}\left\{{\mathsf{P}}_{w_k}\in \mathcal{W}:{\mathsf{P}}_{w_k}\notin [-\tilde{\omega },\tilde{\omega }]=\mathbb{W}\right\}\le \epsilon .$$

(23)

Proof. 

The proof is based on the results in [33] by noting that the proposed optimization problem Equation (22) is a convex program with the number of decision variables equivalent to the number of uncertain variables $n_w=1$. We therefore omit the proof for the sake of brevity. □

Using the proposed robust reformulation in Equation (21), one can obtain a tractable formulation similar to (Proposition 1, [34]). The obtained solution of the optimization program in Equation (19), where Equation (19i) is replaced with Equation (21), is the optimal input sequence ${\mathsf{U}}_k^{*}$. Following the MPC paradigm, only the first element of the optimal control input is implemented at time k, and we proceed based on the receding horizon principle. This means that at every time step k, the robust reformulation of Equation (19) is solved using the current measurement of the state variables $\{x_t\left(k\right),x_d\left(k\right),x_b\left(k\right),x_s\left(k\right)\}$.

6. Numerical Case Study

To demonstrate the grid regulative capacity of (a) reserve scheduling and (b) building HVAC and storage flexibility on the BtG system under wind-power penetration, we consider the following three cases in the simulation study: (1) BtG system with only reserve power enabled, (2) BtG system with only flexibility power enabled, (3) BtG system with both reserve and flexibility power enabled. Please note that we consider the same cost coefficients for both reserve and flexibility deployments in the cost function Equation (18). Throughout the simulation horizon, the wind-power penetration is around 10–15% of the total electricity generation in all cases.

6.1. Simulation Setup

Consider a network of 1 TSO with 3 generators (GRs), 2 DSOs, 13 buildings, and 1 wind farm (WF), as visualized in Figure 3. The case studies are simulated for $N_{\mathrm{sim}}=24$ h with the prediction horizon of an hour and time resolution $h=5$ min. Thus, the length of the prediction horizon is $N_h=\frac{60}{5}=12$ steps. To generate the number of required scenarios of wind power, we use the Markov chain model in [30]. Following Proposition 1, we choose $\epsilon =0.05$, $\beta ={10}^{-4}$, and $S=1328\ge [\frac{2}{ϵ}(2\times 12+ln\frac{1}{\beta })]$, and then solve Equation (22) to obtain the bounded set $\mathbb{W}$ at each sampling time k. The associated cost parameters are $Q_t=Q_d=1000{\mathrm{rad}}^{-2}$, power generation and consumption are valued at $Q_{\mathrm{GR}}=Q_{\mathrm{hvac}}=0.1 \mathrm{M}{\mathrm{W}}^{-2}$, and reserve and flexibility scheduling at $q_R=q_S=1000 \mathrm{M}{\mathrm{W}}^{-1}.$ Please note that to have a fair comparison, we assumed the same cost coefficients for both reserve and flexibility. Since our simulation study considers large-scale buildings, the limits for each storage unit are set equivalent to four Tesla Model S batteries [35], i.e., $P_{\mathrm{stor}}^{\min }=-2 \mathrm{M}\mathrm{W}$, $P_{\mathrm{stor}}^{\max }=2 \mathrm{M}\mathrm{W}$, $x_s^{\min }=2 \mathrm{k}\mathrm{W} \mathrm{h}$, and $x_s^{\max }=400 \mathrm{k}\mathrm{W} \mathrm{h}$. Simulations were implemented in MATLAB, with Yalmip as interface [36] and Gurobi as solver [37]. Building parameters and miscellaneous (uncontrollable) load profiles for the large-scale, commercial buildings were adopted from [17]. Grid parameters (line length, cable type, line susceptance, etc.) were obtained from the MatPower IEEE 5-bus power system [38], with some parameters (e.g., power line limits) scaled appropriately for the current simulation studies. For further details on the building and grid parameters, we refer the reader to the sources above. We also carried out Monte Carlo simulations with the system under optimal control input simulated for 10,000 different wind-power scenarios. We obtain a-posteriori the violation probability (empirical violation level), by means of counting the number of trajectories that violate any of the constraints, to check if the theoretical maximum violation level is indeed satisfied.

6.2. Simulation Results

The power balance in the TSO network for all cases is presented in Figure 4, showing a reduced production for Case 2 and 3, compared to Case 1. In particular, Case 2 yields a total reduction in conventional generation of 2.7% compared to Case 1, whereas for Case 3, this reduction is 12.8%. The solid black line shows the forecast wind power, and the dark green bars show the actual wind power. As shown in both Figure 4 and Figure 6, between 3:00–11:00 and 15:00–17:00, the wind-power error is negative, $\mathrm{\Delta }{P}_w\le 0$, and from 12:00–14:00 and 18:00–24:00, the error is positive, $\mathrm{\Delta }{P}_w\ge 0$. To restore the power balance, when $\mathrm{\Delta }{P}_w\le 0$, Case 2 employs the decreased-demand flexibility and on the contrary, Case 1 uses the up-spinning reserve power to reduce output power dispatch. Similarly, when $\mathrm{\Delta }{P}_w\ge 0$, Case 2 employs the increased-demand flexibility and on the contrary, Case 1 uses the down-spinning reserve power to reduce output power dispatch.

Figure 5 depicts the frequency deviations in the TSO and DSO networks for Case 3. Frequency deviations for all the cases are kept well below allowable limits. Enabling building-side flexibility results in a reduction of 5.3% in cumulative squared frequency deviation for the total network, and a reduction of 32.5% in the DSO network only.

The wind-power error together with the reserve and flexibility dispatch for Case 3 is presented in Figure 6. As reserve and flexibility are valued at the same cost in the objective function, the distribution between the two indicates that both mechanisms are cooperating well to meet the power balance. The total flexibility dispatch accounts for 44.37% of the total compensation for the wind-power error over the 24h simulation time. It is shown that in the case of a positive wind-power error, flexibility dispatch is promoted over reserve, while for a negative wind-power error, the opposite is observed. Finally, Figure 6 also shows that the wind-power error is perfectly compensated throughout the full day, by the sum of flexibility and reserve deployment.

The reserve and flexibility scheduling is compared to the actual dispatch in more detail in Figure 7. In line with Figure 6, up-spinning reserve is promoted over decreased-demand flexibility, while increased-demand is favored over down-spinning reserve. Although total power demand differs by less than 10% between all cases, this motivates why the power dispatch of the cases is of different shape, as observed in Figure 4. As shown in Figure 6, it is important to note that the reserve and flexibility scheduling capacity is always significantly higher than the actual reserve or flexibility dispatch (i.e., the portion of the capacity that is deployed in the operational decision-making). This suggests that the current reserve and flexibility scheduling capacity is also sufficient for dealing with higher wind-power penetration levels than the current 10%.

The results of our Monte Carlo simulation to determine the empirical constraint violation level are presented in Figure 8 for all 3 cases along with a case in which no reserve or flexibility is used. As the majority of violations are observed during peak load hours, only this portion of the simulation horizon is presented. For the case with no reserve or flexibility, empirical violation levels are extremely high, as any deviation of the wind power from the forecast trajectory almost exclusively results in violation of the power line limits. For Case 1 and 3, maximum empirical violation levels are around 0.1%, whereas for Case 2, no violations were observed in any of the 10,000 iterations. In all three cases, the empirical violation level is well below the theoretical limit of $\epsilon =0.05$.

7. Conclusions

In this paper, we have presented two new developments for the BtG integration framework. First, we extended the existing models by introducing uncertain wind-power generation and by explicitly formulating the interactions between TSO, DSOs, and buildings. Second, we developed a unified BtG framework to handle uncertain generation, by integrating demand-side flexibility provided by individual buildings into the traditional process of reserve scheduling. We provided explicit expressions to determine the available amount of building demand-side flexibility from HVAC and electrical storage units. Using the unified BtG model, we formulated a finite-horizon stochastic control problem and provided a tractable robust reformulation with probabilistic feasibility guarantees. As the main outcome of our proposed BtG framework, we conclude that the demand-side flexibility can substitute the traditional reserve scheduling services in power systems in the presence of wind-power generation, without losing stability properties of the power grid and violating the buildings thermal comfort of occupants.

At the same time, the proposed framework is still rather idealistic, and relies on multiple major assumptions. As a final note, we discuss some of these assumptions, and provide possible directions for future research:

• First of all, we did not consider the impact of imperfect communication between the TSO, DSOs, and buildings, and other sources of uncertainty than the wind power, such as demand uncertainty, were left out of scope. Hence, in a more sophisticated BtG integration framework, multiple sources of uncertainty should be incorporated.
• Second, the BtG framework in the current work is formulated as a centralized MPC framework. Although the current centralized implementation runs in reasonable time for hundreds of buildings, we believe that it is worth exploring decentralized control frameworks, e.g., [39], in order to reduce the computational complexity of the problem.
• Finally, part of our current work focuses on integrating the psychological impact of end-users for participating in the ancillary service market, by providing demand-side flexibility to the grid. We are interested in constructing models to simulate the willingness of end-users to participate in the ancillary service market, in order to study how different stimuli can influence the psychological behavior of consumers.