Buildings-to-Grid Integration Framework

TL;DR

This paper introduces a mathematical Buildings-to-Grid integration framework that optimizes combined building and power grid operations, improving efficiency, reducing costs, and providing effective frequency regulation in smart city networks.

Contribution

It develops a novel MPC-based framework that explicitly couples building and grid control actions, accounting for different time-scales and dynamic constraints, which was not previously addressed.

Findings

Reduces total system costs compared to decoupled control strategies.

Achieves building energy savings and frequency regulation in simulations.

Maintains performance under weather and load forecast uncertainties.

Abstract

This paper puts forth a mathematical framework for Buildings-to-Grid (BtG) integration in smart cities. The framework explicitly couples power grid and building's control actions and operational decisions, and can be utilized by buildings and power grids operators to simultaneously optimize their performance. Simplified dynamics of building clusters and building-integrated power networks with algebraic equations are presented---both operating at different time-scales. A model predictive control (MPC)-based algorithm that formulates the BtG integration and accounts for the time-scale discrepancy is developed. The formulation captures dynamic and algebraic power flow constraints of power networks and is shown to be numerically advantageous. The paper analytically establishes that the BtG integration yields a reduced total system cost in comparison with decoupled designs where grid and building operators determine their controls separately. The developed framework is tested on standard power networks that include thousands of buildings modeled using industrial data. Case studies demonstrate building energy savings and significant frequency regulation, while these findings carry over in network simulations with nonlinear power flows and mismatch in building model parameters. Finally, simulations indicate that the performance does not significantly worsen when there is uncertainty in the forecasted weather and base load conditions.

Figures

Tables

Table 1. TABLE I: Building-Integrated Power Network Setup.

	Case 9	Case 14	Case 30	Case 57
Number of Buses	9	14	30	57
Number of Generators	3	5	6	7
Number of Buildings	965	1058	376	1822

Table 2. TABLE II: Cost Comparison in 1,000$ for the Decoupled Problems (Scenarios I and II) versus BtG-GMPC.

Test Case	Cost Function	Scenario I	Scenario II	BtG-GMPC
Case 9	$\sum {𝐱}_g^{\top }{\mathrm{𝐐𝐱}}_g$	490.51	205.57	20.33
	$\sum \Delta {𝐮}_g^{\top }𝐑\Delta {𝐮}_g$	60.00	58.95	59.39
	$\sum J\left({\overline{𝐮}}_g\right)$	125.08	121.05	120.26
	$\sum {𝐜}_b^{\top }{𝐮}_b$	656.73	548.40	548.56
	Total Grid Cost	675.59	385.56	199.98
	Total Cost	1332.33	933.96	748.54
Case 14	$\sum {𝐱}_g^{\top }{\mathrm{𝐐𝐱}}_g$	289.19	134.75	3.02
	$\sum \Delta {𝐮}_g^{\top }𝐑\Delta {𝐮}_g$	70.49	69.60	70.39
	$\sum J\left({\overline{𝐮}}_g\right)$	81.32	77.61	77.22
	$\sum {𝐜}_b^{\top }{𝐮}_b$	705.48	588.12	588.24
	Total Grid Cost	441.01	281.96	150.62
	Total Cost	1146.49	870.08	738.86
Case 30	$\sum {𝐱}_g^{\top }{\mathrm{𝐐𝐱}}_g$	49.49	20.69	0.02
	$\sum \Delta {𝐮}_g^{\top }𝐑\Delta {𝐮}_g$	3.19	3.15	3.12
	$\sum J\left({\overline{𝐮}}_g\right)$	1.44	1.31	1.33
	$\sum {𝐜}_b^{\top }{𝐮}_b$	253.02	211.30	211.31
	Total Grid Cost	54.12	25.14	4.46
	Total Cost	307.14	236.45	215.78
Case 57	$\sum {𝐱}_g^{\top }{\mathrm{𝐐𝐱}}_g$	466.91	25.75	0.09
	$\sum \Delta {𝐮}_g^{\top }𝐑\Delta {𝐮}_g$	159.61	156.95	155.40
	$\sum J\left({\overline{𝐮}}_g\right)$	49.31	46.40	46.48
	$\sum {𝐜}_b^{\top }{𝐮}_b$	1230.20	1027.18	1027.18
	Total Grid Cost	675.83	229.10	201.97
	Total Cost	1906.03	1256.28	1229.15

Table 3. TABLE III: Computational Time for Running BtG-GMPC and Algorithm  1 .

	Case 9	Case 14	Case 30	Case 57
Computational Time (Hours)	0.26	0.50	0.82	3.74

Equations

$${\dot{T}}_{\mathrm{wall}}$$

$${\dot{T}}_{\mathrm{zone}}$$

$${\dot{\mathbf{x}}}_b^{(l)}(t)={\mathbf{A}}_b^{(l)}{\mathbf{x}}_b^{(l)}(t)+{\mathbf{B}}_{u_b}^{(l)}{\mathbf{u}}_b^{(l)}(t)+{\mathbf{B}}_{w_b}^{(l)}{\mathbf{w}}_b^{(l)}(t),$$

$${\mathbf{A}}_b^{(l)}$$

$${\mathbf{B}}_{u_b}^{(l)}$$

$${\dot{\mathbf{x}}}_b(t)={\mathbf{A}}_b{\mathbf{x}}_b(t)+{\mathbf{B}}_{u_b}{\mathbf{u}}_b(t)+{\mathbf{B}}_{w_b}{\mathbf{w}}_b(t),$$

$${\gamma }_{k,m}$$

$${\pi }_{k,l}$$

$$M_k{\ddot{\delta }}_k(t)+D_k{\dot{\delta }}_k(t)$$

$$P_{L_k}(t)=P_{{\mathrm{BL}}_k}(t)+D_k^{\prime }{\dot{\delta }}_k(t)+\sum _{l=1}^{n_b}{\pi }_{k,l}{P}_{\mathrm{bldg}}^{(l)}(t).$$

$$M_k{\ddot{\delta }}_k(t)$$

$${\dot{\delta }}_k(t)$$

$$M_k{\dot{\omega }}_k(t)$$

$${\mathbf{E}}_g{\dot{\mathbf{x}}}_g(t)$$

$${\mathbf{E}}_g$$

$$\mathbf{M}$$

$${\Phi }_k$$

$$\mathbf{D}$$

$${\mathbf{B}}_{w_g}$$

$$\text{LOPF:}\underset{{\bar{\mathbf{u}}}_g=\{{\bar{u}}_{g_i}{\}}_{i=1}^{n_g}}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$${\mathbf{x}}_g(k_g{h}_g)$$

$${\mathbf{x}}_b(k_b{h}_b)$$

$$f(\cdot )$$

$$\underset{\begin{array}{c}{\mathbf{U}}_b,\Delta {\mathbf{U}}_g,{\bar{\mathbf{u}}}_g\\ {\mathbf{X}}_b,{\mathbf{X}}_g\end{array}}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$$\Delta {\mathbf{u}}_g^{\min }\le \Delta {\mathbf{u}}_g(t+k_g{h}_g)\le \Delta {\mathbf{u}}_g^{\max }$$

$${\mathbf{x}}_g^{\min }\le {\mathbf{x}}_g(t+k_g{h}_g)\le {\mathbf{x}}_g^{\max }$$

$$\forall k_g\in \{1,\dots ,T_p/h_g\}$$

$${\mathbf{x}}_b(t+k_b{h}_b)={\mathbf{f}}_b({\mathbf{x}}_b,{\mathbf{u}}_b,{\hat{\mathbf{w}}}_b|t,s)$$

$${\mathbf{u}}_b^{\min }\le {\mathbf{u}}_b(t+k_b{h}_b)\le {\mathbf{u}}_b^{\max }$$

$${\mathbf{x}}_b^{\min }\le {\mathbf{x}}_b(t+k_b{h}_b)\le {\mathbf{x}}_b^{\max }$$

$$\forall k_b\in \{1,\dots ,T_p/h_b\}$$

$$\text{eqref}equ:DCOPF1,\text{eqref}equ:DCOPF3$$

$${\mathbf{u}}_b(t+k_g{h}_g)={\bar{\mathbf{u}}}_b={\mathbf{u}}_b(t+k_b{h}_b)$$

$$\forall k_g{h}_g\in [k_b{h}_b,(k_b+1)h_b).$$

$$f_{\mathrm{BtG}}^{\ast }=\underset{\begin{array}{c}{\bar{\mathbf{u}}}_g,\Delta {\mathbf{U}}_g,{\mathbf{X}}_g,\\ {\mathbf{U}}_b,{\mathbf{X}}_b\end{array}}{\mathrm{minimize}}$$

$$subject\; to$$

$$({\bar{\mathbf{u}}}_g,\Delta {\mathbf{U}}_g,{\mathbf{X}}_g)\in {\mathcal{F}}_g$$

