# Identifying Weakly Connected Subsystems in Building Energy Model for
Effective Load Estimation in Presence of Parametric Uncertainty

Arpan Mukherjee Corresponding author: Email<EMAIL_ADDRESS>Anna Kuechle
Szweda  Andrew Alegria  Rahul Rai  Tarunraj Singh

(Department of Mechanical and Aerospace Engineering, University at Buffalo-
SUNY)

###### Abstract

It is necessary to estimate the expected energy usage of a building to
determine how to reduce energy usage. The expected energy usage of a building
can be reliably simulated using a Building Energy Model (BEM). Many of the
numerous input parameters in a BEM are uncertain. To ensure that the building
simulation is sufficiently accurate, and to better understand the impact of
imprecisions in the input parameters and calculation methods, it is desirable
to quantify uncertainty in the BEM throughout the modeling process.
Uncertainty quantification (UQ) typically requires a large number of
simulations to produce meaningful data, which, due to the vast number of input
parameters and the dynamic nature of building simulation, is computationally
expensive. Uncertainty Quantification (UQ) in BEM domain is thus intractable
due to the size of the problem and parameters involved and hence it needs an
advanced methodology for analysis. The current paper outlines a novel Weakly-
Connected-Systems (WCSs) identification-based UQ framework developed to
propagate the quantifiable uncertainty in the BEM. The overall approach is
demonstrated on the physics-based thermal model of an actual building in
Central New York.

keywords: Building Energy Simulation, Uncertainty Quantification, Thermal Load
Estimation, Weakly Connected Subsystem

##  1 Introduction

In 2016, the building sector used approximately 40% of the energy produced in
the United States [1], and in 2009, buildings contributed 2,184.6 million
metric tons of carbon dioxide equivalent of greenhouse gases due to emissions
[2]. As a result, governments and municipalities have committed to reducing
the energy consumption of buildings. New York State, for example, enacted
Executive Order 88, which requires government buildings to reduce energy
consumption by 20% [3], and the Public Service Commission in New York has
established a public fund to target energy efficiency measures [4].

It is necessary to estimate the expected energy usage of a building to
determine how to reduce energy usage. The expected energy usage of a building
can be reliably simulated using a Building Energy Model (BEM). Currently,
there are several popular options for software to estimate the dynamic energy
use, such as EnergyPlus, eQUEST, TRACE 700, and Carrier HAP [5, 6, 7, 8], each
of which utilizes complex numerical methods to calculate building loads. The
BEMs incorporate energy load estimation in the presence of HVAC systems,
lights, and receptacle loads, and requires specification of a substantial
number of input parameters.

Many of the numerous input parameters in a BEM are uncertain. Some of these
parameters are assumed in the model, while some of these parameters are
measured. Additionally, the owner, designer, or modeling professional must
estimate others. In each case, uncertainty is brought into the modeling
results. These uncertainties in the BEM can be categorized into two types:
uncertainties in the measured parameters and the uncertainty in modeling [9,
10]. To ensure that the building simulation is sufficiently accurate, and to
better understand the impact of imprecisions in the input parameters and
calculation methods, it is desirable to quantify uncertainty in the BEM
throughout the modeling process. Uncertainty quantification (UQ) typically
requires a large number of simulations to produce meaningful data, which, due
to the vast number of input parameters and the dynamic nature of building
simulation, is computationally expensive [11, 12]. We must, therefore,
simplify the BEM to a less computationally expensive model. To enable
computationally efficient UQ of BEM, the following main contributions are made
in this paper:

  * •

Adaptation of lumped resistor-capacitance network model to generate reduced
order differential equation for BEM simulation

  * •

Usage of a gray-box method in conjunction with black-box Kalman filter to
enable BEM parameter estimation

  * •

Application of a novel WCSs identification-based UQ framework developed [13]
to propagate the quantifiable uncertainty in the BEM.

Next, all the three main contributions are further elaborated. The well-
established technique of a lumped resistor-capacitor (RC) network is used to
model the reduced-order zone heat balance equations. The RC method assumes
that the components of the zone load can be estimated by a discrete number of
resistances and capacitances, and the system is treated as equivalent to an
electrical circuit [14]. In particular, the heat transfer due to the zone
envelope can be reduced to a three resistance, two capacitance (3R2C) thermal
network. With the RC network, one can reduce the load calculations to first-
order differential equations that provide a suitable framework to carry out UQ
in BEMs.

The BEM simulation can be carried out using one of the three main models [15,
16]:

  1. 1.

White box models: In white box models detailed information about the known
physical process is used to predict the future states.

  2. 2.

Black-box models: In black-box models measured data is used to estimate non-
physical parameters that abstractly represent the BEM performance.

  3. 3.

Gray-box models: Gray-box models combine both the white-box and black-box
methods through the use of a reduced order model with known building
properties for modeling the physical process. Such models can also be easily
combined with Kalman filtering to enable parameter estimation. Usage of gray-
box models reduces computation time while improving the accuracy of
predictions.

Because of a large number of input parameters, despite the reduced order
model, the resulting filtering problem is computationally expensive. The novel
method of identifying Weakly Connected Subsystems (WCSs) [17, 18, 13] is used
to address this problem of performing UQ in the BEM. WCS-based UQ approach
allows us to group coupled parameters to minimize the associated computation
time required while continuing to propagate the quantifiable uncertainty
present in the BEM. The optimal estimation problem to estimate parameters for
UQ of BEM is the heart of this research.

A representative case study, a well-documented building, located in Central
New York, is used for modeling purposes. A separately simulated model of the
represented building provides input data for the reduced order model, which we
treat as actually measured information, with some initial uncertainty.

The paper is organized as follows: Section 2 describes the related works in
two distinct but related domains: (1) the lumped capacitance RC network and
(2) uncertainty quantification. Section 3 explains the methodology followed to
formulate the state space matrix and the associated input parameters for the
model. Section 3 also outlines the Weakly Connected Subsystems based solution
methodology. Section 4 provides details related to the case study building
model. Section 5 shows the results of the simulation, the cluster analysis,
estimated temperatures, and non-envelope loads. Section 6 discusses the main
conclusions and potential avenues for future work.

##  2 Related Works

###  2.1 Lumped capacitance RC network

The lumped capacitance RC network reduced-order model method has been reliably
used in a diverse set of research work to estimate the heat transfer due to
the building envelope [19, 20, 21]. Further precision has been incorporated
into the RC network by using dynamically estimated capacitances [22] and
second-order thermal network models [23]. The physics-based Three resistor-Two
Capacitance (3R2C) model has been shown to be sufficiently accurate for
approximating building models [24] based on parameters with physical meaning.
In the current work, we build upon 3R2C models to enable our modeling
framework.

The lumped capacitance method has been used to estimate other factors
affecting the thermal zone heating and cooling requirements. The thermal mass
of the zone, in particular, can be modeled with an RC network [16, 25]. When
these values are significant, long-wave radiation [15] and convective heat
between zones [26] can be modeled as black-box data-driven estimations.
Internal loads, especially, have been estimated using a filtering method in
conjunction with the RC method for the surfaces [27].

The existing literature related to the RC method is inadequate for modeling
large-scale problems due to high computational cost. This computational cost
is exacerbated due to increase in the number of inputs. Thus, an essential
contribution is in the development of a computational method to expedite the
UQ using the idea of divide and conquer. We outline such approach in this
paper.

###  2.2 Uncertainty Quantification in BEMs

Uncertainty Quantification (UQ) is becoming more prevalent in the BEM domain,
as building simulation software and methods continue to evolve [28]. Research
in the field of UQ in BEM is growing at a fast pace. As explained earlier,
uncertainty in a BEM can arise due to the modeling process, or uncertainty in
the parameters. All BEMs necessarily make assumptions to simplify the modeled
building. Simplifying the assumptions in the building energy simulation domain
causes the associated uncertainty to be often ignored. Few existing works have
focused on UQ in BEMs. For example, Sun et al. [10] have focused on
incorporating the uncertainty in the solar irradiation calculation and its
effects on the results of BEM simulation. Additionally, post-processing
techniques have been developed for incorporating UQ in BEMs [9].

Recent studies have used parameter estimation for UQ in BEMs. Both physics-
based and surrogate models have been studied to optimize building performance
under uncertainties by simulation methods [29]. It has been shown that even
moderately variable parameters can have a significant effect on the overall
BEM uncertainty, especially, when the small-scale models are combined to
generate a large-scale BEM, as demonstrated in multi-building residential
district models [30, 31]. Sensitivity analysis has often been used to
determine the impact of uncertainty in parameters on the performance of the
building energy model [32, 33]. Most of the work in sensitivity analysis
focuses on identifying a few most influential parameters that can be used for
further investigation during UQ. Additionally, most of the reported works are
often limited to evaluating only particular aspects of a building and not the
whole building [31].

UQ in BEMs is also performed in Stochastic Model Predictive Control (SMPC)
frameworks that utilize the dynamic state-space equations. Oldewurtel et al.
have discussed how SMPC framework can be used to account for the uncertainty
in weather predictions for building energy modeling purposes [34]. In another
related work, an SPMC has been developed to optimize building energy usage and
has been shown to outperform existing Rule-based Control (RBC) [35] framework.
SPMC has also been used for large-scale building (a large number of state
variables) in the presence of uncertainties in weather prediction [36].
Privara et al. [37] have used an identified state-space model of a real
building to estimate the optimal energy consumption. The model developed by
Privara et al. does not include the effect of internal loads or solar gains.
An adaptive MPC has been developed incorporating uncertainties in a wide range
of parameters including building materials, their thermal properties, and HVAC
parameters [38].

A simulation over an extensive range of values is required with a minimum
computational expense to enable UQ. To reduce the computation time, techniques
such as quasi-random sampling [12] and pre-processing historical data in
modeling predictive controls [39] have been used. Similarly, other sampling-
based methods or quadrature-based methods can also be applied [40]. When
applied to large-scale problems with many input parameters, the performance of
existing methods is computationally inefficient. We address these shortcomings
in this paper through the use of Weakly Connected Subsystems based UQ method
[13, 18, 17].

##  3 Methodology

###  3.1 Proposed Framework

Figure 1: Schematic of the WCS identification-based UQ to estimate thermal
load in a large-scale BEM. Weather Map [41]

Figure 1 depicts the overall UQ framework for a large-scale BEM. Given an
office/school building (Figure 1(a)), the geometry and the thermal properties
of the building are assessed to formulate the state-space equation model of a
dynamical system (Figure 1(c)). This state-space equation is formulated using
the concept of RC network (Figure 1(b)). The output equation is framed
depending on the available measurement. A graph-theoretic representation is
adopted to model the thermal network as an undirected graph (Figure 1(d)). The
state-space equation and the initial uncertainty information is used to
quantify the adjacency information for the undirected graph (Figure 1(e)). A
suitable graph clustering algorithm (Louvain modularity optimization) is then
implemented to identify the Weakly Connected Subsystems or WCSs (Figure 1(f)).
Subsequently, the input formulation involving the weather (ambient and soil
temperature) is obtained, and the solar gain for the exterior surfaces and
HVAC airflow and temperatures are calculated, using eQUEST (Figure 1(g)). The
statistical properties of the state variables are propagated through each
WCSs, which are updated based on the measurement availability. The updated
statistical properties give us the estimated parameters such as non-envelope
load and the solar gain for internal surfaces, along with the zonal
temperatures, and the surface temperatures. In the subsequent sections, the
individual components of the overall framework are described in detail.

###  3.2 Parameters of BEM: Geometry and Thermal Properties

The simulation requires the specification of the building geometry and surface
properties. A BEM is divided into thermal zones, based upon the actual HVAC
system. Each area with an individual temperature sensor is considered as a
thermal zone. Within the zone, the following properties of each surface is
determined:

  1. 1.

Exterior Surface: Indicates surface is exposed to outside air conditions.

  2. 2.

Interior Surface: Indicates surface is adjacent to another zone and exposed to
the conditions of that zone.

  3. 3.

Underground Surface: Indicates surface that is below the ground level and
exposed to soil temperatures.

  4. 4.

Adjacent Zone: If interior surface, the identifier of the adjacent zone
connected to the interior surface.

  5. 5.

Area ($A^{j}$): The net area of the surface, for windows $w=1$ to $c$ on
surface $j$ is derived as:

$A^{j}=A_{gross}^{j}-\sum_{w=1}^{c}{A_{win}^{w}}$ (1)

  6. 6.

Resistance ($R^{j}$): A measurement of the heat flow across a surface at
steady state conditions, calculated as follows:

$R^{j}=\sum_{g=1}^{h}\frac{l^{g}}{k^{g}A^{j}}$ (2)

where $l$ is material thickness and $k$ is thermal conductivity for material
$g=1$ to $h$ in the construction of surface $j$.

  7. 7.

Capacitance ($C^{j}$): A measurement of the amount of heat input required to
raise the temperature of a material. Capacitance encapsulates the capability
of a material to store heat and delay heat transfer across the material.
$C^{j}$ is calculated as follows:

$C^{j}=\sum_{g=1}^{h}\frac{\rho^{g}C_{p}^{g}l^{g}A^{j}}{2}$ (3)

where $l$ is material thickness, $\rho$ is density, and $C_{p}$ is specific
heat for material $g=1$ to $h$ in the construction of surface $j$. This
formula assumes half the capacitance is applied to each the inside and the
outside faces of surface $j$.

  8. 8.

Window Resistance ($R_{win}^{w}$): Similar to $R^{j}$. Windows are modeled
using only the specified resistance. The thermal capacitance effects are
negligible [24].

  9. 9.

Inside Surface Conductance ($h_{i}^{j}$): A measure of heat flow at the
interior side of a surface $j$. Assumed to be static based on surface location
and position [24, 42]

  10. 10.

Outside Surface Conductance ($h_{o}^{j}$): Similar to $h_{i}$, but at the
exterior side of a surface $j$.

  11. 11.

Zone Air Mass ($M_{k}$): The air mass of zone $k$ at design conditions.

  12. 12.

Specific Heat of Air ($c_{p}$): The specific heat of air at design conditions.

To calculate input parameters, explained later, the following properties of
the building are also assumed:

  * •

Location: The geographic coordinates of the overall building.

  * •

Orientation: The cardinal direction and angle of the surfaces.

  * •

Shading devices: The size and position of any permanent sun shades.

Surface properties are critically important for the BEM since the heat
transfer occurring through the surfaces is a fundamental concept for these
calculations. Boundary conditions, in particular, determine the extent that
the weather conditions impact the zone. In the case of underground surfaces,
there is no adjacent airspace, so there is no outside surface conduction
$h_{o}$. Instead, a fictitious conductance ($1/{R_{fic}}$) is used (explained
later in Section 3.4).

The values of the material properties are determined by testing. In a
building, the properties are assumed based upon either specific data published
from the manufacturer, or, if unknown, typical values of materials are
compiled in subject references [42]. In some cases, values of typical surface
constructions as a whole are published.

###  3.3 Heat Balance Method

The heart of the BEM problem is the fundamental heat balance equation for zone
$k$ [43]:

$\displaystyle
C_{k}\frac{dT_{k}}{dt}=\sum\dot{T}^{int}+\sum_{j}h_{j}A_{j}(T_{amb}-T_{k})+\sum_{i}\dot{m}_{mix}c_{pa}(T_{k,i}-T_{k})+\dot{m}_{inf}c_{pa}(T_{\infty}-T_{k})+\dot{m}c_{pa}(U_{sa}-T_{k})$

(4)

where:

  * •

$C_{k}\frac{dT_{j}}{dt}$ = energy stored in the zone air

  * •

$\sum\dot{T}^{int}$ = sum of the convective non-envelope internal loads

  * •

$\sum_{j}h_{j}A_{j}(T_{amb}-T_{z})$ = convective heat transfer from the zone
surfaces $j$

  * •