$$({\mathbf{U}}_b,{\mathbf{X}}_b)\in {\mathcal{F}}_b.$$

$$f_b^{\mathrm{MPC}}=\underset{{\mathbf{U}}_b,{\mathbf{X}}_b}{\mathrm{minimize}}$$

$$subject\; to$$

Full text

I Introduction, Prior Art, Paper Contributions

By 2050, a staggering 70% of the world’s population is bound to live and work in cities [1]. A recent assessment from the World Bank suggested that two-thirds of global energy consumption can be attributed to cities, leading to 71% of global direct energy-related greenhouse gas emissions [2]. Smart cities consist of sustainable and resilient infrastructures, where buildings are a major constituent. Building energy consumption contributes to more than 70% of electricity usage—profoundly impacting power grid’s operation. Futuristic cities equipped with optimized building designs have the auspicious potential to play a pivotal role in reducing global energy consumption while maintaining stable electric-grid operations. As buildings are physically connected to the electric power grid, it is natural to understand their coupling and develop a framework for Buildings-to-Grid (BtG) integration. To understand the role and impact of BtG integration, the authors in [3] provide relevant research questions for BtG integration.

The installation of smart meters in buildings and across the power grids enables the BtG integration, which can transform passive buildings into active dispatchable demand resources. The U.S. Department of Energy has highlighted the multiple benefits and opportunities of BtG integration [4]: 1) Buildings can enjoy significant energy savings; 2) the grid’s resources are more efficiently utilized, as peak demand is curbed; 3) the grid can become more stable with fewer frequency excursions; 4) the need for bulk generation and transmission investments is deferred; and 5) with a BtG integration platform in place, distributed energy resources at the buildings’ premises, such as photovoltaic units and electric vehicles, can be more efficiently integrated with the power grid, turning into significant assets. Besides the aforementioned technical motivating factors, the mathematical intuition behind the BtG integration is that when the grid and the buildings jointly optimize their control decisions, they have the potential to yield larger system benefits than when they make these decisions separately. This paper aspires to develop a BtG integration framework, and make this mathematical intuition precise.

Various studies address a breadth of computational and experimental aspects of BtG integration. An overview of demand response potential from smart buildings is presented in [5]. An experimental architecture that enables smart buildings is proposed in [6] with a focus on heating, ventilation, and air conditioning (HVAC) systems and grid integration. A bi-level optimization framework for commercial buildings integrated with a distribution grid is proposed in [7]. Detailed dynamical models for buildings with multiple zones (upper level) and an operational model for the distribution grid with voltage and current balance equations (lower level) are included; nonetheless, a dynamical model of the power grid, suitable for modeling frequency excursions, is missing.

The regulation service provision by smart buildings is investigated in [8], where price signals are exchanged between grid and building operators to alter building energy consumption. Other BtG integration studies have shown that grid-aware building HVAC controls can provide frequency regulation or other ancillary services to the grid [9, 10, 11, 12, 13], largely without sacrificing the occupants’ comfort. The load-shifting capability of buildings has also been explored [14]. Explicit account of the grid dynamics and power flows is on the other hand missing from the previously mentioned works.

HVAC controls and building dynamics are typically modeled as linearized dynamical systems [15]. The modeling particulars depend on the size of the building cluster. As the number of buildings involved in the analysis increase, the dynamical models tend to become simpler—for obvious computational purposes. A typical thermal resistance and capacitance circuit model can be used to represent heat transfer and thermodynamical properties of the building envelope, and is widely used in building control studies [16, 17, 18, 19]. Given these models for building dynamics, various control routines have been developed for building controls. Currently, many commercial buildings use PID controllers for HVAC systems [20]. However, model predictive control (MPC) has proved to be advantageous with respect to PID controllers [21, 22, 23].

Building HVAC control via MPC has been investigated under a wide range of scenarios and setups. For example, uncertainty in the building MPC formulation is considered in [24, 25, 26]. Centralized building MPC routines are proposed in [27, 28], and the works in [21, 22, 29] investigate decentralized or distributed solvers to building MPC problems. In addition, explicit MPC routines have also been developed in the context of building control studies [30]. Other works focus on integrating occupancy behavior and its impact on indoor temperature variations, while still attempting to obtain optimal control laws [23, 31]. The majority of the aforementioned works show significant energy savings given different system dynamics, forecast and parametric uncertainty, and computational limitations.

While the aforementioned research investigated different challenging problems related to BtG integration and building MPC routines, none of these studies produces a high-level mathematical framework that buildings and power grids operators can simultaneously utilize to optimize their performance—a framework that explicitly couples power grid and building control actions and operational decisions. In addition, the majority of the previously mentioned studies focus on one or a group of buildings and the corresponding impact on the power grid, rather than clusters of thousands of buildings in smart cities. The main challenges associated with creating such a framework that addresses the aforementioned research gaps are as follows.

• •

Building control systems are neither connected to each other, nor integrated with the grid. Consequently, a unified optimal energy control strategy—even if it is decentralized—cannot be achieved unless there is a framework that facilitates this integration, in addition to the willingness of building operators to contribute to this framework.

• •

Grid and building dynamics and control actions clearly operate at two different time-scales. While the grid controls and states are often in seconds, the building state dynamics and controls are much slower, often in minutes. Coupling the two dynamic systems together entails addressing this time-scale discrepancy.

• •

Existence of algebraic equations in grid dynamics, resulting in differential algebraic equations (DAE), coupled with the different time-scales, complicates modeling and analysis of BtG integration. In fact, these algebraic equations depict the interdependence between grid and building dynamics.

The chief contribution of this paper is a novel mathematical framework for BtG integration that addresses the aforementioned challenges in a structured and principled way. The MPC-based framework couples building dynamics, grid dynamics that include the network frequency, and the power flow equations. The objective is to generate local control actions for buildings and power generators such that the overall performance is optimized in terms of stability, energy savings, and other socio-economic metrics. The time scale discrepancy between the grid and building dynamics is explicitly accounted for in the developed MPC-based optimization formulation. The building-integrated power network dynamics are modeled by DAEs. In order to include the DAEs into the optimization framework, the DAEs are discretized using Gear’s method. While this discretization method has been the basis for the numerical solution of DAEs for several decades [32], it is the first time that it is brought in to facilitate the development of a BtG integration framework. The paper analytically establishes that the BtG integration yields a reduced total system cost with respect to decoupled designs in which grid and building operators determine their controls separately.

The developed framework is tested in standard IEEE networks that include hundreds to thousands of buildings modeled using ASHRAE data. The building HVAC load is driven by realistic ambient weather patterns and typical temperature requirements for commercial buildings. The simulations are performed for various power networks, and a reduction of up to 20% and 43% in total system cost—with respect to two other decoupled designs—is demonstrated. Finally, simulations indicate that the performance does not significantly worsen when there is uncertainty in the building parameters or forecasted weather and base load conditions.

The remainder of this paper is organized as follows. In Sections II and III, we present the dynamics of building clusters and of the building-integrated power network. Optimal power flow is also integrated in these models. In Section IV, we propose our approach for BtG integration, while addressing the aforesaid challenges. Then, the optimal control problem that models the BtG integration is formulated. A customized algorithm is also developed to seamlessly include the optimal power flow into the integrated framework. Section V produces an analytical discussion on the advantages of the BtG framework over decoupling the optimization of buildings and power grids. Case studies with realistic building parameters and grid constraints are given in Section VI. Finally, future work is outlined in Section VII.

II Building Clusters Dynamics

The patterns of energy usage in buildings are impacted by local climate, heat transfer through the building envelope, daily operation, and occupancy behaviors. Detailed energy models have been developed based on physics and statistics to simulate heat transfer in buildings.

For a large-scale application such as BtG integration, it is unrealistic to consider every thermal zone of each building—thousands of buildings will generate millions of zones. This approach would produce a highly intractable BtG integration problem. Hence, at a BtG integration level, the amount of cooling energy needed which is optimized to minimize the total operation cost is allocated to each building; we define this quantity as $P_{\mathrm{HVAC}}^{(l)}$ for building $l$. Then, at the local level of building $l$, the decision variables of the air-side system (setpoints for air-handling units, damper opening for terminal systems) and the water-side system (flows for pumps, chiller temperature setpoints) can be optimized to maintain the preferred zone temperature, while not exceeding the cooling load limits set by $P_{\mathrm{HVAC}}^{(l)}$. This approach of solving for the cooling loads and then feeding the setpoints to local lower level problems is common in recent building studies; see [22].

In this paper, we focus on the high-level commercial buildings problem described above. We use a typical thermal resistance and capacitance (RC) network to model heat transfer and the thermodynamics of the building envelope, which has been widely used in building control studies [16, 17, 18, 19]. The RC network model assumes a steady-state heat transfer through the building envelope. Considering that the building dynamics have time constant of hours, this model is sufficient for a high-level BtG integration study. A typical three-resistance and two-capacitance (3R-2C) model is shown in Fig. 1.