$\sum_{i}\dot{m}_{mix}c_{pa}(T_{z,i}-T_{k})$ = heat transfer due to inter-zone
air mixing with zones $i$

  * •

$\dot{m}_{inf}c_{pa}(T_{\infty}-T_{k})$ = heat transfer due to infiltration of
outside air

  * •

$\dot{m}c_{pa}(U_{sa}-T_{k})$ = heat transfer due to the output of the HVAC
system into the zone

The definitions of the relevant variables are listed in the Nomenclature and
discussed later. Equation 4 includes a number of simplifying assumptions [42,
44]. These assumptions are typical of most modeling techniques, and allow the
use of a linear equation:

  1. 1.

Air temperature in the zone is well-mixed. The variation of the zonal
temperatures is zero.

  2. 2.

The surface temperature and irradiation is uniform.

  3. 3.

Surface radiation is diffuse.

  4. 4.

One-dimensional heat conduction through the surface construction.

Some other assumptions made for simplicity, which are typical in many BEMs,
are as follows:

  1. 1.

No inter-zone mixing. The supply and return airflows are balanced per zone, so
the overall pressure differential of a building is approximately zero. The
majority of spaces are separated by doors. Therefore, there is minimal air
transfer between thermal zones, and this expression can be eliminated from the
heat balance equation.

  2. 2.

Simplified infiltration. As discussed later, infiltration is not calculated
directly and is included with other internal loads in the space as part of the
lumped non-envelope load, eliminating this expression from the heat balance
equation.

  3. 3.

No long-wave radiation. The air is assumed to be transparent. Long-wave
radiation exchanges between surfaces can be ignored [42].

Latent loads, such as moisture due to the occupants, are not addressed in the
BEM. Humidity is not often directly controlled, but instead is handled as a
side-effect of the sensible (dry-bulb) cooling. Humidity primarily affects
occupant comfort and the HVAC system efficiency, although a small amount of
heat transfer occurs due to moisture in the air, with minimal impact on the
zonal temperatures [42]. Any simplification due to the model itself will be
captured in the lumped non-envelope load estimation. Thus latent loads are
ignored.

The reduced-order model as described in Section 3.4 is based upon the
simplified heat equation as follows:

$C_{k}\frac{dT_{k}}{dt}=\sum\dot{T}^{int}+\sum_{j}h_{j}A_{j}(T_{amb}-T_{k})+\dot{m}c_{pa}(U_{sa}-T_{k})$
(5)

###  3.4 State Space Equation Using RC Network

As explained previously, the resistor-capacitor estimation method reduces the
building surfaces into a discrete number of resistances and capacitances. A
three resistance, two capacitance model (3R2C) is the most used one. This type
of thermal load estimation model has been found to estimate the effects of the
building surfaces on a thermal zone with sufficient accuracy [24]. Using this
method, the heat balance equation can be decomposed into distinct parts that
enable us to perform the load optimization. The RC network for the surface
constructions used for modeling purposes is shown in Figure 2.

Figure 2: RC network for surface constructions

Summarizing the assumptions in Section 3.3 with the RC network model, the
relevant heat balance equations are as follows for surface $j$ and window $w$
in zone $k$ [27]:

$\begin{array}[]{c}\displaystyle
C^{j}\frac{dT_{o}^{j}}{dt}=h_{o}^{j}A^{j}(T_{amb}-T_{o}^{j})+\frac{T_{i}^{j}-T_{o}^{j}}{R^{j}}+U_{o}^{j}\\\
\displaystyle
C^{j}\frac{dT_{i}^{j}}{dt}=h_{i}^{j}A^{j}(T_{k}-T_{i}^{j})+\frac{T_{o}^{j}-T_{i}^{j}}{R^{j}}+T^{j}_{q}\\\
\displaystyle
Q_{structure}=h_{i}^{j}A^{j}(T_{i}^{j}-T_{k})+\frac{T_{amb}-T_{k}}{R_{win}^{w}}\end{array}$
(6)

Underground zones are treated as special cases, as their exterior surfaces are
not directly exposed to the outdoor temperatures. Instead, the majority of the
heat transfer in the underground surfaces takes place at the exposed perimeter
region of the surface. The constructions of the exterior surfaces of
underground zones are modeled a no-mass R-value layer to avoid any
overestimation of heat loss [45]. The overall effective R-value is based on
published data regarding the heat transfer through underground surfaces,
taking into account the depth of the surface and the location and thickness of
the insulation [46].

$R_{fic}^{j}=R_{eff}^{j}-R^{j}$ (7)

where $R_{eff}$ is the effective R-value, and $R$ is calculated as described
in Section 3.2, including a 0.3 meter layer of soil as part of the surface
construction.

The modified RC network for an underground surface is shown in Figure 3.

Figure 3: Modified RC network for underground surfaces

Using this method, one can use consistent equations for any surface type. For
the interior zones, the adjacent zone temperature $T_{out}$ is a calculated
state variable. The $T_{out}$ for exterior surfaces is the ambient temperature
$T_{amb}$ that is available from weather data. The ground or the soil
temperature is used as the ambient temperature $T_{amb}$ of underground
surfaces.

Combining Equations 6, the heat balance equation is modified as a system of
Ordinary Differential Equations (ODE). The state space matrix derived from the
ODE for the BEM can be written as:

$\begin{array}[]{ll}\dot{\textbf{x}}_{t}&=H_{t}\textbf{x}_{t}+B\textbf{u}_{t}\\\
\textbf{z}_{t}&=C\textbf{x}_{t}\end{array}$ (8)

The state-space Equation 8 comprises of the zonal temperature $T_{k}$ and the
inner and outer surface temperatures $T_{i}^{m}$ and $T_{o}^{m}$’s. In
addition, for each zone, the solar gain $T_{q}^{j}$ for each inner surface and
the non-envelope load $T_{k}^{ne}$ are modeled. Thus,
$\textbf{x}\in\mathbb{R}^{N}$ can be decomposed into the collection of zonal
variables as
$\textbf{x}=\\{\textbf{T}_{1},\textbf{T}_{2},\ldots\textbf{T}_{N}\\}$. Each
$\textbf{T}_{k}\in\mathbb{R}^{n_{k}}$, $k=1$ to $N$ represents the collection
the zonal temperature $T_{k}$, inner and outer surface temperatures for the
$m_{k}$ surfaces, and the three load variables in a particular zone. Hence,
$\textbf{T}_{k}$ can be written as:

$\textbf{T}_{k}=\\{T_{k},T_{o}^{1},T_{i}^{1},T_{q}^{1},\ldots,T_{o}^{m_{k}},T_{i}^{m_{k}},T_{q}^{m_{k}},T_{k}^{int}\\}\hskip
14.22636ptn_{k}=2+3m_{k}$ (9)

And,

$N=\sum_{k}^{n}n_{k}=\sum_{k}^{n}3m_{k}+2=2n+3\sum_{k}m_{k}$ (10)

$T_{k}^{ne}$ and $T_{q}^{j}$ for each surface are modeled as state variables
and are estimated with available measurement. Each zone temperature is
determined by the conduction through the surface constructions, as well as the
HVAC airflows and temperatures, surface solar gains, and internal non-envelope
loads. Surfaces adjacent to another zone share properties with that zone.
Zones with differing occupancies and internal loads are expected to have the
most zonal interactions. The components of the zonal variable $\textbf{T}_{k}$
are modeled as [27]:

$\begin{array}[]{ll}\dot{T}_{k}&=\displaystyle\left[-\frac{\dot{m}_{k}}{M_{k}}-\frac{\sum_{m}A^{m}h_{i}^{m}}{M_{k}c_{pa}}-\frac{\sum_{w}\frac{1}{R_{win}^{w}}}{M_{k}c_{pa}}\right]T_{k}+\frac{\sum_{m}A^{m}h_{i}^{m}T_{i}^{m}}{M_{k}c_{pa}}+\frac{1}{M_{k}c_{pa}}T^{int}_{k}+\frac{\frac{1}{R_{win}^{w}}}{M_{k}c_{pa}}T_{amb}+\frac{\dot{m}_{k}U_{k}^{sa}}{M_{k}c_{pa}}\\\
\dot{T}^{j}_{o}&=-\displaystyle\left[\frac{h_{o}^{j}A_{j}}{C^{j}}+\frac{1}{R^{j}C^{j}}\right]T_{o}^{j}+\frac{1}{R^{j}C^{j}}T_{i}^{j}+\frac{h_{o}^{j}A_{j}}{C^{j}}T^{j}_{out}+\frac{U^{j}_{o}}{C^{j}}\\\
\dot{T}^{j}_{i}&=\displaystyle\frac{h_{o}^{j}A_{j}}{C^{j}}T^{j}_{k}+\frac{1}{R^{j}C^{j}}T_{o}^{j}-\displaystyle\left[\frac{h_{i}^{j}A_{j}}{C^{j}}+\frac{1}{R^{j}C^{j}}\right]T_{i}^{j}+\frac{T^{j}_{q}}{C^{j}}\\\
\dot{T}_{k}^{int}&=0,\hskip 8.53581pt\dot{T}^{j}_{q}=0\hskip
42.67912ptj=1,\ldots,m_{k}\end{array}$

(11)

where:

  * $\dot{m}_{k}$ = mass flow of supply air from HVAC system to zone $k$

  * $U_{k}^{sa}$ = temperature of supply air from HVAC system to zone $k$

  * $M_{k}$ = mass of air in zone $k$

  * $c_{pa}$ = specific heat of air

  * $T_{k}^{ne}$ = heat gain from non-envelope internal loads

The measurement to this system $\textbf{z}\in\mathbb{R}^{n}$ are the zonal
temperatures $T_{k}$’s and is characterized by the observation matrix
$C\in\mathbb{R}^{N\times n}$.

###  3.5 Input Formulation

The majority of the inputs for Equation 11 can be obtained directly for zone
$k$ at each time step.

For the outdoor conditions, Typical Meteorological Year (TMY) dataset is used
for the nearest city [47]. The TMY data is not representative of any
particular year but is developed to represent the typical conditions of a
given location. Estimated average monthly ground temperatures at different
geographic locations are also available in the TMY dataset. For a building
with a Building Management System, the $T_{amb}$ conditions would be available
with the HVAC trend data.

Solar heat gain on opaque exterior surface $U_{o}^{j}$, $j=1$ to $m_{k}$ is
calculated based on the buildings location and orientation in conjunction with
the solar radiation intensity [43]:

$U_{o}^{j}=\alpha^{j}(I_{b}\frac{A_{s}^{j}}{A^{j}}\cos\theta^{j}+I_{s}F^{j}_{ss}+I_{g}F^{j}_{sg})$
(12)

The intensities of radiation are not measured directly but can be calculated
from measured radiation (direct normal and diffuse horizontal) using the
luminous efficacy models [48]. The Perez model is the basis of many modern
solar models, including the approximate values given in the TMY dataset [47].
Since the real-time solar data is not available, the calculated input values
for the surface solar heat gain are not dependent on actual weather
conditions. Thus, for simplicity, the solar load values $U_{o}^{j}$ calculated
by eQUEST are used in the model. Only exterior surfaces are subject to direct
solar radiation and are assumed to have outer solar heat gain.

Calculating the HVAC system inputs for the heat balance equation requires
specific information on fan curves, coil properties, pump curves, pressure
losses, control strategies, capacity curves, etc. Furthermore, there is a
feedback loop relationship between the HVAC output and the building load
components [15]. Because of the complexity in calculating the air-side HVAC
parameters and because the data required for this analysis are typically
readily available with a robust building controls system, the hourly air
temperatures and flows calculated by eQUEST are used in our state equations.

###  3.6 Estimation of Solar Gain and Non-Envelope Loads

In this work, the solar heat gains $U_{o}^{j}$ on surface $j$ and Non-Envelope
Load $T^{j}_{q}$ are estimated through the model (as explained in Section
3.7). In a fully-functional BEM, such as simulated with eQUEST, both must be
assigned approximate values to perform the simulation.

In such a BEM, the solar effect on the interior walls $U_{o}^{j}$ due to
windows is calculated using the solar heat gain coefficient (SHGC). First, the
total heat gain through the windows $Q_{isol}^{k}$ for windows $w=1$ to $c$ in
zone $k$ is calculated [42].

$Q_{isol}^{k}=\sum_{w=1}^{c}(SHGC^{w})A^{w}E_{t}^{w}$ (13)

where:

  * $SHGC^{w}$ = Solar Heat Gain Coefficient of window $w$

  * $A^{w}$ = area of window $w$

  * $E_{t}$ = incident total irradiance of window $w$

The SHGC of the window represents the portion of solar radiation transmitted
directly through the window, as well as the absorbed solar radiation. This
value is typically provided by the window manufacturer or can be assumed to be
based on the window properties [42]. From $Q_{isol}^{k}$, the solar flux
$T_{q}^{j}$ on the interior surface $j$ can be estimated, using one of several
approaches. One such approach, is to assume that all radiation first hits the
floor of zone $k$, and is reflected evenly across all the surfaces $r=1$ to
$s$ in zone $k$ [27]. In the current work, this internal solar gain $T_{q}$ is
estimated as a part of the UQ framework.

Non-envelope internal loads $T_{k}^{ne}$, which refer to loads in the space
not related to the building exterior walls, are also unknown and are estimated
through the WCSs-based UQ framework. in a full BEM, load schedules are
assumed, such as used by He at al. [15]. It is unlikely, however, that any
physical building would follow the exact prescribed schedules with precision.
Thus an optimal estimation method becomes necessary, especially due to the
involved high-dimensional system.

In a full simulation, the airflow due to infiltration must be assumed. The
infiltration rate at any given time is based upon the pressure differential
between the building and the exterior environment, as well as the effective
leakage area of the building. These parameters depend on air temperature, air
density, wind speed, and wind direction, and the resulting differential is
modified by the nature of the building openings. To obtain an accurate,
effective leakage area, a blower door test is required, pressurizing the space
to a specified value [42]. For a large building, this is a costly process.
Consequently, the estimated peak infiltration rate is usually simply based on
an assumed construction tightness. Because of the difficulty obtaining
accurate infiltration data, the effect due to infiltration is also estimated
with the other internal loads as part of the non-envelope load estimation.

Most thermal zones also have furniture that provides thermal mass and
additional surfaces for radiation. Specific information on the furniture in
the zone, especially thermal properties, is difficult to obtain. Wang and Xu
[25] use a 2C2R model to estimate the internal mass of a zone using a genetic
algorithm; however, the resulting parameters have no physical meaning other
than an assumed lumped mass for estimation purposes. Therefore, the state
space equations in this work do not directly include any effects due to
internal mass. Internal mass effects are implicitly included in the non-
envelope load estimation.

In the presented case study, there are zone-level HVAC units such as hot water
baseboard radiation. These baseboard units provide additional radiant or
convective heat to the zone. The Building Management System (BMS) controls and
monitors these baseboard units that work in conjunction with the supply air.
The baseboard units are included as a special case of the internal load. Like
the HVAC parameters, the eQUEST hourly heat output values for the baseboard
has been used. In the case study, the BMS tracks the operation of the control
valve, which can be used with the design flow rates and actual hot water loop
temperatures to approximate the unit heat output. This value is then simply
subtracted from the calculated internal load.

###  3.7 Identification of Weakly Connected Subsystems and Optimal Estimation

Optimal estimation of the unknown variables $T_{i},T_{o},T^{int}$ and $T_{q}$
in the problem detailed in Equation 11 refers to solving the following
minimization problem

$\displaystyle\min_{\hat{\textbf{x}}_{k}}E\left[||\textbf{x}_{k}-\hat{\textbf{x}}_{k}||_{2}\right]$
(14)