In this model, a building is treated as a super-zone where the resistance parameters represent the thermal resistance of the building structure, the external facades’ convection, and the internal walls’ convection. Building dynamics with temperatures $T_{\mathrm{wall}}(t)=T_{\mathrm{wall}}$ and $T_{\mathrm{zone}}(t)=T_{\mathrm{zone}}$ are written as

[TABLE]

where $R_{\mathrm{win}}$, $R_{2}$, and $R_{1}$ are physical parameters of the building envelope; $C$ is the lumped thermal capacity of all walls and the roof; $C_{\text{zone}}$ is the thermal capacity of the zone; $\dot{Q}_{\mathrm{sol}}(t)$ is the total absorbed solar radiation on the external wall; $\dot{Q}_{\mathrm{int}}(t)$ is the total internal heat gain from space heat sources such as desktops, people, and lights; $T_{\mathrm{zone}}(t)$ and $T_{\mathrm{wall}}(t)$ are respectively the zone (space) and wall temperatures; and $T_{\mathrm{amb}}(t)$ is the outside ambient temperature. The cooling load can be calculated as $\dot{Q}_{\mathrm{HVAC}}(t)=\mu_{\mathrm{HVAC}}P_{\mathrm{HVAC}}(t),$ where $P_{\mathrm{HVAC}}$ is the power consumed by the HVAC system, and $\mu_{\mathrm{HVAC}}$ is the coefficient of performance of the HVAC system. The dynamics of building $l$ are described by

[TABLE]

where $\mathbf{x}_{b}^{(l)}=[T_{\mathrm{wall}}\;T_{\mathrm{zone}}]^{\top}_{l}$ is the state of building $l$; $\mathbf{u}_{b}^{(l)}=\begin{bmatrix}P_{\mathrm{HVAC}}\end{bmatrix}_{l}$ is the control input variable; $\mathbf{w}_{b}^{(l)}=[T_{\mathrm{amb}}\;\dot{Q}_{\mathrm{sol}}\;\dot{Q}_{\mathrm{int}}]^{\top}_{l}$ is a random uncontrollable input. Various methods have been developed to provide an estimate of $\mathbf{w}_{b}^{(l)}$, denoted by $\hat{\mathbf{w}}_{b}^{(l)}$, for each building [34]. The system state-space matrices in (1) are defined as

[TABLE]

The notation $[\cdot]_{l}$ implies that each building $l$ has a different set of RC parameters. In this paper, we consider clusters of buildings with each cluster connected to a power grid node. Since we aim to understand the impact of buildings’ contribution to frequency regulation and overall energy consumption costs, we present the dynamics of building clusters

[TABLE]

where $n_{b}$ is the total number of buildings in the network; $\mathbf{x}_{b}\in\mathbb{R}^{2n_{b}},\mathbf{u}_{b}\in\mathbb{R}^{n_{b}},$ and $\mathbf{w}_{b}\in\mathbb{R}^{3n_{b}}$. In the absence of communication between buildings, the state-space matrices $\mathbf{A}_{b}$, $\mathbf{B}_{u_{b}}$, and $\mathbf{B}_{w_{b}}$ will all be block diagonal matrices: $\mathbf{A}_{b}=\operatorname{diag}(\mathbf{A}_{b}^{(1)},\ldots,\mathbf{A}_{b}^{(n_{b})})\in\mathbb{R}^{2n_{b}\times 2n_{b}},\;\mathbf{B}_{u_{b}}=\operatorname{diag}(\mathbf{B}_{u_{b}}^{(1)},\ldots,\mathbf{B}_{u_{b}}^{(n_{b})})\in\mathbb{R}^{2n_{b}\times n_{b}},\;\mathbf{B}_{w_{b}}=\operatorname{diag}(\mathbf{B}_{w_{b}}^{(1)},\ldots,\mathbf{B}_{w_{b}}^{(n_{b})})\in\mathbb{R}^{2n_{b}\times 3n_{b}}.$

Remark 1**.**

Here, we assume that the variables to be solved for are the high-level, total cooling load setpoints for individual buildings $\mathbf{u}_{b}(t)$. From this cooling load, building operators solve low-level control and optimization problems. This is customary in building control studies; see for comparison [22].

III Building-Integrated Power Network: Dynamics And Connection to OPF

In this section, we present the dynamics of the building-integrated power network and define the main variables involved in the BtG integration framework. In addition, we discuss the connection of the BtG integration model with the optimal power flow (OPF).

III-A DAE Dynamics of a Power Network with Building Loads

Let $\mathcal{B}=\{1,\ldots,n\}$ and $\mathcal{G}=\{1,\ldots,n_{g}\}$ denote the sets of buses and generators in a power network. Also, let $\mathcal{N}_{k}$ be the neighborhood set of adjacent nodes connected to the $k$th bus. Generators are indexed by $m\in\mathcal{G}$. The mechanical input power to the $m$th generator is denoted by $P_{m}$. Define generator-to-node and building-to-node incidence matrices $\mathbf{\Gamma}\in\mathbb{R}^{n\times n_{g}}$ and $\mathbf{\Pi}\in\mathbb{R}^{n\times n_{b}}$ with entries given by

[TABLE]

The transients of the $k$th bus in a power network can be modeled by the swing equation which relates the rotor angle $\delta$ with the angular velocity $\dot{\delta}$ and the angular acceleration $\ddot{\delta}$ [35]. Define $M_{k}$ and $D_{k}$ as the inertia and damping coefficients of the generator located $k$th bus; if the $k$th bus does not have a generator, then $M_{k}=0$ and $D_{k}=0$. The swing equation for the $k$th bus ($k\in\mathcal{B}$) can be written as

[TABLE]

The load at bus $k$, $P_{L_{k}}(t)$ is described by

[TABLE]

In (5), the first two terms represent uncontrollable loads, while the last one defines the controllable load. Specifically, $P_{\mathrm{BL}_{k}}(t)$ denotes the frequency-insensitive uncontrollable base load at bus $k$, which is typically available via forecasts. The term $D_{k}^{\prime}\dot{\delta}_{k}$ denotes the frequency-sensitive portion of the uncontrollable load at bus $k$. It is assumed specifically that a portion of the load at bus $k$ responds linearly to the frequency variations, which is a classical model [36]; the linear coefficient is $D_{k}^{\prime}$, and the frequency is the derivative of the angle $\dot{\delta}_{k}(t)$. It must be emphasized that the uncontrollable loads are not optimization variables. The term $\sum_{l}\pi_{k,l}P^{(l)}_{\mathrm{bldg}}(t)$ defines the load from buildings indexed by $l$ and attached to bus $k$ participating in regulation. The building load is further decomposed as $P^{(l)}_{\mathrm{bldg}}(t)=P_{\mathrm{HVAC}}^{(l)}(t)+P_{\mathrm{misc}}^{(l)}(t),$ where $l\in\{1,2,\ldots,n_{b}\}$ is the index of buildings. The quantity $P_{\mathrm{HVAC}}^{(l)}(t)$ denotes the portion of controllable power consumption of building $l$, while $P_{\mathrm{misc}}^{(l)}(t)$ represents the uncontrollable miscellaneous power consumption of building $l$ such as lighting, computers, equipment, elevators—amounting to a building’s base-load. The quantity $P_{\mathrm{HVAC}}^{(l)}(t)$ is an optimization variable, and $P_{\mathrm{misc}}^{(l)}(t)$ is typically available via forecasts.

We can now rewrite (4) as

[TABLE]

In (6), $P_{m}(t)$ for generator bus $m$ can be written as $P_{m}(t)=\bar{P}_{m}+\Delta P_{m}(t),$ where $\bar{P}_{m}$ is a solution of an optimal power flow problem—computed every 15 minutes—and $\Delta P_{m}(t)$ is the deviation from the setpoint $\bar{P}_{m}$, which will be furnished by the proposed BtG integration framework.

The angular velocity is $\dot{\delta}_{k}=\omega_{k}$, where $\omega_{k}=\omega^{\mathrm{true}}_{k}-\omega_{0}$, $\omega^{\mathrm{true}}_{k}$ is the actual frequency of the $k$th bus, and $\omega_{0}$ is the synchronous frequency, e.g., $2\pi 60\,\,\mathrm{rad/sec}$ in North America. Given (6), we obtain two first-order differential equations representing the dynamics of the $k$th bus

[TABLE]

The resulting state-space model is a nonlinear system, and is formulated as

[TABLE]

where $\mathbf{x}_{g}(t)=[\delta_{1}\;\ldots\;\delta_{n}\;\omega_{1}\;\ldots\;\omega_{n}]^{\top}=[\boldsymbol{\delta}^{\top}(t)\;\boldsymbol{\omega}^{\top}(t)]^{\top}$ is the state of the grid; $\boldsymbol{\Phi}(\boldsymbol{\delta}(t))$ is the vectorized nonlinear power flow equations in (6); $\mathbf{u}_{b}(t)=[P_{\mathrm{HVAC}}^{(1)}\;\ldots\;P_{\mathrm{HVAC}}^{(n_{b})}]^{\top}$ is the control input vector of the buildings, as defined in (2), and $\mathbf{u}_{g}(t)=\bar{\mathbf{u}}_{g}+\Delta\mathbf{u}_{g}(t)=[\bar{P}_{1}+\Delta P_{1}(t)\;\ldots\;\bar{P}_{n_{g}}+\Delta P_{n_{g}}(t)]^{\top}$ is the power network’s control variable; $\mathbf{w}_{g}(t)=[\mathbf{w}_{\mathrm{BL}}^{\top},\mathbf{w}_{\mathrm{misc}}^{\top}]^{\top}=[P_{\mathrm{BL}_{1}}\;\ldots\;P_{\mathrm{BL}_{n}},\,P_{\mathrm{misc}}^{(1)}\;\ldots\;P_{\mathrm{misc}}^{(n_{b})}]^{\top}$ is a random vector collecting the nodal base loads and miscellaneous building loads. Load forecasting is a very mature area; in the sequel, the forecast of $\mathbf{w}_{g}(t)$, denoted by $\hat{\mathbf{w}}_{g}(t)=[\hat{\mathbf{w}}_{\mathrm{BL}}^{\top},\hat{\mathbf{w}}_{\mathrm{misc}}^{\top}]^{\top}$, is assumed to be available. The state-space matrices in (7) are obtained as follows

[TABLE]

where $\mathbf{E}_{g}\in\mathbb{R}^{2n\times 2n}$ is a singular matrix, $\mathbf{A}_{g}\in\mathbb{R}^{2n\times 2n}$, $\mathbf{A}_{u_{b}}\in\mathbb{R}^{2n\times n_{b}}$, $\mathbf{B}_{u_{g}}\in\mathbb{R}^{2n\times n_{g}}$, and $\mathbf{B}_{w_{g}}\in\mathbb{R}^{2n\times(2n+n_{b})}$.

Dynamic systems of the form (7) are called differential algebraic equations (DAE), and the systems following such equations are called descriptor systems. The chief difference between descriptor systems and standard dynamical systems is that the former have a matrix $\mathbf{E}_{g}$ that multiples $\dot{\mathbf{x}}_{g}$ and may be singular. This is indeed the case here, where entries in the $\mathbf{M}$ matrix corresponding to non-generator buses are zero, giving rise to entire rows of zeros. The next section describes how the generation setpoints in (7) are computed.

III-B Connection to the Optimal Power Flow

Recall that $\mathbf{u}_{g}(t)$ is written as $\mathbf{u}_{g}(t)=\bar{\mathbf{u}}_{g}+\Delta\mathbf{u}_{g}(t),$ where $\bar{\mathbf{u}}_{g}$ contains the generator setpoints, and $\Delta\mathbf{u}_{g}(t)$ is the real-time deviation from these setpoints that automatically drives the power grid to stability after load deviations. Typically, the setpoints are computed every 5–15 minutes through solving economic dispatch or OPF routines [35]. A linearized OPF (LOPF) problem is described by the following program

[TABLE]

In (8), $J(\bar{\mathbf{u}}_{g})$ is a convex cost function that represents the generators’ cost curves. Constraint (8d) represents the safety upper and lower bounds on the generator’s active power. Vectors $\mathbf{u}_{b}$, ${\mathbf{w}}_{\mathrm{misc}}$, and ${\mathbf{w}}_{\mathrm{BL}}$ were introduced in the previous subsection and represent respectively the building HVAC loads, building miscellaneous loads, and nodal base loads. We use the notations $\hat{\mathbf{w}}_{\mathrm{misc}}$ and $\hat{\mathbf{w}}_{\mathrm{BL}}$ to emphasize that the respective forecasted versions of the building miscellaneous loads and nodal base loads enter the OPF.

Matrix $\boldsymbol{\Gamma}$ has entries defined in (3c), and thus $\mathbf{\Gamma}\bar{\mathbf{u}}_{g}$ is an $n\times 1$ vector that gives the generation for each bus of the network. Likewise, matrix $\boldsymbol{\Pi}$ has entries defined in (3f), and $\mathbf{\Pi}(\mathbf{u}_{{b}}+\hat{\mathbf{w}}_{\mathrm{misc}})$ is an $n\times 1$ vector that gives the building loads per bus. Therefore, vector $\mathbf{\Gamma}\bar{\mathbf{u}}_{g}-\mathbf{\Pi}(\mathbf{u}_{{b}}+\hat{\mathbf{w}}_{\mathrm{misc}})-\hat{\mathbf{w}}_{\mathrm{BL}}$, which appears in both (8d) and (8d), represents the net nodal power injections. With $\mathbf{1}_{n}\in\mathbb{R}^{n}$ defined as a vector of all ones, $\left(\mathbf{\Gamma}\bar{\mathbf{u}}_{g}-\mathbf{\Pi}(\mathbf{u}_{{b}}+\hat{\mathbf{w}}_{\mathrm{misc}})-\hat{\mathbf{w}}_{\mathrm{BL}}\right)^{\top}\mathbf{1}_{n}$ gives the sum of all net nodal injections, and constraint (8d) ensures the supply-demand balance. In addition, $\mathbf{F}^{\mathrm{max}}\in\mathbb{R}^{n_{l}}$ is the vector containing the thermal limits for real power flow on the $n_{l}$ branches of the network; and $\mathbf{L}_{\mathrm{ptdf}}\in\mathbb{R}^{n_{l}\times n}$ is a matrix of power transfer distribution factors [37]. This matrix maps net nodal injections to line power flows, and thus (8d) guarantees the satisfaction of line flow limits. Formulation (8) is useful in the next sections.

IV How Can Buildings Impact Power Grids?

Addressing the Main Challenges and BtG-GMPC

In the previous section, we formulate the dynamics of the buildings-integrated power network. The presence of $\mathbf{u}_{b}(t)$ in (7), exemplifies the control potential that buildings have on power system operation and control, and hence the integration advocated in this paper. In this section, we investigate the discrepancies in time-scales between the building (2) and power network dynamics (7) and discuss a formulation of the joint optimal control problem that addresses the time-scale discrepancies, while seamlessly incorporating objectives and constraints from the power grid and building clusters.

IV-A The Not-So-Cruel Curse of Time-Scales

The formulated dynamics in Sections II and III clearly operate in two different time-scales. While grid regulation problems and mechanical input power variations are often in seconds, the building dynamics and controls are much slower. For example, temperatures in buildings change slowly in comparison with frequencies and voltages in power networks.

To overcome this limitation, we design local optimal control laws that operate at different scales. Specifically, the time-step for application of building optimal control laws is $h_{b}$; and the time step for application of grid optimal control laws is $h_{g}$, where $h_{g}<<h_{b}$. This approach reflects the physical realities for these systems, and this consideration can be imposed via constraints in the optimal control problem, whose construction is the objective of this section. Since buildings possess slower dynamic behavior, we restrict the controls of buildings to be fixed for the faster time-scale of the power network.

Given this integration scheme, the discrepancy in time-scales between building and grid dynamics, the natural existence of algebraic equations in the power network model, and the necessity of including hard constraints such as tight frequency and temperature bounds, model predictive control (MPC) is the natural solution to solve the joint optimal control problem. Other control methods such as rule-based control or PID control can still be used for individual buildings, but these techniques provide inferior results in comparison to MPC as discussed in [21, 22, 23]. In addition, analytical optimal control techniques that are based on deriving a closed form solution of the optimal control law cannot be computed due to the reasons outlined above.

IV-B Discretization of the Dynamics via Gear’s Method

Another challenge with BtG integration is the presence of algebraic equations in (7) emerging from power flows of load nodes. Here, we present a simple, yet high-fidelity discretization routine for two dynamical systems with different time-scales and algebraic constraints.

First, we assume that the sampling times for the power grid [cf. (7)] and building [cf. (2)] dynamics are respectively $h_{g}$ and $h_{b}$; note that $h_{b}>>h_{g}$. The discretization we utilize in this paper is based on Gear’s method—a backward differentiation routine—for DAE (descriptor) systems [38]. The discretization of (7) can be written as follows:***For the discretization purposes, we use the linearized power flows by assuming that $\sin(\delta_{i}-\delta_{j})=\delta_{i}-\delta_{j}$. This implies that the state-space matrix $\mathbf{A}_{g}$ in (9) is different than $\mathbf{A}_{g}$ in (7), where the former includes the Laplacian matrix of the network with weights equal to line inductances. However, we simulate the system with the nonlinear power flows in Section VI.

[TABLE]