The above term refers to the error in the a posteriori state estimation. For a
given measurement $\textbf{z}_{t}$, solution to the problem is same as solving
the well known Filtering problem. Consider the solution at time $t-1$ as
$\hat{\textbf{x}}_{t-1}$ and covariance $\Sigma_{t-1}$, the estimate
$\hat{\textbf{x}}_{t}$ and covariance $\Sigma_{t}$ are given as,

$\begin{array}[]{l}\hat{\textbf{x}}_{t}=\hat{\textbf{x}}_{t|t-1}+K\textbf{y}_{t}\\\
\Sigma_{t}=\Sigma_{t|t-1}-K(R+C\Sigma_{t|t-1}C^{T})K^{T}\end{array}$ (15)

where $K$ is known as the Kalman gain and
$\textbf{y}_{t}=\textbf{z}_{t}-\hat{\textbf{x}}_{t|t-1}$ is known as the
measurement residual. The a priori estimates $\hat{\textbf{x}}_{t|t-1}$ and
$\Sigma_{t|t-1}$ are one step solution to the Equation 8 depending on
$\hat{\textbf{x}}_{t-1}$ and covariance $\Sigma_{t-1}$. The expression for the
Kalman gain $K$ is given as [49],

$K=\Sigma_{t|t-1}C\left(R+C\Sigma_{t|t-1}C^{T}\right)^{-1}$ (16)

Due to the high dimensionality of the problem involving a large number of
state variables $N$, the series of matrix operations becomes computationally
expensive. The above filtering problem is solved through Identification of
Weakly Connected Subsystems (WCSs) [13, 17, 18] to increase the computational
efficiency. The WCS-based method is effective in solving high-dimensional UQ
problems involving linear and non-linear filtering problems.

WCSs for the state variable $\textbf{x}\in\mathbb{R}^{N}$ are the countable,
mutually exclusive and exhaustive partitions
$\textbf{y}_{j}\in\mathbb{R}^{n_{j}}$, such that the following relation holds

$\begin{array}[]{l}\displaystyle
P(\textbf{x}_{t})=\prod_{j}P(\textbf{y}_{j_{t}})\\\
\sum_{j}n_{j}=N\end{array}$ (17)

The index $t$ represents that the relation is invariant under the
transformation given by Equation 8 for a time-period $t\in[0,T)$. Performing
such decomposition of $\textbf{x}_{t}$ enables faster UQ by solving parallel
subproblems given as:

$\dot{\textbf{y}}_{j}=H_{j}(t)\textbf{y}_{j}+B_{j}\textbf{u}_{j}$ (18)

The solutions to each subproblem in Equation 18 is given by the following
continuous time Kalman filter [50]:

$\begin{array}[]{l}\displaystyle\dot{E(\textbf{y}_{j})}=H_{j}(t)E(\textbf{y}_{j})+B_{j}\textbf{u}_{j}+K_{j}(z_{j}-C_{j}E(\textbf{y}_{j_{t}}))\\\
\displaystyle\dot{\Sigma_{j}}=H_{j}(t)\Sigma_{j}+\Sigma_{j}H_{j}(t)^{T}-K_{j}R_{j}K_{j}^{T}\\\
\displaystyle K_{j}=\Sigma_{j}H_{j}^{T}R_{j}^{-1}\end{array}$ (19)

This invokes the use parallel computation for both the one-step a priori
estimation and as well the use of Kalman filter for the a posteriori
estimation. The state estimates $\textbf{x}_{t}$ and $\Sigma_{t}$ are computed
from the WCSs using direct sum of the vector spaces as:

$\begin{array}[]{l}E(\textbf{x}_{t})=E(\textbf{y}_{1_{t}})\oplus
E(\textbf{y}_{2_{t}})\oplus\ldots\oplus E(\textbf{y}_{m_{t}})\\\
\Sigma_{t}=\displaystyle\bigoplus_{j}\Sigma_{j_{t}}=\text{diag}\left(\Sigma_{1_{t}},\ldots,\Sigma_{m_{t}}\right)\end{array}$
(20)

The normalized symmetrized adjacency matrix derived from the state-space
matrix $H=\left(h_{ij}\right)$ in Equation 8 [13] is given as:

$W=0.5(D^{-1}A_{abs}+A_{abs}^{T}D^{-1})$ (21)

where, $A_{abs}=\left(|h_{ij}|\right)$ and $D$ is the corresponding degree
matrix of $A_{abs}$. The clusters are identified using Louvain method of
community detection [51]. The method identifies weakly connected components in
a weighted graph by maximizing modularity function defined as:

$Q=\frac{1}{2m_{q}}\displaystyle\sum_{i,j}\left[W_{i,j}-\frac{k_{i}k_{j}}{2m_{q}}\right]\delta(c_{i},c_{j})$
(22)

where, $m_{q}=\sum_{i,j}W_{i,j}=N/2$, $k_{i}=\sum_{i}W_{i,j}$ and $c_{i}$ is
the partition to which $i^{\text{th}}$ state belongs. The delta function
$\delta(c_{i},c_{j})$ is

$\delta(c_{i},c_{j})=\begin{cases}1&c_{i}=c_{j}\\\
0&\text{otherwise}\end{cases}$ (23)

Maximizing $Q$ gives the values of $c_{i}$’s, $i=1$ to $N$ and hence
determines the cluster structure. In the next section, details pertaining to a
building used for BEM modeling purposes are outlined.

##  4 BEM Details

###  4.1 Building Description

The case study building used in this work to illustrate the efficacy of the
outlined UQ framework is a College/University building in Central New York,
United States (see Figure 4). The facility is an existing 4-story 5,050 square
meter building with a mechanical penthouse. It is comprised of primarily
classrooms and offices, student lounges, conference rooms, observation rooms
and as well as other support spaces. The lowest floor comprises of underground
zones, and the northeast portion of the building is attached to an adjacent
structure.

Figure 4: Schematic of College/University building in Central New York used
for creating the case study BEM

The design of the building HVAC system informs the eQUEST calculations for the
inputs for $U^{sa}$ and $\dot{m}$. The individual spaces are conditioned at
the zone level by fan-coil units with fin-tube radiation on the perimeter. A
dedicated outdoor air system supplies tempered ventilation air directly to the
fan-coil units. The hot and chilled water coils are supplied from the campus
plants. Network data rooms are conditioned separately with a variable
refrigerant flow heat pump system.

Surface constructions determine the model parameters $R$ and $C$. The brick
and concrete block envelope has a combination of rigid and spray-applied
insulation at the exterior walls, and the reflective membrane roof is concrete
deck topped with rigid insulation. The windows are tinted high-performance
glazing with sunshades on the southern exposure.

To understand $T^{int}$, the lighting, occupancy, and equipment must be
considered. The building uses high-efficiency LED lighting with occupancy
sensors and daylighting controls. The building operates Monday through Friday
from 7 am to 10 pm, with the expected use of a typical university building.

There are a total of 132 thermal zones in the building, including
unconditioned plenum spaces above the ceilings, for 61 directly conditioned
zones. In the modeled building, there are 668 interior surfaces, 69
underground surfaces (including slab-on-grade floors), 124 exterior surfaces
(including roofs), and 80 windows. The second-floor thermal zones are shown in
Figure 5.

Figure 5: Second Floor Thermal Zones of the BEM

Four zones are specifically considered to illustrate the details of our
methodology. These are a ground floor perimeter zone with underground walls
(Zone 1), a first floor exterior zone with west facing windows (Zone 47), a
second floor zone with no exterior surfaces (Zone 80), and a third floor zone
with a roof and south and east exterior walls (Zone 97). The number of input
parameters needed per zone to calculate only conduction through building
surfaces is listed in Table 1. Each zone also has dynamic parameters
$T_{int}^{k}$ for lighting, occupancy, equipment load, and infiltration, as
well as the inputs $\dot{m}_{k}$ and $U_{sa}^{k}$ for the HVAC system and
assumptions by the BEM program. Each of these inputs has uncertainty
associated with them. As these zones represent only 4 out of the 132 zones, it
is clear that this is a high-dimensional problem that requires a scalable UQ
method.