where $\bar{\mathbf{A}}_{g}=(\mathbf{E}_{g}-h_{g}\beta_{0}\mathbf{A}_{g})^{-1},\mathbf{B}_{0}=h_{g}\beta_{0}\bar{\mathbf{A}}_{g},\beta_{0}=({\sum_{i=1}^{s}1/i})^{-1}\text{,}\,$ $\alpha_{i}=(-1)^{i+1}\beta_{0}\sum_{j=i}^{s}j^{-1}\binom{j}{i}$; $k_{g}$ is the time-step for the grid dynamics; and $s$ is called the order of the method. This method requires a set of $s$ initial conditions; see Remark 2. Similarly, the discrete form of (2) can be written as follows

[TABLE]

where $\bar{\mathbf{A}}_{b}=(\mathbf{I}_{2n_{b}}-h_{b}\beta_{0}\mathbf{A}_{b})^{-1},\mathbf{B}_{1}=h_{b}\beta_{0}\bar{\mathbf{A}}_{b},\beta_{0}=({\sum_{i=1}^{s}1/i})^{-1},\alpha_{i}=(-1)^{i+1}\beta_{0}\sum_{j=i}^{s}j^{-1}\binom{j}{i}$, and $k_{b}$ is time-step for dynamic operation of buildings.

Gear’s discretization amounts to a backward Euler-like implicit method. Gear’s method is applicable to DAEs where the matrix $\mathbf{E}_{g}$ is allowed to be singular; notice that (9) does not rely on the inverse of $\mathbf{E}_{g}$. Interestingly, when the matrix $\mathbf{E}_{g}$ is identity and $s=1$, Gear’s method reduces to the standard backward Euler’s method. The principal merit of implicit methods is that they are typically more stable for solving systems with a larger step size $h$, while still performing well for systems with faster time-constants [38].

Remark 2** (Convergence of Gear’s Method).**

The states of the discretized descriptor system in (9) and (10) converge to the actual ones in a finite number of time-steps, even if the $s$-initial conditions are arbitrarily chosen [38]. A method to compute the correct initial conditions is also provided in [38].

IV-C Joint Optimal Control Problem: BtG-GMPC

The joint optimal control problem, Building-to-Grid Gear MPC (BtG-GMPC), is formulated as in (12). The variables, cost function, and constraints of BtG-GMPC are as follows:

• •

$T_{p}$ is the prediction horizon and $t$ is the initial starting point of the MPC. The formulation shows the MPC for one prediction horizon.

• •

$\mathbf{U}_{b}=\{\mathbf{u}_{b}(t+h_{b}),\mathbf{u}_{b}(t+2h_{b}),\ldots,\mathbf{u}_{b}(t+T_{p})\},\Delta\mathbf{U}_{g}=\{\Delta\mathbf{u}_{g}(t+h_{g}),\Delta\mathbf{u}_{g}(t+2h_{g}),\ldots,\Delta\mathbf{u}_{g}(t+T_{p})\},$ and $\bar{\mathbf{u}}_{g}$ are the three sets of optimization variables that we defined previously. In addition, the two sets of states defined as $\mathbf{X}_{b}=\{\mathbf{x}_{b}(t+h_{b}),\mathbf{x}_{b}(t+2h_{b}),\ldots,\mathbf{x}_{b}(t+T_{p})\}$ and $\mathbf{X}_{g}=\{\mathbf{x}_{g}(t+h_{g}),\mathbf{x}_{g}(t+2h_{g}),\ldots,\mathbf{x}_{g}(t+T_{p})\}$ are also optimization variables.

• •

The cost function $f(\Delta\mathbf{U}_{g},\bar{\mathbf{u}}_{g},\mathbf{U}_{b},\mathbf{X}_{g},\mathbf{X}_{b})$ is defined as the weighted summation of the building costs, the steady-state LOPF costs, the penalties on the deviation from the steady-state generation, and the deviation cost from the nominal synchronous frequency:

[TABLE]

where

• –

$J(\bar{\mathbf{u}}_{g})$ is the LOPF cost function (8a). The parameters of this cost function are widely available in the power systems literature [37].

• –

$\mathbf{c}_{b}(t+k_{b}h_{b})$ is a time-varying vector representing the cost of electricity at time $t+k_{b}h_{b}$. These prices are the wholesale price of electricity for commercial building operators.

• –

The third term in $f(\cdot)$ penalizes the deviations in the mechanical power setpoints of generators through a quadratic cost function, with matrix $\mathbf{R}\in\mathbb{R}^{n_{g}\times n_{g}}$ being the quadratic penalty matrix, which is assumed to be positive semidefinite.

• –

The fourth term in $f(\cdot)$ penalizes the deviations of the generator frequencies from their nominal value using matrix $\mathbf{Q}\in\mathbb{R}^{2n\times 2n}$. The reader is referred to [39, 35] for related constructions. This cost function is similar to the linear quadratic regulator, which is used in dynamical systems and power network stability studies.

• –

The terms multiplying the summations are meant to average the building and grid costs across the planning horizon $T_{p}$.

• •

Constraints (12b)–(12d) depict the dynamics of the building-integrated power grid, as well as lower and upper bounds on the states and inputs of the grid states and controls. Note that ${\mathbf{x}}_{g}(t+k_{g}h_{g})=\mathbf{f}_{g}(\mathbf{x}_{g},\mathbf{u}_{g},\mathbf{u}_{b},\hat{\mathbf{w}}_{g}\,|\,t,s)\overset{\Delta}{=}\bar{\mathbf{A}}_{g}\sum_{i=1}^{s}\alpha_{i}\mathbf{E}_{g}{\mathbf{x}}_{g}(t+h_{g}(k_{g}-i))+\mathbf{B}_{0}({\mathbf{A}}_{\mathbf{u}_{b}}\mathbf{u}_{b}(t+k_{g}h_{g})+{\mathbf{B}}_{u_{g}}\mathbf{u}_{g}(t+k_{g}h_{g})+\mathbf{B}_{w_{g}}\hat{\mathbf{w}}_{g}(t+k_{g}h_{g}))$, where $s$ corresponds to the order of Gear’s method.

• •

Constraints (12e)–(12g) represent the building cluster dynamics and the bounds on the states and inputs of the individual buildings, while constraint (12h) imposes the constraints of the LOPF as discussed in the previous section.

• •

The final constraint (12i) represents the idea of the time-scales integration whereby the building control variables are kept constant between two consecutive building instances. Since $h_{b}>h_{g}$, we assume that between two consecutive building sampling instances (i.e., $k_{b}h_{b}$ and $(k_{b}+1)h_{b}$), the building controls $u_{b}(k_{b}h_{b})$ are all constant variables to be found. Hence, for all $\forall\,k_{g}h_{g}\in[k_{b}h_{b},(k_{b}+1)h_{b})$, $\mathbf{u}_{b}(k_{b}h_{b})=\mathbf{u}_{b}(k_{g}h_{g})=\bar{\mathbf{u}}_{b}.$

BtG-GMPC:

$\displaystyle\operatorname*{minimize}_{\begin{subarray}{c}\mathbf{U}_{b},\Delta\mathbf{U}_{g},\bar{\mathbf{u}}_{g}\\ \mathbf{X}_{b},\mathbf{X}_{g}\end{subarray}}\,\,\,$ $\displaystyle f(\Delta\mathbf{U}_{g},\bar{\mathbf{u}}_{g},\mathbf{U}_{b},\mathbf{X}_{g},\mathbf{X}_{b})$

(12a)

$\displaystyle\mathrm{subject\,to}\hskip 8.5359pt\;$ $\displaystyle\,{\mathbf{x}}_{g}(t+k_{g}h_{g})=\mathbf{f}_{g}(\mathbf{x}_{g},\mathbf{u}_{g},\mathbf{u}_{b},\hat{\mathbf{w}}_{g}\,|\,t,s)$

(12b)

$\displaystyle\Delta\mathbf{u}_{g}^{\mathrm{min}}\leq\Delta\mathbf{u}_{g}(t+k_{g}h_{g})\leq\Delta\mathbf{u}_{g}^{\mathrm{max}}$

(12c)

$\displaystyle\mathbf{x}_{g}^{\mathrm{min}}\leq\mathbf{x}_{g}(t+k_{g}h_{g})\leq\mathbf{x}_{g}^{\mathrm{max}}$

(12d)

$\displaystyle\forall\;k_{g}\in\{1,\ldots,T_{p}/h_{g}\}$

$\displaystyle\mathbf{x}_{b}(t+k_{b}h_{b})=\mathbf{f}_{b}(\mathbf{x}_{b},\mathbf{u}_{b},\hat{\mathbf{w}}_{b}\,|\,t,s)$

(12e)

$\displaystyle\mathbf{u}_{b}^{\mathrm{min}}\leq\mathbf{u}_{b}(t+k_{b}h_{b})\leq\mathbf{u}_{b}^{\mathrm{max}}$

(12f)

$\displaystyle\mathbf{x}_{b}^{\mathrm{min}}\leq\mathbf{x}_{b}(t+k_{b}h_{b})\leq\mathbf{x}_{b}^{\mathrm{max}}$

(12g)

$\displaystyle\forall\;k_{b}\in\{1,\ldots,T_{p}/h_{b}\}$

$\displaystyle\eqref{equ:DCOPF1},\eqref{equ:DCOPF3}\text{}$

(12h)

$\displaystyle\mathbf{u}_{b}(t+k_{g}h_{g})=\bar{\mathbf{u}}_{b}=\mathbf{u}_{b}(t+k_{b}h_{b})$

(12i)

$\displaystyle\forall\,k_{g}h_{g}\in[k_{b}h_{b},(k_{b}+1)h_{b}).$

Algorithm 1 illustrates a routine that implements BtG-GMPC’s rolling horizon window along with the integration of the LOPF problem. Given the BtG-GMPC parameters (including the first $s$-initial steps of the discretized dynamics), the algorithm computes the optimal solutions to the LOPF problem and the joint MPC. For simplicity, we assume that the prediction horizon $T_{p}$ is equivalent to the time-scale in which the optimal dispatch is solved, i.e., 5 to 15 minutes.

We also assume that $h_{g}<h_{b}<T_{p}<<T_{\mathrm{final}}$ and $h_{b}/h_{g}$, $T_{p}/h_{b}$, $T_{p}/h_{g}$ are all positive integers. The algorithm starts by finding the solution to the generator’s operating points $\bar{\mathbf{u}}_{g}$ for any multiple of the prediction horizon $T_{p}$, as well as the deviation from this setpoint $\Delta\mathbf{u}_{g}(t)$ and $\mathbf{u}_{b}(t)$ up until the next planning horizon, and so on. As in classical MPC routines, only the first instance of the optimal control trajectory is applied, while the rest are discarded. Note that the BtG-GMPC with LOPF is only solved for when $t$ (the counter) is a multiple of $T_{p}$. If $t$ is not a multiple of $T_{p}$, but a multiple of the building’s sampling time $h_{b}$, the building and grid controls are computed. The final case captures the gap between the two time-scales: where the building and grid controls are applied, the building controls are kept constant from the previous optimal computations, while grid controls are computed in the meantime for every grid sampling time.

Remark 3** (Tractability of BtG-GMPC).**

Problem (12) is a quadratic program. Even for large-scale systems, this optimization routine is tractable, and can be solved by off-the-shelf solvers such as CPLEX, MOSEK, or Matlab’s QuadProg.

Remark 4** (Fast MPC and Time-Complexity).**

The BtG-GMPC optimization is applied online as predictions for the uncontrollable inputs might not be available for times greater than the prediction horizon $T_{p}$. However, given that prediction for uncontrollable inputs are available prior to the start of the day, this problem can be solved offline. If solved online, fast online MPC algorithms for quadratic programs have been developed in [40] and can be immediately applied to BtG-GMPC. Otherwise, the problem can be solved offline, which eases the communication requirement of exchanging optimal solutions. Note that in BtG-GMPC, the maximum total number of variables at each time-step is equal to $3n_{b}+2n+2n_{g}=N$. As reported in [40], MPC formulations take $\mathcal{O}(T_{p}\cdot N^{3})$ at each time-step. This is based on novel interior-point-based implementations.

V Comparisons with Decoupled BtG Designs

The BtG integration framework developed in this paper enables building and grid operators to jointly optimize their decisions. This section analytically formalizes and compares the BtG framework with decoupled designs, in which the grid and the building operators schedule generation and the building power consumption separately.

Problem (12) jointly optimizes over two groups of variables: grid decisions $(\bar{\mathbf{u}}_{g},\Delta\mathbf{U}_{g},\mathbf{X}_{g})$ and building decisions $(\mathbf{U}_{b},\mathbf{X}_{b})$. The problem can be written in the following way, which brings out the coupling between the grid and building decisions.

[TABLE]

The objective (13a) corresponds to (12a). All constraints of problem (12) are captured by one of the constraints of problem (13). Specifically, set $\mathcal{C}$ represents the coupling between grid and building decisions, that is, $\mathcal{C}$ is the set of all decisions $(\bar{\mathbf{u}}_{g},\Delta\mathbf{U}_{g},\mathbf{X}_{g},\mathbf{U}_{b},\mathbf{X}_{b})$ that satisfy (12b) and (8d). Set $\mathcal{F}_{g}$ represents all the constraints pertaining to grid decisions only, that is, constraints (8d), (12c), and (12d). Set $\mathcal{F}_{b}$ represents all constraints pertaining to building decisions only, that is, (12e), (12f), (12g), and (12i).

Grid and building controls are optimized separately in traditional power systems. In particular, building operators control the building HVAC load $\mathbf{U}_{b}$ based on electricity prices $c_{b}$, forecasted weather conditions, and occupancy behaviors [with the latter two captured by $\mathbf{w}_{b}(t)$]. Bang-bang control or more sophistacted MPC methods may be used to determine the building HVAC load. In turn, grid operators forecast the grid load, which comprises the building load $\mathbf{u}_{b}+\hat{\mathbf{w}}_{\mathrm{misc}}$ and the remaining base load $\hat{\mathbf{w}}_{\mathrm{BL}}$, and determine generator setpoints and mechanical power adjustments.

The previously described process can be formalized within the proposed framework as follows. Supposing the building HVAC controls are optimized via MPC, the building operators solve the following optimization problem:

[TABLE]

That is, the objective is to minimize the cost of building operation, subject to the dynamical constraints of the buildings as well as the state and control bounds. It is actually not difficult to see that the previous optimization can be performed by each building operator separately. Let $\mathbf{U}_{b}^{\mathrm{MPC}},\mathbf{X}_{b}^{\mathrm{MPC}}$ be the solution of problem (14).

The grid operator optimizes the generator setpoints and mechanical power adjustments based on the predicted grid load, which includes building loads $\mathbf{U}_{b}^{\mathrm{MPC}}$. The grid operator thus solves the following optimization problem:

[TABLE]

Let $\bar{\mathbf{u}}_{g}^{\mathrm{MPC}},\Delta\mathbf{U}_{g}^{\mathrm{MPC}},\mathbf{X}_{g}^{\mathrm{MPC}}$ be the solution of (15).

The total cost of operation for the previously mentioned decoupled design is $f_{g}^{\mathrm{MPC}}+f_{b}^{\mathrm{MPC}}$, where $f_{g}^{\mathrm{MPC}}$ and $f_{b}^{\mathrm{MPC}}$ are respectively the optimal values of (15) and (14). The relationship with the cost from BtG integration $f^{*}_{\mathrm{BtG}}$ [cf. (13)] is provided in the following proposition.

Proposition 1**.**

It holds that

[TABLE]

Proof:

Consider the decisions $(\bar{\mathbf{u}}_{g}^{\mathrm{MPC}},\Delta\mathbf{U}_{g}^{\mathrm{MPC}},\mathbf{X}_{g}^{\mathrm{MPC}},\mathbf{U}_{b}^{\mathrm{MPC}},\mathbf{X}_{b}^{\mathrm{MPC}})$. These are feasible for problem (13), because they are feasible for problems (15) and (14). Since $f^{*}_{\mathrm{BtG}}$ is the optimal value of (13), it holds for any feasible point that

[TABLE]

But $\mathbf{U}_{b}^{\mathrm{MPC}},\mathbf{X}_{b}^{\mathrm{MPC}}$ is the solution of (14), and therefore, $f_{b}(\mathbf{U}_{b}^{\mathrm{MPC}},\mathbf{X}_{b}^{\mathrm{MPC}})=f_{b}^{\mathrm{MPC}}$ holds. Likewise, $\bar{\mathbf{u}}_{g}^{\mathrm{MPC}},\Delta\mathbf{U}_{g}^{\mathrm{MPC}},\mathbf{X}_{g}^{\mathrm{MPC}}$ is the solution of (15), and $f_{g}(\bar{\mathbf{u}}_{g}^{\mathrm{MPC}},\Delta\mathbf{U}_{g}^{\mathrm{MPC}},\mathbf{X}_{g}^{\mathrm{MPC}})=f_{g}^{\mathrm{MPC}}$ holds. Utilizing the latter two optimal values in (17), (16) follows. ∎

The previous proposition asserts that the decoupled design incurs a total system cost that is no smaller than the one of the proposed BtG integration scheme. Intuitively, the BtG design allows to jointly look for grid and building control actions in set $\mathcal{C}$, as opposed to fixing the building controls first, and then solving the grid optimization problem.

Attention is now turned to the case where bang-bang control is used to determine the HVAC loads. Bang-bang control is the simplest type and most common type of HVAC control where the controller follows a strict temperature set point (e.g., $22.22\degree$C). The HVAC control system is switched on (or off) as soon as the zone temperature exceeds (or is below) the dead band which is generally $\pm 0.5\degree\mathrm{C}$. Bang-bang control thus does not optimally solve (14), because it restricts when the control system is turned on. To make a fair comparison, it is supposed that the dead band is tuned so that the resulting temperatures do not go outside of the intervals specified by (12g). Likewise, the control action is not allowed to exceed the bounds dictated by (12f). By design, bang-bang control adheres to the building dynamics described by (12e).