Table 1: BEM Inputs per Zone - Conduction through Surfaces Only

| Number of | Number of | Number of |  | Number of | Number of | Number of
---|---|---|---|---|---|---|---
| Exterior | Interior | Underground | Number of | Construction | Adjacent | BEM Input
| Surfaces | Surfaces | Surfaces | Windows | Materials | Zones | Parameters
Zone 1 | 0 | 6 | 4 | 0 | 11 | 2 | 52
Zone 47 | 2 | 12 | 0 | 2 | 11 | 3 | 80
Zone 80 | 0 | 6 | 0 | 0 | 10 | 4 | 48
Zone 97 | 2 | 8 | 0 | 6 | 14 | 4 | 81

The main floors of the modeled building are partially attached to an adjacent
building. However, the adjacent building is controlled independently of the
modeled building, and the properties of the adjacent building are unknown.
Therefore, the surface boundary conditions with the adjacent building are
assumed to be adiabatic, with no heat transfer of any kind between the two
buildings.

The modeled building is used to demonstrate our proposed methodology largely
due to its size and complexity. With 61 directly conditioned thermal zones,
one can explore the differences between and interconnectivity of perimeter
zones and core zones. There are several exterior wall construction types, as
the first floor differs from the upper floors, and the underground walls
require a modified approach to calculate the heat transfer. Complex lighting
controls create a further challenge, which compels the internal load
estimation and UQ. There are a large number of input parameters that
necessitates the use of a scalable UQ method.

###  4.2 Full-Scale BEM Details

To provide data to our reduced-order model, this building has also been
modeled in eQUEST (version 3.65) [6] following the methodology laid out in
ASHRAE 90.1-2013 (Appendix G [52]). All known parameters of the building, as
discussed in Section 4, is provided as an input to the software program.
eQUEST uses a white-box approach to simulate building energy. At first, it
calculates the zone load at each time step. This is followed by converting the
zone load into the required load and conditions of the HVAC system. The
calculated load is fed back into the load calculation [53]. The process for
calculating building energy use is depicted in Figure 6.

Figure 6: Schematic for calculation of building energy

Since the modeled building has a digital building management system, one has
access to the HVAC airflows and temperatures, water flow rates and
temperatures, zone temperatures, as well as a variety of other control points.
A list of actual measurements available for the fan-coil units in this
building from the BMS is listed in Table 2. The input to Equation 8 is
obtained directly from the simulation of the eQUEST model.

Table 2: Fan Coil Unit Control Points | Analog | Analog | Digital | Digital
---|---|---|---|---
Description | Input | Output | Input | Output
Fan Enable/Disable |  |  |  | X
Fan Status |  |  | X |
Space Temperature | X |  |  |
Heating Valve Position |  | X |  |
Cooling Valve Position |  | X |  |
Leaving Air Temperature | X |  |  |
Filter Pressure Differential Sensor |  |  | X |
Condensate Alarm |  |  | X |

In the eQUEST model, internal loads are assumed using hourly schedules.
Lighting peak loads are calculated directly from the known lighting layout.
However, the occupancy and receptacle loads are estimated. Building code
provides an approximation for the maximum occupancy for each space type [54].
Receptacle loads can be highly variable, especially with today’s proliferation
of electronics. Typical representative load values have been compiled in
COMNET [55]. These peak load conditions are then scaled hourly based on
estimated schedules. The ASHRAE 90.1 User’s Manual provides typical fractional
schedules for occupancy, zone lighting, and receptacle loads [46] that have
been modified in this work. The building lighting in the modeled building is
assumed to be controlled by a combination of daylight harvesting controls,
occupancy sensors, dimming switches, and manual on/off controls. Using simple
reduction factors, ASHRAE 90.1 provides guidance on the expected effect of
some of the lighting controls [52]. However, outcomes are not representative
of a particular building, and it is likely that in practice, the actual hourly
lighting power density would be substantially different than the predicted.

##  5 Results

###  5.1 Adjacency and Cluster Analysis

The BEM described in Section 4 contains 2649 state variables (Involving
$\\{T_{k},T_{o},T_{i},T_{q},T^{int}\\}$), with 1125 input variables. The
state-space equation (in Equation 8) of the BEM is based on the thermal
properties of the materials of the building and the interaction of the zonal
and surface temperatures. The adjacency matrix $W$ is a reflection of coupling
along with the presence of any inter-zone coupling. The variable $T_{out}^{j}$
in Equation 11 can be both the ambient temperature $T_{amb}$ or the
temperature of adjacent zone. Thus, the inter-zone coupling (if present)
depends on the magnitude of $h^{j}_{o}A_{j}/C_{j}$. Additionally, the
adjacency matrix is highly sparse. This makes a visual representation of the
whole $W$ very difficult to interpret. Results aggregated from the four
thermal zones listed in Table 1 representing diverse zonal conditions: Zone 1,
Zone 47, Zone 80 and Zone 97 are discussed next.

(a) Zone 1

(b) Zone 47

(c) Zone 80

(d) Zone 97

Figure 7: Adjacency Information for (a) Zone 1, (b) Zone 47, (c) Zone 80 and
(d) Zone 97

Figure 7 shows the sub-matrices of the Normalized Adjacency matrix $W$
corresponding to the four zones. The zones show the involved variables
$T_{k},T_{o},T_{i},T^{int}$. Zone 1 and 80 do not have exterior surfaces. Thus
the variable $T_{q}$ is absent from the corresponding adjacency figures
(Figures 7(a) and (c)). Figure 8 demonstrates the inter-zone coupling between
these four zones. The plot shows no interaction between these four zones.
However, this result does not conclude that there is zero interaction between
other zones.

Figure 8: Adjacency Matrix for Zone 1, Zone 47, Zone 80 and Zone 97 combined

The cluster matrix is depicted in Figure 9.

Figure 9: Cluster Matrix

The application of the Louvain modularity optimization-based clustering
algorithm on the derived $W$ matrix for 2649 states resulted in identification
of $m=97$. The number of identified clusters ($m=97$) is less than the number
of zones(=132). However, the cluster structure works in accordance with the
adjacency information identified in Figure 7 and 8. It can be observed in
Figure 9 that consecutive states participate in clusters. Intuitively it can
be concluded that the clustering algorithm either recognizes a single zone or
a group of zones as a cluster. In most cases it ignores the inter-zone
coupling. In the next section, the whole system is analyzed along with the
estimation of the thermal loads via these identified clusters or WCSs. To
compare the performance of the WCS-based UQ, the whole BEM is also analyzed
without any clustering.

###  5.2 Estimated Temperatures and Thermal Loads

For the simulation, a measurement noise is assumed for each of the zonal
temperature as $R=\text{diag}\\{r_{1},r_{2},\ldots,r_{132}\\}$ where
$\\{r_{1},r_{2},\ldots,r_{132}\\}$ are randomly generated numbers between 0
and $0.5$. This corresponds to the accuracy of typical room temperature sensor
[56]. Before displaying the estimated thermal loads, the accuracy of the WCS-
based UQ algorithm is demonstrated first by a plot of error in estimation
metric vs time. This metric is computed at each time as:

$e_{\mu}=\frac{1}{n}\begin{Vmatrix}\frac{E(\textbf{x}_{t_{c}})-E(\textbf{x}_{t_{f}})}{E(\textbf{x}_{t_{f}})}\end{Vmatrix}_{2}$
(24)

where, $E(\textbf{x}_{t_{c}})$ is the mean of the state variable
$\textbf{x}_{t}\in\mathbb{R}^{n}$ obtained from the Equation 20.
$E(\textbf{x}_{t_{f}})$ is obtained by running the full model (Equation 15).
Figure 10 shows the plot of $e_{\mu}$ vs time for a span of one year.

Figure 10: Plot of error in estimation vs time

(a) Zone 1

(b) Zone 47

(c) Zone 80

(d) Zone 97