Let $\mathbf{U}_{b}^{\mathrm{BB}},\mathbf{X}_{b}^{\mathrm{BB}}$ be the control actions and resulting states of bang-bang control; the previous discussion implies that $(\mathbf{U}_{b}^{\mathrm{BB}},\mathbf{X}_{b}^{\mathrm{BB}})$ is feasible for problem (14). Let $f^{\mathrm{BB}}_{b}$ be the resulting cost of building operation. Also, let $\check{f}_{g}^{\mathrm{MPC}}$ be the optimal value of (15) where $\mathbf{U}_{b}^{\mathrm{MPC}}$ is replaced by $\mathbf{U}_{b}^{\mathrm{BB}}$; the resulting system cost is $\check{f}_{g}^{\mathrm{MPC}}+f^{\mathrm{BB}}_{b}$. The following proposition relates the costs derived from bang-bang control with the costs of MPC based operation.

Proposition 2**.**

Suppose that $(\mathbf{U}_{b}^{\mathrm{BB}},\mathbf{X}_{b}^{\mathrm{BB}})$ is feasible for problem (14). Then it holds for the resulting building operation cost that

[TABLE]

and for the system cost that

[TABLE]

Proof:

Eq. (18) follows from the fact that $(\mathbf{U}_{b}^{\mathrm{BB}},\mathbf{X}_{b}^{\mathrm{BB}})$ is feasible for (14), while $f^{\mathrm{MPC}}_{b}$ is the optimal value of the same problem.

To prove (19), let $\check{\bar{\mathbf{u}}}_{g}^{\mathrm{MPC}},\Delta\check{\mathbf{U}}_{g}^{\mathrm{MPC}},\check{\mathbf{X}}_{g}$ be the solution of (15) where $\mathbf{U}_{b}^{\mathrm{MPC}}$ is replaced by $\mathbf{U}_{b}^{\mathrm{BB}}$. It follows that $(\check{\bar{\mathbf{u}}}_{g}^{\mathrm{MPC}},\Delta\check{\mathbf{U}}_{g}^{\mathrm{MPC}},\check{\mathbf{X}}_{g}^{\mathrm{MPC}},\mathbf{U}_{b}^{\mathrm{BB}},\mathbf{X}_{b}^{\mathrm{BB}})$ is feasible for problem (13). Thus it holds that

[TABLE]

from which (19) follows. ∎

The previous proposition asserts that MPC for building HVAC controls incurs smaller costs than bang-bang control, as long as the bang-bang control adheres to the same constraints as MPC—a fact that has previously been demonstrated in the building literature. But more importantly, similarly to Proposition 1, it is concluded that fixing the building controls to the particular scheme cannot improve the system costs over jointly designing the building and grid controls.

For easier reference, we refer to the decoupled design where the building HVAC loads are computed via bang-bang control and subsequently the grid is optimized via MPC as Scenario I (with optimal value $\check{f}_{g}^{\mathrm{MPC}}+f^{\mathrm{BB}}_{b}$, cf. Proposition 2). The respective design where the building HVAC loads are optimized via MPC is referred to as Scenario II (with optimal value $f_{g}^{\mathrm{MPC}}+f_{b}^{\mathrm{MPC}}$, cf. Proposition 1). The developed BtG framework is referred as Scenario III (with optimal value $f^{*}_{\mathrm{BtG}}$) . The next section provides numerical simulations that test the previously mentioned designs, and corroborate the analytically derived comparisons.

VI Case Studies

In this section, we investigate the impact of the proposed BtG integration on the performance of grid’s stability and the cost-effectiveness of building control systems.

VI-A Experimental Setup and Parameters

The case studies are performed on various power networks for $T_{\mathrm{final}}=24$ hours. In particular, we use casefiles (case9,case14,case30,case57) from Matpower [37] to test different power networks, parameters, and total number of buildings. Table I documents the total number of buses, generators, buildings in the aforementioned casefiles. All of the data and the codes needed to reproduce the results are made available on this page: https://github.com/ahmadtaha1/BtG. The codes are simulated using CPLEX’s quadratic program solver [41], and are written in the Matlab environment. The simulations are performed on a PC running Windows 10 Enterprise, Intel(R) Xeon(R) CPU E3-1271 V3 with a 3.60-GHz processor, and 32 GB of RAM. The code allows the testing of any power network with custom-defined total number of buildings at predefined load buses.

The parameters, exact constraints, weather data, electricity prices, and other details of the problem are carefully chosen to reflect reasonable conditions and IEEE/ASHRAE practices. The parameters are chosen as follows.

• •

The grid’s base load forecast and the miscellaneous loads of all buildings are chosen carefully. The base-load variations are based on the National Grid New York electricity company’s posted demand curve from the Standard Service in New York: https://www9.nationalgridus.com/niagaramohawk/business/rates/5_load_profile.asp. In addition, the disturbances to the buildings, hourly data for $T_{\mathrm{amb}},\dot{Q}_{\mathrm{sol}},\dot{Q}_{\mathrm{int}},P_{\mathrm{BL}},$ and $P_{\mathrm{misc}}$ are all included in the Github link.

• •

Building loads are modeled based on one reference commercial building located at the main campus of the University of Texas at San Antonio. For this reference building, the construction materials are known and further determined by ASHRAE standard 90.1-2016. The building size is calculated based on design documents.

• •

The RC-parameters for all buildings are obtained using a normal distribution around the following reference (mean) building. The mean parameters are $R_{1}=R_{2}=1.16\times 10^{-4}$ ($\mathrm{\degree C/W}$), $R_{\mathrm{win}}=6.55\times 10^{-3}$ ($\mathrm{\degree C/W}$), $C_{\mathrm{zone}}=7.033\times 10^{9}$ ($\mathrm{J/\degree C}$), and $C=1.133\times 10^{9}$ ($\mathrm{J/\degree C}$). The average building size is around $10,000\,\mathrm{m}^{2}$.

• •

The prices of electricity for the HVAC loads are reproduced from [22]. The cost functions for generator mechanical power setpoints $\bar{u}_{g}$ are extracted from Matpower [37]. The same quadratic cost is used for the variations ($\Delta u_{g}(t)$). The deviation in frequency is penalized with $Q_{k}=50000$ \/(\mathrm{rad/sec})^{2}$($Q_{k}$is the$k$th diagonal entry in$\mathbf{Q}$), and the angles are left without any penalties in the$\mathbf{Q}$-matrix.

• •

The parameters of the power network, including line parameters $b_{kj},M_{k}$, and $D_{k}$ for all buses $k\in\mathcal{B}$, are obtained from Matpower [37] and the power system toolbox [42].

• •

We choose $h_{g}=10\sec$, $h_{b}=300\sec$, and a prediction horizon $T_{p}=900\sec$. For simplicity, 1st order Gear’s method is used in the simulations. The value of $h_{g}$ is consistent with the discrete time interval at which automatic generation control commands are dispatched [43, Sec. 12.3].

• •

The bound-constraints in (12) are as follows: (a) $59\leq f_{k}=\omega_{k}^{\mathrm{true}}/2\pi\leq 61$ (Hz), (b) $21.5\leq T_{\mathrm{zone}}\leq 23$ (°C) for time-periods between 8AM and 8PM, (c) $22\leq T_{\mathrm{zone}}\leq 25$ (°C) for time-periods between 8PM and 8AM, and (d) $0\leq P_{\mathrm{HVAC}_{l}}\leq 800$ (KW). The limits on the output power of generators can be found in Matpower [37].

VI-B *Impact on Frequency Regulation & Energy Savings *

In this section, we present the numerical results for the BtG-GMPC and Algorithm 1 (also named Scenario III in the sequel), in comparison with solving the optimal control of buildings and power grids separately via MPC (Scenario II). In addition, we compare the results of BtG-GMPC with bang-bang control of HVAC systems, which is still very common in today’s industries, combined with grid-only MPC (Scenario I). The three scenarios are analytically discussed in Section V, have increasing sophistication, but use the same parameters, initial conditions, constraints, and costs. The reader is referred to Section V for a comparison between the decoupled designs (Scenarios I and II) and BtG-GMPC.

For brevity, we only show the plots for case57 from [42] with 1822 commercial buildings, but present a cost comparison for all other casefiles. The uploaded Github codes contain the data for other simulations with the corresponding figures for building and power network states and optimal control inputs. In addition, Table II shows the cost comparison between the different scenarios. The cost functions are defined as follows: $\mathbf{x}_{g}^{\top}\mathbf{Q}\mathbf{x}_{g}$ denotes the frequency deviation cost; $\Delta\mathbf{u}_{g}^{\top}\mathbf{R}\Delta\mathbf{u}_{g}$ represents the mechanical input power deviations cost; $J(\bar{\mathbf{u}}_{g})$ is the LOPF cost; $\mathbf{c}_{b}^{\top}\mathbf{u}_{b}$ depicts the HVAC cooling load costs. Note that these costs are all multiplied by $1,000. In the next section, we compute the perfect cost reduction and compare the different scenarios. The percent reduction in cost is computed as follows:

[TABLE]

where the ‘Cost’ refers to any of the reported cost functions in Table II and X,Y correspond to any of the three scenarios.

VI-B1 Scenario I

The resulting HVAC power consumptions and corresponding zone temperatures for all buildings are shown in Fig. 2-(c,d). As expected, the bang-bang building control maintains the temperature in the aforementioned band.

After simulating this case over a period of 24 hours, these bang-bang HVAC inputs are provided to the grid optimal control problem, which is formulated as an MPC based on Gear’s method. Figs. 2-(c,d) show the frequencies, zone temperatures, total generation and HVAC loads for all 1822 buildings and 57 buses in the network. Due to the intermittent nature of the HVAC load, the grid frequency experiences significant deviations from its nominal value (60 Hz); see Fig. 2-(c). The frequency variations are more prominent at durations when the total load is at peak values. The total costs for Scenario I are provided in Table II.

VI-B2 Scenario II

In this scenario, we solve the MPC problems for buildings and the power grid separately. First, the building optimal controls are computed via the same MPC formulated in (12), while eliminating the power grid constraints and variables. The MPC solution for the building’s HVAC loads is then fed into a grid-only MPC. This scenario is useful in the sense that grid operators can model the building’s load via a classical building MPC model—this can be viewed as a decoupled BtG-GMPC. The numerical results of this case show an improvement in the HVAC power consumption of buildings from Scenario I. Specifically, the cooling load costs decreased (from Scenario I to Scenario II) by 42.9, 16.6, 16.4, and 16.5% for Case 9, Case 14, Case 30, and Case 57. In addition, the overall system costs decreased by an average of 29.7% for all casefiles (from Scenario I to II).

VI-B3 Scenario III (BtG-GMPC)

In this scenario, we test the performance of Algorithm 1 and the corresponding optimization problem (12). Scenario I and II are separately compared to Scenario III.

As asserted by the analytical discussion in Section V, the BtG-GMPC results show significant improvement in grid’s frequency deviations and the overall costs as shown in Fig. 2-(a) and Table II. The frequency variations are notably lower in comparison with the previous two scenarios, and this is clear from a comparison between Figs. 2-(a) and (c). Also, there is a significant reduction in the overall cost of operation. The results show around 43, 35, 30, and 36% total cost reduction between Scenarios I and III (for Case 9, Case 14, Case 30, and Case 57), and 20, 15, 9, and 2% total cost reduction between Scenarios II and III. Also, the grid frequency deviations are reduced, leading to a decrease in the cost of grid operation by an average of 74.53% from Scenario I to III and 46.75% from Scenario II to III for all casefiles. Finally, the plots in Fig. 2 show that the zone temperatures are well-maintained within the required range, with very little fluctuation (in comparison with Scenario I), and the total HVAC consumption is reduced by an average of 16.52% which is also clear from Figs. 2-(b,d). For brevity, we do not show the plots for Scenario II and other scenarios for other power networks as they are included in the Github link.

VI-B4 Computational Speed

The computational times for running a 24-hour time-horizon are shown in Table III for the four power networks. The power network Case57 requires 3.74 hours to run for the entire 24-hour simulation horizon. The computational times project that if a larger power network or a micro-grid is considered, with potentially 10,000 buildings or more (and thousands of buses), Algorithm 1 can still be implemented in real time. This is due to two reasons. Firstly, a system operator or a utility company (see Remark 5) will have more computational power at their disposal. Secondly, a fast MPC routine as the one described in Remark 4 can also be implemented to speed up the computations.

VI-B5 Robustness to Forecast and Model Uncertainty: Offline, Day-Ahead Solutions

In the previous section, Fig. 2 and Table II present the results from the MPC solution ($\mathbf{X}_{g},\mathbf{X}_{b},\mathbf{U}_{b},\mathbf{U}_{g}$) using the predicted unknown inputs ($\hat{\mathbf{w}}_{g},\hat{\mathbf{w}}_{b}$) for the linearized power system using Gear’s method. In this section, we 1) extract the optimal MPC control variables $(\mathbf{U}_{b},\mathbf{U}_{g})$ for the entire simulation horizon; 2) feed these inputs to the nonlinear continuous-time DAE solver for (7) and linear ODE solver for (1); 3) add zero-mean Gaussian noise with 10% standard deviation from the unknown inputs $\hat{\mathbf{w}}_{g}$ in (7) and $\hat{\mathbf{w}}_{b}$ in (1); and 4) perturb the 3R-2C building model (1) with zero-mean Gaussian noise with 10% standard deviation from the nominal matrices $\mathbf{A}_{b},\mathbf{B}_{u_{b}}$, and $\mathbf{B}_{w_{b}}$.

Four major reasons justify this numerical simulation, namely, 1) to validate Gear’s discretization method; 2) to assess the performance of the integration framework on the nonlinear continuous-time DAE model for the grid under mismatch between the forecasted and true disturbances; 3) to test the temperature behavior of buildings under model uncertainties that could be a result of parametric misidentification from building operators; and 4) to examine a scenario where the demand response signals are communicated a day prior to the schedules. The latter essentially alleviates the real-time communication burden of the BtG integration framework, by allowing the grid operators to send the demand response schedules way ahead in time and thereby avoiding the necessity to communicate the schedules in real time.

Fig. 3 illustrates that the frequency deviations and zone temperatures are kept within reasonable ranges, even under significant parametric and load uncertainty (10% load mismatch is relatively large in power networks), and the results hence depict that BtG-GMPC is robust to significant disturbances. Note that the MPC scheme is solved offline given the prediction of the loads and the temperatures, and the resulting controls were used as inputs to the nonlinear DAE solver, which demonstrated the good performance of the developed framework. Specifically, the zone temperatures for all buildings are still within the acceptable range, although some buildings experience zone temperatures of $25^{\circ}$C and $20^{\circ}$C, which is due to the parametric mismatch. This can be compared to Fig. 2-b. In addition, the load mismatch does not destabilize the grid’s frequencies; see Fig. 3 and Fig. 2-a for comparison.

Remark 5** (Who Solves BtG-GMPC?).**

BtG-GMPC assumes the knowledge of various parameters such as building RC-constants and generator cost curves. We consider the following: 1) A system operator or a large utility ideally solves BtG-GMPC; 2) commercial building operators contributing to this routine are required to provide modeling parameters for their buildings; and 3) the global signals computed are communicated to the now-contributing operators of individual buildings and generators. The added value of this coupling is two-fold. First, the theoretical impact of expanding the feasible space of two separate problems ensures that the coupled problem’s solution will be superior to the decoupled one as illustrated in the previous section. Second, the coupling translates into tangible impact for buildings and the grid, as has been demonstrated in this section. With that in mind, the communication costs are not considered here, and it is assumed that the computed optimal setpoints are communicated instantly to individual buildings and generators.

VII Paper Summary and Limitations; Future Work

This paper introduces the first explicit building-to-grid integration dynamical model with optimal power management formulations and different time-scales. The paper considers realistic, high-level building models and frequency-focused grid dynamics, in addition to algebraic equations modeling the nodes without generation. We also introduced Gear’s method as a high fidelity DAE discretization routine that is leveraged to model BtG integration. The developed framework and optimization problem BtG-GMPC provides setpoints for individual buildings and power grid generators, as well as buildings-aware, optimized optimal power flow setpoints for generators. The formulated problem can be solved efficiently using any quadratic program solver. Case studies have demonstrated the impact of BtG-GMPC on reducing overall energy costs and minimizing frequency deviations.

We have kept the dynamical models simple as the focus is on energy consumption and frequency deviations. However, the framework is general and interested researchers can seamlessly extend BtG-GMPC to include advanced models. Given that, three main challenges are not addressed here.

• •

We do not consider the problem of controlling or adjusting the reactive power of buildings, and its impacts on regulating the grid’s voltages. Since the focus of this work is on the framework with specific impacts on frequency regulation, we leave this natural extension to future work. Note that Gear’s method and the discretization still holds for models with reactive powers and voltages.

• •

The BtG framework developed in this paper pertains to normal grid operation. Unplanned incidents and contingencies can occur and require appropriate response. In order to deal with plausible contingencies, system operators include reserve scheduling in OPF [44]. Incorporating reserve scheduling is thus an interesting future direction.

• •

The impact of slow communications between grid and building operators, load prediction errors, and mismatch in building parameters are all investigated in Section VI, and shown to have little impact on the system states. A more sophisticated BtG-GMPC (12) that incorporates uncertainty in loads and building models is an important improvement that yields a BtG routine tolerant to these unknown inputs.