Figure 11: Adjacency Information for (a) Zone 1, (b) Zone 47, (c) Zone 80 and
(d) Zone 97

Figure 10 shows the accuracy of the WCS-based UQ approach. The identified WCSs
ignore the inter-zone mixing of the surface temperatures. This mixing has very
little effect on the results. The error metric $e_{\mu}$ converges to a value
between 0.02 and 0.20 (2-20$\%$) towards the later part of the year. The slow
convergence might be due to the presence of large number of state variables
compared to fewer available observations. Also, the individual zonal
temperatures are not an explicit function of the surface temperatures. This
causes difficulty in the estimation of the surface temperatures $T_{o}$ and
$T_{i}$ from $T_{k}$. Due to this convergence rate, subsequent analysis has
been shown skipping the first month of the year.

Figure 11 shows the accuracy in estimating the zonal temperatures for the four
zones; Zone 1, Zone 47, Zone 80 and Zone 97.

Figures 12, 13, and 14 show the estimated non-envelope loads for the four
specific zones for the month of February, June, and October respectively. The
loads follow the trend of occupancy, as expected, and have daily peaks midday
and drop at night. The thermal load values typically increase daily during
weekdays, before decreasing over the weekends. This is likely due to the high
thermal mass associated with the zone, particularly the floor and the exterior
walls, which allows the zone to retain heat at night, and provides an
increased baseline non-envelope load each day. Zone 1, for example, has an
especially high thermal capacitance due to the surrounding soil, which
minimizes the temperature drop in the evenings.

Because the model has approximately 2700 state variables, and only 132
outputs, the resulting estimations are equally distributed based on the
coefficients of the state space coefficients, and the estimated results follow
a pattern similar to the estimated zonal temperatures. The assumption that the
thermal loads are slow-changing with a $\dot{T}_{k}^{int}=0$ minimizes the
estimated nightly reduction in non-envelope load, since it does not permit a
drastic change. Note that these values include the uncertainty due to the
baseboard heat.

(a) Zone 1

(b) Zone 47

(c) Zone 80

(d) Zone 97

Figure 12: Non-Envelope Loads for (a) Zone 1, (b) Zone 47, (c) Zone 80 and (d)
Zone 97 for the month of February

(a) Zone 1

(b) Zone 47

(c) Zone 80

(d) Zone 97

Figure 13: Non-Envelope Loads for (a) Zone 1, (b) Zone 47, (c) Zone 80 and (d)
Zone 97 for the month of June

(a) Zone 1

(b) Zone 47

(c) Zone 80

(d) Zone 97

Figure 14: Non-Envelope Loads for (a) Zone 1, (b) Zone 47, (c) Zone 80 and (d)
Zone 97 for the month of October

The solar gains for some of the surfaces for Zone 47 and Zone 97 are displayed
in Figure 15, 16, and 17 for the three months (February, June, and October).
The graph shows an equal distribution of the solar load into the interior
surfaces. This is an expected result, as eQuest models the radiation as
uniformly distributed [43]; actual observed results would be likely to differ
somewhat [15] due to the differences in direct and reflected solar radiation
on the surfaces. Note that the surface heat flux somewhat follows the trend of
occupancy; this is likely due to the radiant heat output of the lighting in
the zone. This is especially apparent at night, when the solar gain is
expected to drop substantially, but some lighting remains on. As with the
thermal loads, the solar gains follow a similar pattern to the estimated
outputs due to the high number of estimated state variables.

(a) Zone 47

(b) Zone 97

Figure 15: Solar Gains for (a) Zone 47 and (b) Zone 97 for the month of
February

(a) Zone 47

(b) Zone 97

Figure 16: Solar Gains for (a) Zone 47 and (b) Zone 97 for the month of June

(a) Zone 47

(b) Zone 97

Figure 17: Solar Gains for (a) Zone 47 and (b) Zone 97 for the month of June

##  6 Conclusion and Future Works

In this work, a reduced-order thermal BEM has been developed for a large scale
office/school building. The lumped capacitance RC network model is used to
calculate the non-envelope load and solar gain for interior surfaces. Using
the simplified BEM in conjunction with the Weakly Connected Subsystems optimal
estimation method allows for estimating these parameters for each thermal zone
along with associated uncertainties.

One significant advantage of the methodology is the ability to reduce the
computation expense of large-scale dynamical systems while maintaining
accuracy and providing uncertainty information at each step. Thus, the BEM
demonstrated in this work can be built upon, adding detail to offer more
precision. For instance, we included infiltration as part of the non-envelope
load estimator; instead, this component could be broken out separately to
explore the impact of the envelope leakage on the building model. In the same
vein, the BEM can be used to compare expected building behavior against the
actual building states. This in turn could detect anomalies in the building
operation. Used in conjunction with Model Predictive Control, one can also
optimize energy usage and occupant comfort.

The flexibility offered with this technique allows us to vary parameters over
time, which is not available in some of the more comprehensive simulation
programs. For example, foam insulation is known to degrade over time, due to
the diffusion of thermal gases as the material ages [57]. Most modeling
programs, such as eQUEST, require that the thermal conductivity of material
remain static over the duration of the simulation. The conductivity at each
time step can be used to provide a more precise representation of the dynamic
BEM using the lumped capacitance model.

The technique discussed in this work can be applied to any building and is
especially suited for large buildings with diverse occupancy, due to the
inter-zone effects, and variable internal loads. With the noise filtered out
and the zonal temperatures being estimated, the solar gains and the non-
envelope loads can be determined.

## Acknowledgement

This material is based upon work supported by the National Science Foundation
under Grant NSF CMMI #1301235. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the author(s) and do
not necessarily reflect the views of the National Science Foundation.

## References

  * [1] U. E. I. Administration, “Energy consumption estimates by sector,” 2017.
  * [2] U. E. I. Administration, “Emissions of greenhouse gases in the united states 2009,” 2011.
  * [3] N. Y. State, “Directing state agencies and authorities to improve the energy efficiency of state buildings,” 2010.
  * [4] N. Y. S. D. of Public Service, “System benefits charge,” 2011.
  * [5] U.S Department of Energy, “Energyplus,” 2014.
  * [6] J. J. Hirsch et al., “equest, the quick energy simulation tool,” DOE2. com, 2006.
  * [7] TRANE, “Trace 700,” 2017.
  * [8] Carrier, “Hap,” 2017.
  * [9] Y. Ding, Y. Shen, J. Wang, and X. Shi, “Uncertainty sources and calculation approaches for building energy simulation models,” Energy Procedia, vol. 78, pp. 2566–2571, 2015.
  * [10] Y. Sun, H. Su, C. J. Wu, and G. Augenbroe, “Quantification of model form uncertainty in the calculation of solar diffuse irradiation on inclined surfaces for building energy simulation,” Journal of Building Performance Simulation, vol. 8, no. 4, pp. 253–265, 2015.
  * [11] A. Rysanek and R. Choudhary, “Optimum building energy retrofits under technical and economic uncertainty,” Energy and Buildings, vol. 57, pp. 324–337, 2013.
  * [12] B. Eisenhower, Z. O’Neill, V. A. Fonoberov, and I. Mezić, “Uncertainty and sensitivity decomposition of building energy models,” Journal of Building Performance Simulation, vol. 5, no. 3, pp. 171–184, 2012.
  * [13] A. Mukherjee, R. Rai, P. Singla, T. Singh, and A. K. Patra, “Comparison of linearization and graph clustering methods for uncertainty quantification of large scale dynamical systems,” International Journal for Uncertainty Quantification, vol. 7, no. 1, pp. 23–56, 2017.
  * [14] J. Vivian, A. Zarrella, G. Emmi, and M. De Carli, “An evaluation of the suitability of lumped-capacitance models in calculating energy needs and thermal behaviour of buildings,” Energy and Buildings, 2017.
  * [15] L. He, B. Lei, H. Bi, and T. Yu, “Simplified building thermal model used for optimal control of radiant cooling system,” Mathematical Problems in Engineering, vol. 2016, 2016.
  * [16] D. W. U. Perera, C. Pfeiffer, and N.-O. Skeie, “Modeling and simulation of a multi-zone building for better control,” 2014.
  * [17] A. Mukherjee, R. Rai, P. Singla, T. Singh, and A. Patra, “Laplacian graph based approach for uncertainty quantification of large scale dynamical systems,” in American Control Conference (ACC), 2015, pp. 3998–4003, IEEE, 2015.
  * [18] A. Mukherjee, R. Rai, P. Singla, T. Singh, and A. Patra, “Non-negative matrix factorization based uncertainty quantification method for complex networked systems,” in ASME 2015 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, pp. V02AT03A048–V02AT03A048, American Society of Mechanical Engineers, 2015.
  * [19] A. P. Ramallo-González, M. E. Eames, and D. A. Coley, “Lumped parameter models for building thermal modelling: An analytic approach to simplifying complex multi-layered constructions,” Energy and Buildings, vol. 60, pp. 174–184, 2013.
  * [20] J. Širokỳ, F. Oldewurtel, J. Cigler, and S. Prívara, “Experimental analysis of model predictive control for an energy efficient building heating system,” Applied energy, vol. 88, no. 9, pp. 3079–3087, 2011.
  * [21] A. Li, Y. Sun, and X. Xu, “Development of a simplified resistance and capacitance (rc)-network model for pipe-embedded concrete radiant floors,” Energy and Buildings, 2017.
  * [22] E. Á. R. Jara, F. J. S. de la Flor, S. Á. Domínguez, J. L. M. Félix, and J. M. S. Lissén, “A new analytical approach for simplified thermal modelling of buildings: Self-adjusting rc-network model,” Energy and Buildings, vol. 130, pp. 85–97, 2016.
  * [23] C. Underwood, “An improved lumped parameter method for building thermal modelling,” Energy and Buildings, vol. 79, pp. 191–201, 2014.
  * [24] K. J. Kircher and K. M. Zhang, “On the lumped capacitance approximation accuracy in rc network building models,” Energy and Buildings, vol. 108, pp. 454–462, 2015.
  * [25] S. Wang and X. Xu, “Parameter estimation of internal thermal mass of building dynamic models using genetic algorithm,” Energy conversion and management, vol. 47, no. 13, pp. 1927–1941, 2006.
  * [26] S. Goyal, C. Liao, and P. Barooah, “Identification of multi-zone building thermal interaction model from data,” in Decision and Control and European Control Conference (CDC-ECC), 2011 50th IEEE Conference on, pp. 181–186, IEEE, 2011.
  * [27] Z. O’Neill, S. Narayanan, and R. Brahme, “Model-based thermal load estimation in buildings,” IBPSA-USA Journal, vol. 4, no. 1, pp. 474–481, 2010.
  * [28] M. Woloszyn and I. Beausoleil-Morrison, “Treating uncertainty in building performance simulation,” 2017.
  * [29] A.-T. Nguyen, S. Reiter, and P. Rigo, “A review on simulation-based optimization methods applied to building performance analysis,” Applied Energy, vol. 113, pp. 1043–1058, 2014.
  * [30] M. Kavgic, A. Summerfield, D. Mumovic, and Z. Stevanovic, “Application of a monte carlo model to predict space heating energy use of belgrade’s housing stock,” Journal of Building Performance Simulation, vol. 8, no. 6, pp. 375–390, 2015.
  * [31] R. Baetens and D. Saelens, “Modelling uncertainty in district energy simulations by stochastic residential occupant behaviour,” Journal of Building Performance Simulation, vol. 9, no. 4, pp. 431–447, 2016.
  * [32] G. C. Rodríguez, A. C. Andrés, F. D. Muñoz, J. M. C. López, and Y. Zhang, “Uncertainties and sensitivity analysis in building energy simulation using macroparameters,” Energy and Buildings, vol. 67, pp. 79–87, 2013.
  * [33] W. Tian, “A review of sensitivity analysis methods in building energy analysis,” Renewable and Sustainable Energy Reviews, vol. 20, pp. 411–419, 2013.
  * [34] F. Oldewurtel, A. Parisio, C. N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, and M. Morari, “Use of model predictive control and weather forecasts for energy efficient building climate control,” Energy and Buildings, vol. 45, pp. 15–27, 2012.
  * [35] F. Oldewurtel, A. Parisio, C. N. Jones, M. Morari, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, and K. Wirth, “Energy efficient building climate control using stochastic model predictive control and weather predictions,” in American control conference (ACC), 2010, pp. 5100–5105, IEEE, 2010.
  * [36] F. Oldewurtel, C. N. Jones, A. Parisio, and M. Morari, “Stochastic model predictive control for building climate control,” IEEE Transactions on Control Systems Technology, vol. 22, no. 3, pp. 1198–1205, 2014.
  * [37] S. Privara, J. Širokỳ, L. Ferkl, and J. Cigler, “Model predictive control of a building heating system: The first experience,” Energy and Buildings, vol. 43, no. 2, pp. 564–572, 2011.
  * [38] S. H. Kim, “Building demand-side control using thermal energy storage under uncertainty: An adaptive multiple model-based predictive control (mmpc) approach,” Building and Environment, vol. 67, pp. 111–128, 2013.
  * [39] M. Maasoumy, M. Razmara, M. Shahbakhti, and A. S. Vincentelli, “Handling model uncertainty in model predictive control for energy efficient buildings,” Energy and Buildings, vol. 77, pp. 377–392, 2014.
  * [40] A. H. Stroud, “Approximate calculation of multiple integrals,” 1971.
  * [41] N. Oceanic and A. Administration, “National weather service,” 2018.
  * [42] ASHRAE(2013), 2013 ASHRAE Handbook: Fundamentals.  ASHRAE, 2013.
  * [43] U. DoE, “Energyplus engineering reference,” EnergyPlus Version 8.7 Documentation, 2016.
  * [44] D. E. Fisher and C. O. Pedersen, “Convective heat transfer in building energy and thermal load calculations,” tech. rep., American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc., Atlanta, GA (United States), 1997.
  * [45] F. Winklemann, “Underground surfaces: How to get better underground surface heat transfer calculations in doe2.1e,” 2003.
  * [46] ASHRAE, ANSI/ASHRAE/IES Standard 90.1-2013. Energy Standard for Buildings Except Low-Rise Residential Buildings.  ASHRAE, 2013.
  * [47] W. Marion and K. Urban, User’s Manual for TMY2s: Typical Meteorological Years: Derived from the 1961-1990 National Solar Radiation Data Base.  National Renewable Energy Laboratory, 1995.
  * [48] R. Perez, P. Ineichen, R. Seals, J. Michalsky, and R. Stewart, “Modeling daylight availability and irradiance components from direct and global irradiance,” Solar energy, vol. 44, no. 5, pp. 271–289, 1990.
  * [49] R. E. Kalman et al., “A new approach to linear filtering and prediction problems,” Journal of basic Engineering, vol. 82, no. 1, pp. 35–45, 1960\.
  * [50] A. H. Jazwinski, Stochastic processes and filtering theory.  Courier Corporation, 2007.
  * [51] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, “Fast unfolding of communities in large networks,” Journal of statistical mechanics: theory and experiment, vol. 2008, no. 10, p. P10008, 2008.
  * [52] ASHRAE, 90.1-2013 User’s Manual.  ASHRAE, 2014.
  * [53] J. J. Hirsch et al., “Doe 2.2 building energy use and cost analysis program,” See http://www. doe2. com, 2003.
  * [54] I. C. Council, “International mechanical code,” 2015.
  * [55] RESNET, “Commercial buildings energy modeling guidelines and procedures,” 2010\.
  * [56] I. Siemens Industry, “Hvac components and controls catalog,” 2013.
  * [57] P. de Wilde, W. Tian, and G. Augenbroe, “Longitudinal prediction of the operational energy use of buildings,” Building and Environment, vol. 46, no. 8, pp. 1670–1680, 2011.