WGMwgm QEqe EPep PMSpms BECbec DEde

[orcid=0000-0002-3295-412X] [1] Conceptualization, Methodology, Software,
Validation, Formal analysis, Data Curation, Visualization, Writing - Original
Draft

Conceptualization, Methodology, Writing - Review & Editing, Supervision

Writing - Review & Editing, Resources, Funding acquisition

Conceptualization, Methodology, Writing - Review & Editing, Supervision

[cor1]Corresponding author<EMAIL_ADDRESS>(L. Di Natale)

# Towards Scalable Physically Consistent Neural Networks: an Application to
Data-driven Multi-zone Thermal Building Models

Di Natale L  Svetozarevic B  Heer P  Jones C.N  Urban Energy Systems
Laboratory, Swiss Federal Laboratories for Materials Science and Technology
(Empa), 8600 Dübendorf, Switzerland  Laboratoire d’Automatique, EPFL, 1015
Lausanne, Switzerland

###### Abstract

With more and more data being collected, data-driven modeling methods have
been gaining in popularity in recent years. While physically sound, classical
gray-box models are often cumbersome to identify and scale, and their accuracy
might be hindered by their limited expressiveness. On the other hand,
classical black-box methods, typically relying on Neural Networks (NNs)
nowadays, often achieve impressive performance, even at scale, by deriving
statistical patterns from data. However, they remain completely oblivious to
the underlying physical laws, which may lead to potentially catastrophic
failures if decisions for real-world physical systems are based on them.
Physically Consistent Neural Networks (PCNNs) were recently developed to
address these aforementioned issues, ensuring physical consistency while still
leveraging NNs to attain state-of-the-art accuracy.

In this work, we scale PCNNs to model building temperature dynamics and
propose a thorough comparison with classical gray-box and black-box methods.
More precisely, we design three distinct PCNN extensions, thereby exemplifying
the modularity and flexibility of the architecture, and formally prove their
physical consistency. In the presented case study, PCNNs are shown to achieve
state-of-the-art accuracy, even outperforming classical NN-based models
despite their constrained structure. Our investigations furthermore provide a
clear illustration of NNs achieving seemingly good performance while remaining
completely physics-agnostic, which can be misleading in practice. While this
performance comes at the cost of computational complexity, PCNNs on the other
hand show accuracy improvements of $17$ – $35\%$ compared to all other
physically consistent methods, paving the way for scalable physically
consistent models with state-of-the-art performance.

###### keywords:

Neural Network
Physical consistency
Building modeling
Deep Learning
Physics-inspired

##  1 Introduction

Under the pressing issue of climate change, there is a worldwide effort to
decrease our global energy consumption. Being responsible for a large part of
the final energy consumption and greenhouse gas emissions [1], buildings are a
primary target in that trend, with space heating and cooling being the main
identified culprits [2]. Interestingly, we can intervene at any stage of the
life of a building to decrease its energy intensity, either designing and
constructing more efficient new structures [3], retrofitting old edifices [4],
or designing smart controllers to minimize the energy consumption of the
current building stock [5]. However, while decreasing the energy consumption
of buildings is the main goal of advanced control algorithms, this cannot be
done at the expense of the comfort of the inhabitants, who require the
temperature to stay within a comfortable range [6]. This calls for accurate
building temperature models to close the sim2real gap of advanced control
algorithms [7, 8]. Indeed, Model Predictive Control (MPC) uses a model to
predict the impact of possible power input sequences and choose the optimal
one [9], for example, and Reinforcement Learning (RL) control policies have to
be trained in simulation prior to their deployment [6, 10].

###  1.1 Towards data-driven methods

Since the evolution of the temperature in a thermal zone is governed by the
laws of thermodynamics, the most natural way to model it is to write down the
corresponding Ordinary Differential Equations (ODEs) and then use custom
solvers or discretization schemes to propagate them through time, such as in
[11, 12]. To alleviate the engineering burden of constructing the ODEs
describing building temperature dynamics, allow more complex structures to be
modeled, and accelerate the entire pipeline, custom modeling tools are often
used in practice, such as EnergyPlus, Modelica, TRNSYS, or IDA ICE [13, 14,
15]. Such detailed simulation tools however still require expert knowledge and
access to many design parameters that are often not directly available [16],
which makes them infamously hard to calibrate [17]. Moreover, solving the
complex underlying ODEs to simulate each time step can entail a significant
computational burden at runtime [18].

In recent years, owing to the growing amount of data collected in buildings,
researchers started to leverage data-driven methods, bypassing the cumbersome
procedures and expert knowledge required to set up classical physics-based
models [19]. This gave rise to so-called gray-box or black-box frameworks,
both of which use historical data for calibration or training purposes, as
pictured in Figure 1.

Figure 1: General pipeline of data-driven building thermal modeling
frameworks.

###  1.2 Gray- and black-box models

When a control-oriented thermal building model is designed, typically for MPC,
data is in most cases used to identify the parameters of a simplified physics-
based model [20], usually a low-order Resistance-Capacitance (RC) model, such
as in [21, 22, 16, 23, 24]. These models are particularly popular due to their
ease of implementation, interpretability, close ties to the underlying
physics, and because they often give rise to linear dynamics. The latter
characteristic is indeed particularly desirable in MPC applications since an
appropriate choice of objective function then renders the optimization problem
to solve at each time step convex and hence tractable. Nonetheless, the
parameter identification procedure of RC models is generally nontrivial and
sensitive to the data quality [25, 26], which partially explains why low-order
models often perform better than complex ones [27, 28]. Alternatively or
additionally, when data is available, the residual errors of an often
simplified model can be fit to improve its performance, such as in [29].

While the aforementioned gray-box methods only require limited domain
knowledge since simplified ODEs are used, which ODEs to choose is not always
clear [30] and wrong choices might hinder the performance of the identified
model. To avoid this pitfall, when enough data is available, fully data-driven
black-box models might be used, such as AutoRegressive models with eXogenous
inputs (ARX) or Neural Networks (NNs). Indeed, they do not rely on any prior
knowledge of the system to model since they directly extract patterns from
data to explain and predict the behavior of the system. Consequently, black-
box methods are often easier and faster to deploy, more flexible, and thus
often more scalable than their gray-box counterparts [31, 32]. Furthermore,
since they do not have to follow a predefined underlying architecture, black-
box methods are generally more expressive, being capable of capturing unknown
nonlinear dynamics, and hence usually perform better [33].

###  1.3 Neural Networks and their inconsistencies

With the recent advances in Deep Learning (DL), NN-based solutions are gaining
in popularity to represent unknown and potentially highly nonlinear dynamics,
making them state-of-the-art solutions for time series modeling [34].
Unsurprisingly, given the broad range of applications of NNs, researchers have
already applied them to model building thermal dynamics, e.g. in [35, 5, 36].
When applying NNs to physical systems, one should however keep their well-
known generalization issue in mind [37], which is partially caused by the
underspecification plaguing DL applications [38]. This can indeed be
particularly problematic for NNs modeling physical processes since they are
physics-agnostic and might hence find unrealistic solutions [4]. If the
training data set does not span all the operating conditions of the system,
there is thus a risk for NNs to fail to generalize in a meaningful manner to
new conditions, an issue that is typically expected for thermal building
models since the collected data sets are generally inherently incomplete [39].

Deep NNs are indeed able to learn shortcuts [40], which means they might fit
the training data well without fundamentally understanding the problem, hence
failing to generalize. They for example attain superhuman performance on image
recognition tasks, and yet fail when undistinguishable noise is added [37] or
the background changes [41, 42]. They are also able to generate captions
without ever looking at the corresponding images [43], or detect pneumonia
from X-ray scans only by looking at hospital-specific tokens and each
hospital’s pneumonia prevalence [44].

While these are only a few examples, they clearly indicate how NNs can find
ways to perform extremely well without fundamentally solving the task at hand.
These flawed models are however unable to generalize and cannot be deployed in
real-world applications since we have no means to know how they will react to
new conditions. While tasks such as object recognition and captioning might be
hard to characterize in general, the case of physical system modeling is
different. Indeed, we often know the underlying physical laws and can hence
impose constraints on NNs that help them understand the task at hand.

###  1.4 Physics-inspired Neural Networks

To incorporate some knowledge of the underlying physics in NN training
procedures and promote desired system properties, counterbalancing the
aforementioned brittleness of classical NNs, researchers recently proposed to
design Physics-informed or Phyiscs-inspired NNs (PiNNs) [45]. While many works
modify the loss function of NNs to steer the learning towards physically
meaningful solutions [46, 47], these schemes cannot provide any guarantee
about the final model respecting the desired constraints. More systematic
approaches hence directly alter the networks’ architecture to ensure the
underlying physical laws are followed by design, i.e. at all times, such as in
[48, 49, 50, 39]. Additionally, since the desired properties are hard-coded in
such models, the loss function does not need to be altered, which avoids
common pitfalls of classical PiNNs, such as the difficult trade-off between
the accuracy and the physical consistency of the model, which can also
increase the amount of data needed [51]. These specific NN architectures
ensuring some physical properties are philosophically related to the
celebrated convolutional [52], recurrent [53], or graph [54] NN architectures,
which were designed to capture spatial, timely, or neighboring relationships
in the input data, respectively.

Despite the recent popularity of the field, to the best of the authors’
knowledge, PiNNs were only applied to thermal building modeling in [55, 56,
57]. Gokhale et al. relied on the classical PiNN framework, augmenting the
loss function of their NNs and creating latent states to include some physical
intuition in otherwise standard networks [55], while Nagarathinam et al.
designed a specific PiNN architecture for building control [56]. On the other
hand, Drgoňa et al. used NNs to replace the matrices in linear models of
building dynamics, which allowed them to enforce the stability and
dissipativity of the learned system by constraining the eigenvalues of one of
the NNs [57]. However, these works cannot provide guarantees about the
physical consistency of their solutions in general beyond stability and
dissipativity.

On the contrary, the Physically Consistent Neural Networks (PCNNs) developed
in previous work were theoretically proven to always yield physically
consistent temperatures predictions, outperformed a classical gray-box model,
and attained an accuracy on par with pure black-box models, but were limited
to single-zone temperature modeling [39].

###  1.5 Contribution

In this work, we propose three different extensions of single-zone PCNNs [39]
to model entire buildings, exemplifying how one can use the modularity of
PCNNs to either expand the physics-inspired or black-box module. We then show
that they all retain the desired physical consistency, both theoretically and
with a numerical analysis, and investigate their performance on a three-zone
building case study. Through extensive comparisons with several classical and
state-of-the-art gray- and black-box methods, we demonstrate the ability of
multi-zone PCNNs to leverage NNs to be more expressive than physically
grounded gray-box methods and even outperform purely black-box NNs in terms of
accuracy. Our investigations also clearly illustrate the phenomenon of
shortcut learning on building temperature data, with NNs able to fit the data
very accurately despite being completely oblivious to the underlying physics,
which can be misleading in practice. Altogether, these experiments prove the
effectiveness of the proposed PCNNs as thermal building models, alleviating
the need for any engineering overhead, following the underlying physical laws,
and reaching state-of-the-art accuracy.

The remainder of this paper is structured as follows. Section 2 firstly
defines the notion of physical consistency and recalls the main principles
behind PCNNs before the proposed multi-zone architectures are introduced and
theoretically analyzed in Section 3. Section 4 then details the case study,
implementation considerations, and baseline models. Finally, the results are
analyzed in Section 5 and Section 6 concludes the paper.

##  2 Background

In this section, we recall a few prerequisites required to understand the
methods proposed in this work, clarifying some definitions, formally
introducing the notion of physical consistency, and recalling the design of
single-zone PCNNs.

###  2.1 Definitions

Two zones are said to be adjacent if they share at least one common wall in a
building, and the collection of zones adjacent to a given zone $z$ form its
neighborhood $\mathcal{N}(z)$. Note that we consider a zone to be included in
its neighborhood, i.e. $z\in\mathcal{N}(z)$. Similarly, a zone is connected
with the outside if it is composed of at least one external wall. To
generalize the notion of neighborhood, we define the $n$-hop neighborhood
$\mathcal{N}^{n}(z)$ as the set of zones that can be reached in $n$ steps from
zone $z$, moving to an adjacent zone at each step. Note that, by definition,
we have $\mathcal{N}^{1}(z)=\mathcal{N}(z)$, and $y\in\mathcal{N}^{n}(z)\iff
z\in\mathcal{N}^{n}(y)$.

Throughout this work, we assume the building to be connected, i.e. there is no
zone (or group of zones) isolated from the rest. This assumption is trivial in
practice as one can easily train several separate models if this condition is
not met.

###  2.2 Physical consistency

In this paper, we deem the temperature model of a building $\mathcal{B}$ with
$m$ thermal zones to be physically consistent if the following conditions are
met for each zone $z\in\mathcal{B}$:

$\displaystyle\frac{\partial T^{z}_{k+i}}{\partial u^{z}_{k+j}}$
$\displaystyle>0$ $\displaystyle\forall$ $\displaystyle 0\leq j<i,$ (1)
$\displaystyle\frac{\partial T^{z}_{k+i}}{\partial T^{out}_{k+j}}$
$\displaystyle>0$ $\displaystyle\forall$ $\displaystyle 0\leq j<i,$ (2)
$\displaystyle\frac{\partial T^{z}_{k+i}}{\partial T^{y}_{k+j}}$
$\displaystyle>0$ $\displaystyle\forall$ $\displaystyle 0\leq j<i,\ \forall
y\in\mathcal{N}^{i-j}(z),$ (3)

where $T^{z}_{k}$ is the temperature in zone $z$ at time step $k$, $u$ its
heating/cooling power input, and $T^{out}$ represents the outside temperature.

For example, (1) implies that applying more heating power $u^{z}_{k+j}$ at
time step $k+j$ leads to higher temperatures $T^{z}_{k+i}$ for all all
subsequent time steps $i>j$. In other words, heating a zone has the expected
and intuitive impact of increasing its temperature, following the laws of
thermodynamics. Similarly, higher ambient temperatures induce higher
temperatures inside, and higher temperatures in zone $y\in\mathcal{N}^{n}(z)$
lead to higher temperatures in zone $z$ after $n$ steps. Note that cooling
powers are defined to be negative in this work, hence inducing lower
temperatures when more cooling is applied, as expected.

######  Remark 1 (Generalization of the approach).

Note that the definition of physical consistency proposed in (1)-(3) can
easily be extended for applications where additional criteria need to be met
by the learned model, to enforce physically consistent temperature predictions
with respect to solar gains, for example. Interestingly, these conditions can
also be seamlessly adapted to other fields beyond building modeling where
simple physical rules can be encoded in a similar fashion. One can then
construct a PCNN architecture following the principles presented in Sections
2.3 and 3 to ensure the learned model respect these desired criteria.

###  2.3 Single-zone PCNNs

Figure 2: Visualization of PCNN predictions at each step, where the inputs are
processed in parallel by a black-box and a physics-based module and then
summed.

Conceptually, PCNNs are composed of a black-box and a physics-inspired module
running in parallel to compute the next output at each step, as depicted in
Figure 2. The former captures potentially complex nonlinearities while the
latter ensures that predefined rules are respected, which typically represent
physical laws and can be encoded by conditions similar to the ones proposed in
(1)-(3). In the case of PCNNs modeling the evolution of the temperature in a
single thermal zone $z$ while respecting the criteria in (1)-(3), they can be
mathematically described as follows [39]:

$\displaystyle D_{k+1}$ $\displaystyle=D_{k}+f(x_{k},D_{k}),$ (4)
$\displaystyle E_{k+1}$
$\displaystyle=E_{k}+a_{h}\max{\\{u_{k},0\\}}+a_{c}\min{\\{u_{k},0\\}}$
$\displaystyle\quad-b(T_{k}-T^{out}_{k})-\sum_{z^{\prime}\in\mathcal{N}(z)}{c_{z^{\prime}}(T_{k}-T^{z^{\prime}}_{k})},$
(5) $\displaystyle T_{k+1}$ $\displaystyle=D_{k+1}+E_{k+1},$ (6)
$\displaystyle D_{k}$ $\displaystyle=T(k),$ $\displaystyle E_{k}$
$\displaystyle=0,$

where $D\in\mathbb{R}$ represents the evolution of the black-box module based
on a freely parametrized function $f:\mathbb{R}^{d+1}\to\mathbb{R}$, typically
composed of NNs, and $E\in\mathbb{R}$ is the energy accumulator, i.e. the
physics-inspired module. The latter is influenced by the power inputs
$u\in\mathbb{R}$ and heat transfers to the outside and adjacent zones. On the
other hand, $x\in\mathbb{R}^{d}$ regroups all the inputs that are not treated
in $E$, typically time information and solar irradiation measurements.
Finally, the constants $a_{h}$, $a_{c}$, $b$, and
$\\{c_{z^{\prime}}\\}_{z^{\prime}\in\mathcal{N}(z)}$ capture the impact of
heating, cooling, and heat losses to the outside and the neighboring zones on
the modeled zone temperature, respectively. Applying (4)-(6) recursively over
the prediction horizon, starting from the measured temperature $T(k)$ at time
$k$, PCNNs can predict the evolution of the temperature while satisfying the
criteria in (1)-(3). We refer the reader to the original paper for additional
details [39].

One important key to the effectiveness and generality of PCNNs comes from the
fact that all the parameters $a_{h}$, $a_{c}$, $b$,
$\\{c_{z^{\prime}}\\}_{z^{\prime}\in\mathcal{N}(z)}$, and $f$ are learned
simultaneously using automatic BackPropagation Through Time (BPTT) [58]
(Section 4.2). As mentioned in Remark 1, the very generic structure of PCNNs
can also be applied to model complex phenomena beyond thermal modeling,
typically where only part of the physics is well understood. Indeed, it is
always possible to adapt the structure of the physics-inspired module, which
might also include nonlinearities, let the black-box module capture completely
unknown dynamics in parallel, and seamlessly learn everything simultaneously
in an end-to-end pipeline.

######  Remark 2 (Consistency with respect to initial conditions).

Since we want to model several zones in this paper, the condition (3) is
different from the one enforced in [39]. In particular, it includes the case
$y=z$, $j=0$, i.e. we want the derivative of any zone temperature with respect
to its initial temperature to be positive to avoid spurious effects. This was
not considered in the original paper on single-zone modeling and requires a
slight modification of (4), as discussed in the following section.

##  3 Modeling the thermal behavior of buildings

This section proposes three different extensions of the single-zone PCNN
architecture to model an entire building and then discusses their physical
consistency.

###  3.1 Extensions to multi-zone PCNNs

While the PCNN architecture described in Section 2.3 was shown to work well
for single-zone modeling, this work proposes three possible extensions of this
framework to simultaneously capture the evolution of the temperature in
several interconnected zones exchanging energy, i.e. in a whole building. The
only additional information required is the topology of the modeled building,
i.e. we assume to know which zones are adjacent and which have an external
wall, and then learn its thermal behavior from data without engineering
overhead.

######  Remark 3 (Topology).

If the topology is unknown, one can also assume each pair of zones to be
adjacent and every zone to have an external wall and then learn to put non-
existing connection parameters to zero from data.

####  3.1.1 Learning several single-zone PCNNs

The most natural and straightforward extension is to separately learn one PCNN
for each of the zones to model, as depicted in Figure 3. Since this method
involves duplicating the original structure for each zone and fitting them
independently, we will refer to the final model of the building as the X-PCNN
architecture. Mathematically, for a given zone $z$, we can write the
corresponding equations as follows:

$\displaystyle D^{z}_{k+1}$ $\displaystyle=D^{z}_{k}+f^{z}(x^{z}_{k}),$ (7)
$\displaystyle E^{z}_{k+1}$
$\displaystyle=E^{z}_{k}+a^{z}_{h}\max{\\{u^{z}_{k},0\\}}+a^{z}_{c}\min{\\{u^{z}_{k},0\\}}$
$\displaystyle\ \ \
-b(T^{z}_{k}-T^{out}_{k})-\sum_{y\in\mathcal{N}(z)}{\tilde{c}^{z}_{y}(T^{z}_{k}-T^{y}_{k})},$
(8) $\displaystyle T^{z}_{k+1}$ $\displaystyle=D^{z}_{k+1}+E^{z}_{k+1},$ (9)
$\displaystyle D^{z}_{k}$ $\displaystyle=T^{z}(k),$ $\displaystyle E^{z}_{k}$
$\displaystyle=0,$

where the superscript $z$ denotes zone-dependent information and all the
variables have the same meaning as in Section 2.3.

Figure 3: Predictions of the proposed X-PCNN architecture, where the
temperature of each of the $m$ zones is predicted separately.

To retain the physical consistency of this model, however, one needs to ensure
that $\tilde{c}^{z}_{y}=\tilde{c}^{y}_{z}$ for each pair of adjacent zones $y$
and $z$, so that the amount of energy flowing from $z$ to $y$ always equals
the amount of energy received by $y$ from $z$, and vice versa. Since each zone
is modeled and trained separately in this case, such a condition cannot be
imposed during the learning phase, and we thus rely on a heuristic to correct
the parameters and enforce this desired property a posteriori. Once the models
have been trained, for every pair of adjacent zones $z$ and $y$, we compute
the average value identified by both PCNNs and define:

$\displaystyle c^{zy}=c^{yz}=\frac{\tilde{c}^{z}_{y}+\tilde{c}^{y}_{z}}{2}.$
(10)

For every zone $z$, we then replace $\tilde{c}^{z}_{y}$ with $c^{zy}$ in (8)
for all $y\in\mathcal{N}(z)$.

######  Remark 4 (Independence of $f^{z}$ from $D^{z}$).

Note the difference between (7) and (4), with $D^{z}$ not appearing in $f^{z}$
in (7). Following Remark 2, this ensures that condition (3) is respected at
all times, as analyzed in Section 3.2 and Remark 6.

####  3.1.2 Sharing the physics-inspired module

Figure 4: Predictions of the proposed M-PCNN architecture, where the physics-
inspired module is shared but multiple black-box modules are learned, one for
each zone.

To avoid the hand-crafted correction (10), which might significantly impact
the parameters learned by each PCNN, one can fuse all the physics-inspired
modules together, again leveraging our prior knowledge of the underlying
physical laws. This gives rise to the so-called M-PCNN architecture, pictured
in Figure 4, where distinct black-box modules are assigned to each zone, but
the physics-inspired module is shared and outputs a vector
$\bm{E}\in\mathbb{R}^{m}$ containing the energy accumulated in each zone at
each step:

$\displaystyle\bm{E}_{k+1}$
$\displaystyle=\bm{E}_{k}+\bm{a}_{h}\odot\max{\\{\bm{u}_{k},\bm{0}\\}}$
$\displaystyle\quad+\bm{a}_{c}\odot\min{\\{\bm{u}_{k},\bm{0}\\}}$ (11)
$\displaystyle\quad-\bm{b}\odot(\bm{T}_{k}-T^{out}_{k})-\bm{\Delta T}_{k},$
$\displaystyle\bm{E}_{k}$ $\displaystyle=\bm{0},$

where the bold notations correspond to vectorized quantities in
$\mathbb{R}^{m}$, one dimension for each zone, i.e.
$\bm{a_{h}}=[a_{h}^{1},...,a_{h}^{m}]^{T}$, and similarly for $\bm{a_{c}}$,
$\bm{b}$, and $\bm{u_{k}}$, and $\odot$ stands for the element-wise product of
two vectors. Finally, $\bm{\Delta T}_{k}\in\mathbb{R}^{m}$ corresponds to
energy transfer between each zone and its neighborhood:

$\displaystyle\bm{\Delta T}_{k}^{z}$
$\displaystyle=\sum_{y\in\mathcal{N}(z)}{c^{zy}(T^{z}_{k}-T^{y}_{k})},\qquad\forall
z\in\mathcal{B}$ (12)

where the subscript $z$ denotes the $z$-th entry of a vector. By definition,
we know $c^{zy}=c^{yz}$ if $y$ and $z$ are adjacent since both represent the
same heat transfer coefficient, which is easily enforced during training since
all the zones are now modeled simultaneously, avoiding the a posteriori
correction (10) required for X-PCNNs.

Each dimension of $\bm{T}\in\mathbb{R}^{m}$, i.e. the temperature in each zone
$z$, is then computed as the sum of the physics-inspired and black-box
modules, as before:

$\displaystyle\bm{T}_{k+1}^{z}$ $\displaystyle=D^{z}_{k+1}+\bm{E}_{k+1}^{z},$
(13) $\displaystyle D^{z}_{k+1}$ $\displaystyle=D^{z}_{k}+f^{z}(x^{z}_{k}),$
(14) $\displaystyle D^{z}_{k}$ $\displaystyle=T^{z}(k).$

####  3.1.3 Sharing both modules

To reduce the computational complexity of the model and introduce parameter
sharing between the zones – which are typically similar in the same building
–, we propose a third architecture, dubbed S-PCNN, where both the black-box
and physics-inspired modules are shared. In practice, this means that the
black-box module, typically consisting of NNs, now has $m$ outputs
corresponding to the main dynamics of each of the zones, as pictured in Figure
5. Using the vectorized notations as before, i.e. $\bm{D}\in\mathbb{R}^{m}$,
we can write the equations of this architecture as follows:

$\displaystyle\bm{D}_{k+1}$
$\displaystyle=\bm{D}_{k}+\bm{\tilde{f}}(\bm{\tilde{x}}_{k}),$ (15)
$\displaystyle\bm{E}_{k+1}$
$\displaystyle=\bm{E}_{k}+\bm{a}_{h}\odot\max{\\{\bm{u}_{k},\bm{0}\\}}$
$\displaystyle\quad+\bm{a}_{c}\odot\min{\\{\bm{u}_{k},\bm{0}\\}}$ (16)
$\displaystyle\quad-\bm{b}\odot(\bm{T}_{k}-T^{out}_{k})-\bm{\Delta T}_{k},$
$\displaystyle\bm{T}_{k+1}$ $\displaystyle=\bm{D}_{k+1}+\bm{E}_{k+1},$ (17)
$\displaystyle\bm{D}_{k}$ $\displaystyle=\bm{T}(k),$ $\displaystyle\bm{E}_{k}$
$\displaystyle=\bm{0},$

where the physics-inspired module is the same as for the M-PCNN but we now
only have one shared nonlinear function
$\bm{\tilde{f}}:\mathbb{R}^{d^{\prime}}\to\mathbb{R}^{m}$ transforming the
inputs $\bm{\tilde{x}}\in\mathbb{R}^{d^{\prime}}$. Throughout this work, we
only consider external inputs that are shared by all the zones, i.e.
$d^{\prime}=d$ and $\bm{\tilde{x}}:=x^{z},\ \forall z\in\mathcal{B}$.

Figure 5: Predictions of the proposed S-PCNN architecture, sharing both the
physics-inspired and black-box modules.

######  Remark 5 (Zone-dependent inputs).

If some measurements differ zone by zone, one can either stack them in a
vector $\bm{\tilde{x}}=[(x^{1})^{\top},...,(x^{m})^{\top}]^{\top}$ and use
(15) as is or for example design a shared function
$\bm{\tilde{f}}:\mathbb{R}^{d}\to\mathbb{R}$ and modify (15) to
$\bm{D}_{k+1}^{z}=\bm{D}_{k}^{z}+\bm{\tilde{f}}(x^{z}_{k})$, $\forall
z\in\mathcal{B}$.

###  3.2 Thermodynamical consistency

Relying on the transformation (10) ensuring that heat transfer coefficients
between each adjacent zones are equal in the corresponding single-zone PCNNs,
one can vectorize the X-PCNN physics-inspired module (8), putting the
parameters of each zone in vectors, to get:

$\displaystyle\bm{E}_{k+1}$
$\displaystyle=\bm{E}_{k}+\bm{a}_{h}\odot\max{\\{\bm{u}_{k},\bm{0}\\}}$
$\displaystyle\quad+\bm{a}_{c}\odot\min{\\{\bm{u}_{k},\bm{0}\\}}$ (18)
$\displaystyle\quad-\bm{b}\odot(\bm{T}_{k}-T^{out}_{k})-\bm{\Delta T}_{k},$

using the definition of $\bm{\Delta T}$ from (12). As can be seen directly,
this expression is the same as the ones describing physics-inspired modules of
the M-PCNN and S-PCNN architectures in (11) and (16). This means all the
proposed multi-zone PCNNs rely on the same physical model at inference time,
with however possibly different parameter values learned during training. This
is intuitively expected since they all model the same thermal effects and
hence have to follow the same physical principles.

Similarly, we can rewrite the black-box modules of the X-PCNN and M-PCNN
architectures in vectorized form as:

$\displaystyle\bm{D}_{k+1}$
$\displaystyle=\bm{D}_{k}+\bm{\bar{f}}(\bm{\bar{x}}_{k}),$ (19)

where $\bm{\bar{f}}=[f^{1}(x^{1}),...,f^{m}(x^{m})]^{T}$ and
$\bm{\bar{x}}=[(x^{1})^{\top},...,(x^{m})^{\top}]^{\top}$ groups the different
inputs.

For the three proposed architectures, putting (18) and (19) together, we hence
get:

$\displaystyle\bm{T}_{k+1}$ $\displaystyle=\bm{D}_{k+1}+\bm{E}_{k+1}$
$\displaystyle=\bm{T}_{k}+\bm{f}(\bm{x}_{k})+\bm{a}_{h}\odot\max{\\{\bm{u}_{k},\bm{0}\\}}$
$\displaystyle\quad+\bm{a}_{c}\odot\min{\\{\bm{u}_{k},\bm{0}\\}}$ (20)
$\displaystyle\quad-\bm{b}\odot(\bm{T}_{k}-T^{out}_{k})-\bm{\Delta T}_{k},$
$\displaystyle\bm{D}_{0}$ $\displaystyle=\bm{T}(k),$ $\displaystyle\bm{E}_{0}$
$\displaystyle=\bm{0},$

where $\bm{f}(\bm{x}_{k})$ stands for $\bm{\tilde{f}}(\bm{\tilde{x}}_{k})$ or
$\bm{\bar{f}}(\bm{\bar{x}}_{k})$ for S-PCNNs, respectively X- and M-PCNNs. The
only structural difference between the three proposed models (once the heat
transfer coefficients of the X-PCNN have been adjusted) hence comes from the
form of $\bm{f}(\bm{x})$. Remarkably, however, this does not impact their
physical consistency, as demonstrated in the following two propositions.

######  Proposition 1 (Heat propagation).

Independently of the structure of $\bm{f}$ and $\bm{x}$, any model of the form
(20) satisfies:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j}^{y}}$
$\displaystyle\geq 0\qquad\forall z,y\in\mathcal{B},\ \forall 0\leq j<i,$ (21)

with equality if and only if $y\notin\mathcal{N}^{(i-j)}(z)$, as long the
following conditions hold:

$\displaystyle b^{z}+\sum_{y\in\mathcal{N}(z)}c^{zy}$ $\displaystyle<1,$
$\displaystyle\forall z$ $\displaystyle\in\mathcal{B},$ (22) $\displaystyle
c^{zy}$ $\displaystyle>0,$ $\displaystyle\forall z$
$\displaystyle\in\mathcal{B},\ \forall y\in\mathcal{N}(z).$ (23)

###### Proof.

See Appendix A.1. ∎

In words, Proposition 1 means that heat propagates from any zone $y$ to all
the other zones $z$ as physically expected, i.e. higher temperatures in a
given zone $y$ will eventually lead to higher temperatures in all the other
zones, i.e. after $(i-j)$ steps.

Proposition 1 can then be used to prove the following proposition stating that
heating or cooling any zone ultimately increases, respectively decreases, the
temperature in the whole building through heat transfers and that higher and
lower ambient temperatures also impact the building as expected.

######  Proposition 2 (Physical consistency with respect to inputs).

Independently of the structure of $\bm{f}$ and $\bm{x}$, any model of the form
(20) satisfies:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{u}_{k+j}^{y}}$
$\displaystyle\geq 0,$ $\displaystyle\forall$ $\displaystyle
z,y\in\mathcal{B},\ \forall 0\leq j<i,$ (24)

with equality if and only if $y\notin\mathcal{N}^{(i-j-1)}(z)$, and

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial T^{out}_{k+j}}$
$\displaystyle>0,$ $\displaystyle\forall$ $\displaystyle z\in\mathcal{B},\
\forall 0\leq j<i,$ (25)

as long as (22)-(23) hold and:

$\displaystyle a^{z}_{h},a^{z}_{c},b^{z}$ $\displaystyle>0,$
$\displaystyle\forall$ $\displaystyle z\in\mathcal{B}.$ (26)

###### Proof.

See Appendix A.2. ∎

######  Corollary 1 (Physical consistency of PCNNs).

Independently of the structure of $\bm{f}$ and $\bm{x}$, any model of the form
(20) respects the physical consistency criteria (1)-(3) if:

$\displaystyle b^{z}+\sum_{y\in\mathcal{N}(z)}c^{zy}$ $\displaystyle<1,$
$\displaystyle\forall z$ $\displaystyle\in\mathcal{B},$ (27) $\displaystyle
a^{z}_{h},a^{z}_{c},b^{z},c^{zy}$ $\displaystyle>0,$ $\displaystyle\forall z$
$\displaystyle\in\mathcal{B},\ \forall y\in\mathcal{N}(z).$ (28)

###### Proof.

Assuming that (27) and (28) hold, we can apply Propositions 1 and 2. Setting
$z=y$ in (24) and recalling that any zone is in its own neighborhood – which
implies positiveness of (24) –, PCNNs satisfy (1). The satisfaction of (2)
directly follows from the second part of Proposition 2. Finally, according to
Proposition 1, (21) is strictly positive if and only if
$y\in\mathcal{N}^{(i-j)}(z)$, satisfying (3). ∎

This corollary thus proves that each of the proposed PCNN architectures
remains physically consistent as long as all the parameters $\bm{a_{h}}$,
$\bm{a_{c}}$, $\bm{b}$, and $\bm{c}$ are small positive constants. Note that
this makes intuitive sense since all these parameters correspond to inverses
of resistances and capacitances, hence small positive numbers, in real
buildings. Interestingly, these conditions can easily be enforced during the
training procedure without modifying the classical BPTT algorithm, hence
allowing us to rely on well-developed tools to train our models, as detailed
in Section 4.2.

Remarkably, the strength of our approach lies in the fact that all the models
will remain consistent whatever the structure of $\bm{f}$ is, being shared or
not111This is the main difference between the M-PCNN and S-PCNN architectures
in our case, for example., composed of NNs or other nonlinearities. This gives
the user complete freedom in the design of the black-box module without
jeopardizing the physical consistency of the model. Similarly, all the
parameters of the physics-inspired module might for example be time-varying or
computed as nonlinear functions of external inputs without impacting the
consistency of the model as long as they stay small and positive at all times.

######  Remark 6 (Inputs of $\bm{f}$).

While the structure of $\bm{f}$ does not impact the validity of Propositions 1
and 2, its inputs do. In particular, $\bm{f}$ has to be independent of
$\\{T,u,T^{out}\\}$ for the first step of the proofs of both propositions to
hold in general (Appendix A). If $\bm{f}=\bm{f}(\bm{x},\bm{D})$ for example,
it would modify the case $z=y$ in equation (40) in Appendix A.1, and the
satisfaction of (22) would then not be sufficient to guarantee the required
nonnegativity of the partial derivatives in (21).

######  Remark 7 (A control perspective).

The multi-zone PCNNs (20) are power input-affine. This makes such models
interesting in control applications, typically for MPC schemes aimed at
decreasing the energy consumption of buildings.

##  4 Case study

To assess the quality of the multi-zone PCNN architectures detailed in Section
3, we carry out an extensive performance analysis on a case study, where the
objective is to predict the temperature dynamics over three-day-long horizons
with $15\text{\,}\mathrm{min}$ time steps. This section presents the building
where the data was collected, implementation details, and then introduces the
other gray- and black-box methods used as benchmarks.

###  4.1 Dataset

Figure 6: NEST, Duebendorf, Switzerland, with UMAR highlighted in white ©
Zooey Braun, Stuttgart.

The data used in this work was collected in NEST, a vertically integrated
district composed of several units [59], located in Duebendorf, Switzerland,
and pictured in Figure 6. In this paper, we focus on the Urban Mining and
Recycling (UMAR) unit, circled in white in Figure 6, an apartment with two
bedrooms and a living room in between them. We are thus modeling three thermal
zones arranged in a line in this work, i.e. Zone $1$ is connected to Zone $2$,
and Zone $2$ is also connected to Zone $3$, and each of them has at least one
external wall. In heating mode, all the zones are heated by letting hot water
flow through ceiling panels. During the cooling season, on the other hand,
cold water can flow through the panels to cool down the zones.

We rely on three years of data collected between May 2019 and May 2022 and
preprocessed as explained in [39, App. C]. This involved downsampling the data
to $15\text{\,}\mathrm{min}$ intervals, smoothing the time series, and
disaggregating the thermal power consumption of UMAR into the consumption of
each zone. Besides these computed thermal power inputs, the dataset also
contains measurements of the temperature in each zone and outside, the
horizontal solar irradiation on-site, and the status of the system, i.e. if it
is in heating or cooling mode. To facilitate the learning process of the NNs
in the black-box modules, we completed the data with additional time
information, i.e. the day of the week and the sine and cosine transformations
of the time of the day and the month [39, App. D]. Finally, we split the data
in a training and a validation set, respectively denoted $\mathcal{D}_{t}$ and
$\mathcal{D}_{v}$, containing possibly overlapping time series of up to three
days of data, as detailed in Appendix B.

###  4.2 Implementation details

To ensure the physical consistency of the proposed multi-zone PCNNs, i.e. to
fulfill the conditions (27) and (28), we parametrize the log-value of each
parameter, i.e. we learn
$\tilde{a}^{z}_{h},\tilde{a}^{z}_{c},\tilde{b}^{z},\tilde{c}^{zy}$, $\forall
y\in\mathcal{B}$, $\forall y\in\mathcal{N}(z)$ and define:

$\displaystyle s$ $\displaystyle=s_{0}\exp{(\tilde{s})},$
$\displaystyle\forall s=\\{a^{z}_{h},a^{z}_{c},b^{z},c^{zy}\\}$ (29)

where $s_{0}$ is the initial value of the parameter, defined using the same
rules of thumb as in [39]. Starting from $\tilde{s}=0$, PCNNs hence learn to
scale the initial value $s_{0}$ instead of modifying it directly, which is
more numerically stable and ensures that $s$ stays small enough, while the
exponential function keeps all the parameters positive at all times.

######  Remark 8 (Upperbound on $s$).

While $\tilde{b}^{z}$ or $\tilde{c}^{zy}$ can in principle grow uncontrollably
and lead to a violation of the necessary condition (27), this was not an issue
in our experiments. Nonetheless, one can always introduce bounds on the
learned values $\tilde{s}$ if required, typically leveraging activation
functions like the sigmoid or hyperbolic tangent to control the range of
learned values.

In this paper, each function $\bm{f}$ has the same encoder-LSTM-decoder
architecture, which is repeated when several modules are required for X-PCNNs
and M-PCNNs. Both the encoder and decoder are feedforward NNs with $32$ hidden
units and the LSTM comprises two layers with dimension $64$ and is followed by
a normalization layer. These numbers of neurons were chosen over larger ones
since we did not observe any significant increase in performance with larger
architectures. The input of each black-box module, $\bm{x}\in\mathbb{R}^{6}$,
gathers the solar irradiation on a horizontal surface and the time
information. NN-based models share a common learning rate of $5e^{-4}$ and
batch size of $4096$, and minimize the Mean Square Error (MSE) over a given
batch of data $B$:

$\displaystyle\mathcal{L}_{data}$ $\displaystyle=\frac{1}{|B|}\sum_{s\in
B}\left[\frac{1}{l_{s}}\sum_{k=0}^{l_{s}-1}\left[\frac{1}{m}\sum_{z=1}^{m}\xi^{z,s}_{k}\right]\right],$
(30) $\displaystyle\xi^{z,s}_{k}$
$\displaystyle=\left(\bm{T}^{z,s}_{k+1}-T^{z,s}(k+1)\right)^{2},$

where the subscript $s$ corresponds to which time series of length $l_{s}$ the
predictions and measured temperatures are taken from. To then validate the
performance of a model, we rely on the Mean Absolute Error (MAE), where

$\displaystyle\xi^{z,s}_{k}$
$\displaystyle=|\bm{T}^{z,s}_{k+1}-T^{z,s}(k+1)|,$

and the Mean Absolute Percentage Error (MAPE), with:

$\displaystyle\xi^{z,s}_{k}$
$\displaystyle=\frac{\bm{T}^{z,s}_{k+1}-T^{z,s}(k+1)}{T^{z,s}(k+1)}.$

The PCNNs are implemented in PyTorch [60] and were trained on NVIDIA P100
Graphics Processing Units (GPUs). The code and data will be made available on
https://gitlab.nccr-automation.ch/loris.dinatale/multi-zone-pcnns.

###  4.3 Benchmark models

To analyze the performance of the proposed PCNN architectures, we perform an
extensive ablation study and compare them to state-of-the-art gray- and black-
box methods. An overview of all the models used in this work, and whether they
are physically consistent, can be found in Table 1. Note that the linear and
LSTM models correspond to learning only the physics-inspired module of
S-PCNNs, respectively the black-box one. Furthermore, Res-cons is equivalent
to fitting both modules of S-PCNNs sequentially, showcasing the importance of
learning all the parameters of PCNNs simultaneously to attain state-of-the-art
accuracy.

|  | Physically
---|---|---
Category | Model | consistent
| Linear | ✓
Gray-box | Res | ✗
| Res-cons | ✓
| ARX | ✗
Black-box | ARX-KF | ✗
| LSTM | ✗
| PiNN | ✗
|  X-PCNN (Ours) | ✓
PCNNs |  M-PCNN (Ours) | ✓
|  S-PCNN (Ours) | ✓
Table 1: Summary of the models compared in this work, and whether they follow
the underlying physical laws.

####  4.3.1 Linear gray-box model

Firstly, it makes intuitive sense to investigate the accuracy of the physics-
inspired module of PCNNs on its own, leading to the following linear gray-box
model, hereafter referred to as the Linear model:

$\displaystyle\bm{T}_{k+1}$
$\displaystyle=\bm{T}_{k}+\bm{a}_{h}\odot\max{\\{\bm{u}_{k},\bm{0}\\}}$
$\displaystyle\quad+\bm{a}_{c}\odot\min{\\{\bm{u}_{k},\bm{0}\\}}$ (31)
$\displaystyle\quad-\bm{b}\odot(\bm{T}_{k}-T^{out}_{k})-\bm{\Delta
T}_{k}+\bm{e}\odot\bm{Q}^{win}_{k},$

where $\bm{Q}^{win}_{k}$ gathers the solar irradiation on the windows of each
zone in a vector, engineered from the measured irradiation on a horizontal
surface. Since there is no black-box module taking care of the impact of the
sun on building temperatures in this model, we indeed need to include it
manually. This can be done efficiently for UMAR but does not generalize to
arbitrary buildings, limiting the applications of such linear models, as
detailed in Appendix C. As for the other heat gains, $\bm{e}$ gathers the
trainable scaling parameters reflecting the impact of solar gains on each zone
temperature in a vector.

Since the classical least squares parameter identification gave rise to
physically inconsistent parameters, we chose to identify
$\bm{a}_{h},\bm{a}_{c},\bm{b},\bm{c},\bm{e}$ for each zone using Bayesian
Optimization (BO), as detailed in Appendix D. As for X-PCNNs, the heat
transfer coefficients between two adjacent thermal zones were then averaged
based on (10).

####  4.3.2 Residual models

A natural extension of the aforementioned linear model is to consider residual
models, where the idea is to fit the errors of the linear model predictions
with a black-box module to improve its performance. Assuming the linear model
in (31) to provide predictions $\bm{\hat{T}}_{k+1}$, a residual model fits a
function $f_{Res}:\mathbb{R}^{d^{\prime}+2m+1}\to\mathbb{R}^{m}$, typically
modeled with NNs, to the residual errors, i.e. it minimizes:

$\displaystyle\mathcal{L}_{Res}$ $\displaystyle=\frac{1}{|B|}\sum_{s\in
B}\left[\frac{1}{l_{s}}\sum_{k=0}^{l_{s}-1}\left[\frac{1}{m}\sum_{z=1}^{m}(\epsilon^{z,s}_{k})^{2}\right]\right],$
(32) $\displaystyle\epsilon^{z,s}_{k}$
$\displaystyle=f_{Res}^{z}(\bm{T}_{k}^{s},\bm{x}_{k}^{s},\bm{u}_{k}^{s},T^{out,s}_{k})$
$\displaystyle\qquad-\left(\bm{T}^{z,s}(k+1)-\bm{\hat{T}}_{k+1}^{z,s}\right),$

and then predicts temperatures as follows:

$\displaystyle\bm{T}_{k+1}$
$\displaystyle=\bm{\hat{T}}_{k+1}+f_{Res}(\bm{T}_{k},\bm{x}_{k},\bm{u}_{k},T^{out}_{k}).$
(33)

This model is dubbed Res in the rest of this paper, and $f_{Res}$ has the same
encoder-LSTM-decoder structure as the proposed PCNNs for fair comparisons.

Remarkably, such residual models cannot be ensured to respect the underlying
physical laws in general since $f_{Res}$ is not independent of zone
temperatures, power inputs, and ambient temperatures. Since they are composed
of a physics-inspired base model and a black-box module running in parallel,
as the proposed PCNN architectures, we can indeed use similar arguments to
prove their physical consistency (see Remark 6). Consequently, we also
investigate the performance of a physically consistent residual model in this
work, dubbed Res-cons, where the black-box function learning the residuals
only depends on $\bm{x}$, as PCNNs. This model hence fits a function $f_{Res-
cons}:\mathbb{R}^{d^{\prime}}\to\mathbb{R}^{m}$ to the residuals, trained
similarly to its physically inconsistent counterpart, with the following
temperature predictions:

$\displaystyle\bm{T}_{k+1}$ $\displaystyle=\bm{\hat{T}}_{k+1}+f_{Res-
cons}(\bm{x}_{k}).$ (34)

Note that residual models first fit the base model to the data and then use
black-box methods to fit the residual errors while PCNNs learn both modules
together. This also implies that the physics-inspired module reflects the main
dynamics of residual models while it only ensures the physical consistency of
PCNN architectures, letting more expressive functions like NNs capture the
main system dynamics.

####  4.3.3 Autoregressive model with exogenous inputs

As a first black-box method, we analyze the performance of an ARX model, where
autoregressive lags of the states and inputs are used to predict the next
state:

$\displaystyle\bm{T}_{k+1}$
$\displaystyle=\bm{\alpha}_{0}\bm{T}_{k}+\bm{\alpha}_{1}\bm{T}_{k-1}+\ldots+\bm{\alpha}_{\delta}\bm{T}_{k-\delta}$
$\displaystyle\quad+\bm{\beta}_{0}\bm{\hat{x}}_{k}+\bm{\beta}_{1}\bm{\hat{x}}_{k-1}+\ldots+\bm{\beta}_{\delta}\bm{\hat{x}}_{k-\delta},$
(35) $\displaystyle\bm{\hat{x}}_{k}$
$\displaystyle=[\bm{u}_{k},T^{out}_{k},Q^{sun}_{k}]^{T},$

where $Q^{sun}_{k}\in\mathbb{R}$ is the solar irradiation measurement on a
horizontal surface, and the parameters
$\bm{\alpha}_{0},\ldots,\bm{\alpha}_{\delta}\in\mathbb{R}^{m\text{x}m}$,
$\bm{\beta}_{0},\ldots,\bm{\beta}_{\delta}\in\mathbb{R}^{m\text{x}(m+2)}$ are
identified through least square regression using the scikit-learn library
[61]. In this work, we set $\delta=11$, i.e. we use information from the last
$3\text{\,}\mathrm{h}$ to define the next temperatures, similarly to the warm
start period of the proposed PCNN architectures.

For comparison purposes, we also implemented an advanced ARX model relying on
the statsmodels package [62], which includes Kalman smoothing and filtering
operations out-of-the-box and sets
$\bm{\beta}_{1},\ldots,\bm{\beta}_{\delta}=0$ so that only current information
on external inputs is used. Since the identification procedure was harder in
that case, we identified each zone $z$ separately, with:

$\displaystyle\bm{\hat{x}}_{k}^{z}$
$\displaystyle=[\bm{u}_{k}^{z},T^{out}_{k},Q^{sun}_{k},\bm{T}^{y_{1}}_{k},\bm{T}^{y_{2}}_{k},\ldots,\bm{T}^{y_{|\mathcal{N}(z)|}}_{k}]^{T},$

where $y_{1},\ldots,y_{|\mathcal{N}(z)|}\in\mathcal{N}(z)$ are the zones
adjacent to $z$, and the final model is dubbed ARX-KF. Note that ARX models
cannot be enforced to be physically consistent in general.

######  Remark 9 (Kalman filtering and smoothing).

Note that these operations could be included for any other model as well,
potentially impacting their performance. In this work, we focus on methods
working on unfiltered data, typically involving NNs, but we also provide this
ARX-KF as an example of what can be achieved with existing toolboxes on a
laptop compared to NN-based methods requiring access to GPUs for training.

####  4.3.4 Neural network models

As another natural ablation of PCNNs, we also investigate the quality of the
black-box module alone. Instead of treating the power inputs and temperatures
in a separate module, all the inputs are fed in the black-box function
$f_{LSTM}:\mathbb{R}^{d^{\prime}+2m+1}\to\mathbb{R}^{m}$, leading to the LSTM
model:

$\displaystyle\bm{T}_{k+1}$
$\displaystyle=\bm{T}_{k}+f_{LSTM}(\bm{T}_{k},\bm{u}_{k},\bm{x}_{k},T^{out}_{k}).$
(36)

As expected, such classical NN-based methods are naturally physically
inconsistent and might fail to capture the underlying physical laws even if
they fit the data well (Section 1).

Finally, we also compare PCNNs to a standard physics-informed NN, hereafter
the PiNN model, again relying on the same architecture as the black-box
modules of PCNNs and the LSTM model. However, as is classically done, its loss
function is modified to steer the learning toward physically meaningful
solutions, with:

$\displaystyle\mathcal{L}_{PiNN}$
$\displaystyle=\mathcal{L}_{data}+\lambda\mathcal{L}_{phys},$ (37)

where $\lambda$ is a tuning hyperparameter. Since the purpose of this
additional loss term $\mathcal{L}_{phys}$ is to capture physical
inconsistencies and penalize them, we naturally design it to bias the model
towards solutions satisfying the desired properties (1) and (2). Consequently,
we penalize negative gradients of the final predicted temperatures, i.e. at
time $l_{s}$, with respect to control inputs and ambient temperatures observed
along the horizon:

$\displaystyle\mathcal{L}_{phys}$ $\displaystyle=\frac{1}{|B|}\sum_{s\in
B}\left[\frac{1}{l_{s}}\sum_{k=0}^{l_{s}-1}\left[\frac{1}{m}\sum_{z=1}^{m}g^{z,s}_{k}\right]\right],$
(38) $\displaystyle g^{z,s}_{k}$
$\displaystyle=\sum_{y=1}^{m}\left[r\left(-\frac{\partial\bm{T}^{z,s}_{l_{s}}}{\partial\bm{u}^{y,s}_{k}}\right)\right]+r\left(-\frac{\partial\bm{T}^{z,s}_{l_{s}}}{\partial
T^{out,s}_{k}}\right),$ (39)

where $r(x)=\max\\{x,0\\}$, also known as the Rectified Linear Unit (ReLU)
function. Since we are interested in physically consistent models in this
work, we empirically fixed $\lambda=100$, which ensures the loss term is
dominated by the physical inconsistencies, thereby steering the PiNN towards
interesting solutions.

######  Remark 10 (Computational complexity).

To ensure a model is following the underlying physical laws at all times, one
should check the gradients throughout the prediction horizon, and not only for
the last predictions, as proposed in (38). However, since each gradient
computation requires one forward and one backward pass of the data, the
computational complexity grows linearly with the number of predictions to
analyze. Consequently, we only compute the gradients of the last predictions
with respect to all the control inputs and ambient temperatures observed along
the horizon to steer PiNNs, alleviating the associated computational burden.

##  5 Results

This section provides a comprehensive analysis of the models presented in
Table 1. In Section 5.1, we start by discussing the performance of the
different models on the validation data. In Section 5.2, we then provide a
visual and qualitative discussion of the NN-based models222The LSTM, PiNN, and
PCNN architectures. Despite also being composed of an NN, residual models are
not considered as NN-based models in this work since their main dynamics are
still captured by the underlying linear model, and not the NN., before
numerically investigating their physical consistency more in-depth in Section
5.3. Finally, Section 5.4 concludes with a brief overview of the computational
complexity associated with all the models. Altogether, this will allow us to
understand the trade-offs between the physical consistency, accuracy, and
computational complexity of the various data-driven building modeling methods
examined in this work.

Note that the NN-based models were run with several random seeds, and, unless
stated otherwise, the results discussed throughout this Section were obtained
using the best-performing seed in each case. Remarkably, however, this does
not impact our conclusions significantly since all the architectures proved to
be robust to the choice of random seed, with standard deviations in the range
of $0.01$ – $0.04$ and $0$ – $0.2\%$ for the MAE and MAPE, respectively, as
detailed in Appendix E.

###  5.1 Performance analysis

The best performance of all the analyzed methods in terms of MAE and MAPE on
more than $750$ three-day-long time series from the validation set is reported
in Table 2. As can be observed, all the proposed PCNN architectures attain
state-of-the-art accuracy, both in terms of MAE and MAPE. They are followed by
physically inconsistent black-box methods, especially the ones relying on very
expressive NNs. As expected, the least expressive class of methods, gray-box
models, performs the worst. Combining these results with the physical
consistency of each method, we can conclude that the proposed PCNN
architectures take the best out of both gray- and black-box methods, attaining
state-of-the-art performance while respecting the underlying physical laws by
construction without trade-off, making them ideal thermal building models.

|  | Phys. |  |
---|---|---|---|---
Category | Model | cons. | MAE | MAPE
| Linear | ✓ | $1.79$ | $7.5\%$
Gray-box | Res | ✗ | $1.79$ | $7.7\%$
| Res-cons | ✓ | $1.50$ | $6.4\%$
| ARX | ✗ | $1.68$ | $7.1\%$
Black-box | ARX-KF | ✗ | $1.35$ | $5.6\%$
| LSTM | ✗ | $1.27$ | $5.5\%$
| PiNN | ✗ | $1.37$ | $5.8\%$
| X-PCNN (Ours) | ✓ | $\bm{1.17}$ | $\bm{4.9\%}$
PCNNs |  M-PCNN (Ours) | ✓ | $1.25$ | $5.3\%$
|  S-PCNN (Ours) | ✓ | $1.22$ | $5.1\%$
Table 2: Physical consistency, MAE, and MAPE of the methods investigated in
this work.

Remarkably, enforcing the physical consistency of LSTMs, as in PCNNs, seems to
improve their accuracy in this case study despite the introduced constraints.
This confirms the ongoing trend in ML research to include prior knowledge in
NN architectures. Even if it might intuitively seem that introducing
structural constraints should hinder the expressiveness of LSTMs, these
results suggest that it can on the contrary be helpful. Moreover, one can draw
similar conclusions with the two residual models investigated in this work,
with Res-cons clearly outperforming its physically inconsistent counterpart
despite both models relying on the same linear basis. While these results are
not reported here, ensuring the black-box modules to be independent of
$D^{z}$, as mentioned in Remark 4, also increased the performance of X-PCNNs.
Altogether, these results point towards performance benefits of grounding NN
architectures in the underlying physics, ensuring that they learn meaningful
solutions.

As shown in Table 2, the residual models (Res and Res-cons), which are
conceptually close to PCNNs333Especially Res-cons, where the only difference
in architecture comes from the solar irradiation processing., are unable to
attain similar performance to PCNNs. This hints towards the benefits of
learning all the parameters together in an end-to-end fashion instead of first
identifying the linear part and then fitting the residual errors.

PCNNs are on average $30$ – $35\%$ and $17$ – $22\%$ more accurate than the
other physically consistent methods, namely the Linear and Res-cons model,
respectively. To visualize the error propagation of these methods over three
days, their MAE is plotted in Figure 7. This shows that the proposed PCNNs not
only perform better on average, but along the entire prediction horizon,
except during the first few hours, where the linear and residual model attain
similar performance. The main reason behind this behavior is the warm start of
PCNNs, which often gives erroneous first predictions, but they quickly make up
for it and show much stronger performance in the long run. At the end of the
horizon, PCNNs indeed show an error $34$ – $41\%$ and $10$ – $18\%$ lower than
Linear and Res-cons, respectively, with the best performance again achieved by
the X-PCNN.

The results of our case study suggest that the black-box modules of PCNNs are
able to process the raw solar irradiation data and infer its impact on the
zone temperatures. Indeed, they outperform gray-box models, which have access
to engineered solar gains (Appendix C). Interestingly, since the impact of the
sun is implicitly computed by their black-box modules, PCNNs could easily be
applied to any building, even when shading comes into play, making the
engineered preprocessing of solar data required for gray-box models much more
complex. Furthermore, while the only nonlinear gains considered in this paper
come from solar irradiation, the flexibility of the black-box module would
also allow it to learn other gains, such as the ones stemming from occupants.

Figure 7: MAE of all the physically consistent methods averaged over the three
zones and the time series of three days in the validation data set.

###  5.2 The necessity of physical consistency

Now that Section 5.1 established that PCNNs attain state-of-the-art
performance in terms of accuracy, even outperforming pure black-box methods,
let us visualize the behavior of one S-PCNN, one PiNN, and one LSTM for a
given random seed in Figure 8. To that end, the thermal power is turned off in
Zone $1$ and $2$ and we examine the impact of heating (red), cooling (blue),
or providing no power input (black) in Zone $3$. Note that the heating pattern
corresponds to the true power inputs measured in Zone $3$ in March $2021$,
which we mirror to create the cooling pattern.

Figure 8: Visualization of heat propagation for the proposed S-PCNN
architecture on the left compared to a PiNN in the middle and an LSTM on the
right. The bottom plots show the heating (red) and cooling (blue) patterns
applied to Zone $3$ while the power is turned off in Zone $1$ and $2$,
compared to the situation when no power is applied (black). The other plots
depict the corresponding temperature predictions of each model in each of the
three zones.

As one can immediately realize, following the laws of thermodynamics, heating
or cooling Zone $3$ increases, respectively decreases its temperature in the
S-PCNN model. This effect is then propagated to the adjacent Zone $2$, and
later to Zone $1$, impacting their temperatures even though they are neither
heated nor cooled, illustrating the effect of Corollary 1 ensuring physically
consistent predictions. Note that while only the temperature predictions of
one S-PCNN are pictured in Figure 8, similar effects were observed for X-PCNNs
and M-PCNNs. This is expected since all of them share the same properties,
i.e. the same physics-inspired module, to ensure they follow the criteria
(1)-(3).

On the other hand, all power inputs lead to almost indistinguishable
temperature predictions for the PiNN and LSTM. Despite achieving a very good
fit of the data (Section 5.1), these models are hence obviously flawed and can
be misleading in practical applications. We can sometimes even observe lower
temperature predictions when heating is turned on than when the zones are
cooled, a clear sign of physical inconsistency, e.g. in Zone $2$ at the end of
March $11$ and $12$. This clearly exemplifies the issue of shortcut learning
in the case of thermal modeling, where NNs manage to fit the data well without
understanding the underlying physics.

Furthermore, this illustrates how PiNNs only steer the learning towards
interesting solutions without providing any guarantees concerning the actual
behavior of the model. In fact, trained PiNNs always gave very similar
predictions to LSTMs in our experiments, as can also be seen in Figure 8,
hinting that modifying the loss function $\mathcal{L}_{PiNN}$ of the model did
not have much impact on the final solution found despite the large $\lambda$
used. While tuning this hyperparameter might lead to better results, it is a
notoriously cumbersome task and would still never guarantee the physical
consistency of the final model [63]. Thus, it was not considered in this work.

Very importantly, these results point out a somewhat counter-intuitive and
often overlooked characteristic of NNs: contrary to physics-based models, a
good fit to the data does not necessarily imply that the quality of the model
is good. In our case, the PiNNs and LSTMs were indeed able to fit the data
well without considering the impact of heating and cooling, i.e. solely
mapping external conditions to building temperatures. One hence has to be
careful when NNs are used to model physical systems and make sure the trained
models do not simply find shortcuts to fit the data well without respecting
the underlying physical laws. This calls for physically grounded
architectures, such as PCNNs, for applications where the physical consistency
of the model is critical.

We suspect that LSTMs and PiNNs were able to fit the data very well without
considering the impact of heating and cooling because of the specific data
used in this case study. Firstly, windows cover the entire East facade of
UMAR, rendering the building especially sensitive to solar gains and external
weather conditions. Secondly, while different controllers have been applied
during the collection period of the data set, all of them had the same
objective of maintaining the building temperature in a comfortable range and
hence reacted similarly to external conditions. Coupling these facts, it seems
indeed plausible to accurately predict building temperatures solely based on
external conditions and without considering heating and cooling inputs. This
might explain why the very expressive LSTMs and PiNNs found non-physical
shortcuts modeling the evolution of inside temperatures well.

Interestingly, Res did capture a much more significant impact of heating and
cooling, but remains completely oblivious to the underlying physics, with
cooling often resulting in higher temperatures than heating. This illustrates
the need to also consider physical consistency when designing residual models,
such as in the proposed Res-cons architecture. In general, all these results
hence suggest that physical consistency should always be considered when
dealing with NNs for physical systems.

Note that the identification of the ARX-KF model assigned a negative scaling
parameter for the power input to Zone $2$, for example, meaning that heating
this zone will lead to lower temperatures. We could observe similar issues
with the parameters of the classical ARX model, indeed confirming that ARX
models might not be physically consistent in practice, as claimed in Section
4.3.

###  5.3 Numerical analysis of physical consistency

This section investigates the gradients of the predictions of NN-based models
numerically, to strengthen the theoretical and visual claims in Table 2 and
Figure 8. Since gradients can be retrieved automatically through the
torch.autograd module [60], it allows us to numerically assess if the models
respect criteria (1) and (2), a necessary condition to ensure physical
consistency. Following Remark 10, we investigate the gradients of the
temperature predictions at the end of the three-day-long horizon with respect
to the power inputs and ambient temperatures observed at each time step. Note
that this corresponds to the gradients used to steer the learning of PiNNs in
(39), except for the X-PCNN, for which fewer gradients can be computed, as
detailed in Appendix F.

Figure 9: Distribution of the gradients of the temperatures at the end of the
prediction horizon with respect to power inputs and external temperatures
observed along the horizon for the NN-based models.

Overall, this gives us access to more than two million gradient values for
each model, except the X-PCNN, with slightly over one million values. The
resulting density histograms are shown in Figure 9, where one can directly
observe negative gradients only for the two black-box models not grounded in
the underlying physics. In fact, penalizing negative gradients in
$\mathcal{L}_{PiNN}$ decreased the magnitude of the PiNN gradients, steering
them to zero, but did not change the proportion of negative ones. In other
words, it did not improve the physical consistency of PiNNs since they still
violate conditions (1) and (2) as often as classical LSTMs. On the other hand,
the proposed PCNN architectures keep all the gradients that require positivity
in $\mathbb{R}_{+}$, as desired, providing a numerical argument supporting
their physical consistency.

###  5.4 Computational complexity

As final comparison metric between the models, Figure 10 presents the time
required by each model per training iteration. Importantly, these numbers are
subject to implementation considerations and hence have to be taken with a
grain of salt since we did not optimize the models. Nonetheless, all of them
used the same backbone architecture, which allows relative comparisons, e.g.
between the three proposed PCNNs, between the two residuals models, or between
LSTMs and PiNNs. Note that the linear and ARX models are not considered here
since their “training” procedure is very different: it does not require access
to a GPU and does not rely on gradient descent

Firstly, as expected, PiNNs take more time to run than classical LSTMs since
each batch has to be forwarded and backwarded through the networks twice, once
to compute the predictions used in $\mathcal{L}_{data}$ and another time to
calculate the gradients in $\mathcal{L}_{phys}$. Secondly, residual models
need access to the predictions of the underlying linear model at each step to
compute the residual errors before fitting them, which also entails a clear
computational overhead compared to LSTMs. Finally, the PCNN architectures all
require to compute both the black-box module output $D_{k}$ and the physics-
inspired module predictions $E_{k}$ at each step $k$ along the horizon, which
also entails additional overhead on top of classical black-box models.
Interestingly, this is comparable to what happens in residual models,
explaining to some extent why the latter and S-PCNNs require similar amounts
of resources.

Compared to S-PCNNs, M-PCNNs and X-PCNNs are significantly more
computationally intensive. This intuitively follows from the shared black-box
module of S-PCNNs reducing the number of parameters to fit compared to
M-PCNNs. On the other hand, X-PCNNs require learning several models separately
instead of everything together, which duplicates the computational overhead
needed to create and move data to the GPU at each iteration and leads to an
increased computational burden compared to M-PCNNs. Stemming from these
remarks, we would expect these differences to grow if we were to apply PCNNs
to larger buildings with more thermal zones.

Figure 10: Training time per iteration of the methods relying on a GPU.

######  Remark 11 (Parallelizing X-PCNNs).

The training times reported here correspond to the total time required to
train each model for one iteration, i.e. the sum of training times of single-
zone PCNNs in the case of X-PCNNs, to represent the total amount of
computations needed. In practice, however, the different single-zone PCNNs can
easily be trained in parallel since they are independent, which can
significantly decrease the effective training time of X-PCNNs (dividing it
approximately by three in our setting with three zones). This would make them
the fastest multi-zone PCNNs to deploy but at the cost of additional
computational complexity.

##  6 Conclusion

This work presented extensions of single-zone PCNNs to the multi-zone setting,
thereby providing fully data-driven control-oriented multi-zone building
thermal models. The main idea of PCNNs is to let a physics-inspired and a
black-box module run in parallel, the former guaranteeing the compliance of
the output with the underlying physical laws – the laws of thermodynamics in
the case of building temperature modeling – and the latter capturing unknown
nonlinear dynamics, typically relying on NNs.

The proposed multi-zone PCNNs respect the underlying physics by design and at
all times despite requiring little engineering, contrary to classical
physically consistent methods. On the other hand, they outperformed state-of-
the-art black-box methods in terms of accuracy on a case study, hinting that
the constrained structure introduced to ensure they follow some ground rules
does not hinder their expressiveness. Our analyses showed little difference
between S-, M-, and X-PCNNs in general, with S-PCNNs entailing the least
computational complexity and X-PCNNs attaining the best accuracy. Remarkably,
all of them showed significantly better performance than classical physically
consistent data-driven methods, with accuracy improvements of $30$ – $35\%$
and $17$ – $22\%$ compared to a linear and a residual model, respectively.
While these results were obtained on a specific building, the performance gap
suggests that this trend would be observed for other applications. PCNNs
should thus remain the best modeling choice in general, even if classical
black-box methods might attain a better accuracy on different data sets.

Our investigations also illustrated a well-known pitfall of classical PiNNs
and LSTMs, which can find shortcuts to fit the data well without respecting
the underlying physical laws. This exemplifies the need to not solely consider
the fit to the data as a measure of the quality of NNs but also ensure that
their predictions make sense from a physical point of view. Our findings hence
support the current trend to incorporate inductive biases, i.e. prior
knowledge, in NNs to alleviate their infamous generalization issues, leading
to more principled architectures like the proposed PCNNs.

In light of these results, PCNNs pave the way for NN-based methods that can
simultaneously provide state-of-the-art performance and physical guarantees.
Furthermore, while only solar irradiation measurements and time information
were fed to the black-box module of PCNNs throughout this study, showcasing
the ability of the proposed approach to handle highly nonlinear effects, other
inputs could be integrated straightforwardly.

Thanks to their flexibility, we hence believe PCNNs to be an essential step
towards the safe deployment of NNs in real-world applications, closing the
sim2real gap of advanced control algorithms, and hope to spark an interest
both in the building modeling community and beyond. It would indeed be
interesting to investigate the capabilities of PCNNs to model different
buildings, incorporate additional nonlinearities in their black-box modules or
rules in the physics-inspired ones, and tackle other complex physical systems.

## Acknowledgements

This research was supported by the Swiss National Science Foundation under
NCCR Automation, grant agreement 51NF40_180545, and in part by the Swiss Data
Science Center, grant agreement C20-13.

## Declaration of competing interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * International Energy Agency (2020) [IEA] International Energy Agency (IEA), Tracking Buildings 2020, https://www.iea.org/reports/tracking-buildings-2020, Accessed: 2021.05.28, 2020.
  * Commission et al. [2022] E. Commission, D.-G. for Energy, L. Kranzl, M. Fallahnejad, R. Büchele, A. Müller, M. Hummel, T. Fleiter, T. Mandel, M. Bagheri, G. Deac, C. Bernath, J. Miosga, C. Kiefer, J. Fragoso, S. Braungardt, V. Bürger, D. Spasova, J. Viegand, R. Naeraa, S. Forthuber, Renewable space heating under the revised Renewable Energy Directive : ENER/C1/2018-494 : final report, Publications Office, 2022.
  * Westermann and Evins [2019] P. Westermann, R. Evins,  Surrogate modelling for sustainable building design–a review,  Energy and Buildings 198 (2019) 170–186.
  * Deb and Schlueter [2021] C. Deb, A. Schlueter,  Review of data-driven energy modelling techniques for building retrofit,  Renewable and Sustainable Energy Reviews 144 (2021) 110990.
  * Svetozarevic et al. [2022] B. Svetozarevic, C. Baumann, S. Muntwiler, L. Di Natale, M. N. Zeilinger, P. Heer,  Data-driven control of room temperature and bidirectional EV charging using deep reinforcement learning: Simulations and experiments,  Applied Energy 307 (2022) 118127.
  * Lei et al. [2022] Y. Lei, S. Zhan, E. Ono, Y. Peng, Z. Zhang, T. Hasama, A. Chong,  A practical deep reinforcement learning framework for multivariate occupant-centric control in buildings,  Applied Energy 324 (2022) 119742.
  * Höfer et al. [2021] S. Höfer, K. Bekris, A. Handa, J. C. Gamboa, M. Mozifian, F. Golemo, C. Atkeson, D. Fox, K. Goldberg, J. Leonard, et al.,  Sim2Real in robotics and automation: Applications and challenges,  IEEE transactions on automation science and engineering 18 (2021) 398–400.
  * Kadian et al. [2020] A. Kadian, J. Truong, A. Gokaslan, A. Clegg, E. Wijmans, S. Lee, M. Savva, S. Chernova, D. Batra,  Sim2real predictivity: Does evaluation in simulation predict real-world performance?,  IEEE Robotics and Automation Letters 5 (2020) 6670–6677.
  * Drgoňa et al. [2020] J. Drgoňa, J. Arroyo, I. C. Figueroa, D. Blum, K. Arendt, D. Kim, E. P. Ollé, J. Oravec, M. Wetter, D. L. Vrabie, et al.,  All you need to know about model predictive control for buildings,  Annual Reviews in Control 50 (2020) 190–232.
  * Di Natale et al. [2022] L. Di Natale, B. Svetozarevic, P. Heer, C. N. Jones,  Near-optimal Deep Reinforcement Learning Policies from Data for Zone Temperature Control,  arXiv preprint arXiv:2203.05434 (2022).
  * Yang et al. [2021] T. Yang, L. Zhao, W. Li, J. Wu, A. Y. Zomaya,  Towards healthy and cost-effective indoor environment management in smart homes: A deep reinforcement learning approach,  Applied Energy 300 (2021) 117335.
  * Dawood et al. [2022] S. M. Dawood, A. Hatami, R. Z. Homod,  Trade-off decisions in a novel deep reinforcement learning for energy savings in HVAC systems,  Journal of Building Performance Simulation 15 (2022) 809–831.
  * Crawley et al. [2001] D. B. Crawley, L. K. Lawrie, F. C. Winkelmann, W. F. Buhl, Y. J. Huang, C. O. Pedersen, R. K. Strand, R. J. Liesen, D. E. Fisher, M. J. Witte, et al.,  EnergyPlus: creating a new-generation building energy simulation program,  Energy and buildings 33 (2001) 319–331, https://doi.org/10.1016/S0378--7788(00)00114--6.
  * Wetter and Haugstetter [2006] M. Wetter, C. Haugstetter,  Modelica versus TRNSYS–A comparison between an equation-based and a procedural modeling language for building energy simulation,  Proceedings of SimBuild 2 (2006).
  * Mazzeo et al. [2020] D. Mazzeo, N. Matera, C. Cornaro, G. Oliveti, P. Romagnoni, L. De Santoli,  EnergyPlus, IDA ICE and TRNSYS predictive simulation accuracy for building thermal behaviour evaluation by using an experimental campaign in solar test boxes with and without a PCM module,  Energy and Buildings 212 (2020) 109812, https://doi.org/10.1016/j.enbuild.2020.109812.
  * Harb et al. [2016] H. Harb, N. Boyanov, L. Hernandez, R. Streblow, D. Müller,  Development and validation of grey-box models for forecasting the thermal response of occupied buildings,  Energy and Buildings 117 (2016) 199–207.
  * Zhang et al. [2019] Z. Zhang, A. Chong, Y. Pan, C. Zhang, K. P. Lam,  Whole building energy model for HVAC optimal control: A practical framework based on deep reinforcement learning,  Energy and Buildings 199 (2019) 472–490.
  * Ascione et al. [2017] F. Ascione, N. Bianco, C. De Stasio, G. M. Mauro, G. P. Vanoli,  Artificial neural networks to predict energy performance and retrofit scenarios for any member of a building category: A novel approach,  Energy 118 (2017) 999–1017.
  * Bourdeau et al. [2019] M. Bourdeau, X. qiang Zhai, E. Nefzaoui, X. Guo, P. Chatellier,  Modeling and forecasting building energy consumption: A review of data-driven techniques,  Sustainable Cities and Society 48 (2019) 101533, https://doi.org/10.1016/j.scs.2019.101533.
  * Foucquier et al. [2013] A. Foucquier, S. Robert, F. Suard, L. Stéphan, A. Jay,  State of the art in building modelling and energy performances prediction: A review,  Renewable and Sustainable Energy Reviews 23 (2013) 272–288.
  * Maasoumy et al. [2014a] M. Maasoumy, M. Razmara, M. Shahbakhti, A. S. Vincentelli,  Handling model uncertainty in model predictive control for energy efficient buildings,  Energy and Buildings 77 (2014a) 377–392.
  * Maasoumy et al. [2014b] M. Maasoumy, M. Razmara, M. Shahbakhti, A. S. Vincentelli,  Selecting building predictive control based on model uncertainty,  in: 2014 American Control Conference, IEEE, pp. 404–411.
  * Li et al. [2021] Y. Li, Z. O’Neill, L. Zhang, J. Chen, P. Im, J. DeGraw,  Grey-box modeling and application for building energy simulations-A critical review,  Renewable and Sustainable Energy Reviews 146 (2021) 111174.
  * Arroyo et al. [2020] J. Arroyo, F. Spiessens, L. Helsen,  Identification of multi-zone grey-box building models for use in model predictive control,  Journal of Building Performance Simulation 13 (2020) 472–486.
  * Lin et al. [2012] Y. Lin, T. Middelkoop, P. Barooah,  Issues in identification of control-oriented thermal models of zones in multi-zone buildings,  in: 2012 IEEE 51st IEEE conference on decision and control (CDC), IEEE, pp. 6932–6937.
  * Shamsi et al. [2021] M. H. Shamsi, U. Ali, E. Mangina, J. O’Donnell,  Feature assessment frameworks to evaluate reduced-order grey-box building energy models,  Applied Energy 298 (2021) 117174, https://doi.org/10.1016/j.apenergy.2021.117174.
  * Shamsi et al. [2019] M. H. Shamsi, U. Ali, J. O’Donnell,  A generalization approach for reduced order modelling of commercial buildings,  Journal of Building Performance Simulation 12 (2019) 729–744, https://doi.org/10.1080/19401493.2019.1641554.
  * Berthou et al. [2014] T. Berthou, P. Stabat, R. Salvazet, D. Marchio,  Development and validation of a gray box model to predict thermal behavior of occupied office buildings,  Energy and Buildings 74 (2014) 91–100, https://doi.org/10.1016/j.enbuild.2014.01.038.
  * Sheng and Jia [2020] F. Sheng, L. Jia,  Short-term load forecasting based on SARIMAX-LSTM,  in: 2020 5th International Conference on Power and Renewable Energy (ICPRE), IEEE, pp. 90–94.
  * Li and Wen [2014] X. Li, J. Wen,  Review of building energy modeling for control and operation,  Renewable and Sustainable Energy Reviews 37 (2014) 517–537.
  * Royer et al. [2016] S. Royer, S. Thil, T. Talbert,  Towards a generic procedure for modeling buildings and their thermal zones,  in: 2016 IEEE 16th International Conference on Environment and Electrical Engineering (EEEIC), IEEE, pp. 1–6, https://doi.org/10.1109/EEEIC.2016.7555567.
  * Tien et al. [2022] P. W. Tien, S. Wei, J. Darkwa, C. Wood, J. K. Calautit,  Machine Learning and Deep Learning Methods for Enhancing Building Energy Efficiency and Indoor Environmental Quality–A Review,  Energy and AI (2022) 100198\.
  * Deb et al. [2017] C. Deb, F. Zhang, J. Yang, S. E. Lee, K. W. Shah,  A review on time series forecasting techniques for building energy consumption,  Renewable and Sustainable Energy Reviews 74 (2017) 902–924.
  * Abbasimehr and Paki [2022] H. Abbasimehr, R. Paki,  Improving time series forecasting using LSTM and attention models,  Journal of Ambient Intelligence and Humanized Computing 13 (2022) 673–691.
  * Zou et al. [2020] Z. Zou, X. Yu, S. Ergan,  Towards optimal control of air handling units using deep reinforcement learning and recurrent neural network,  Building and Environment 168 (2020) 106535.
  * Delcroix et al. [2021] B. Delcroix, J. L. Ny, M. Bernier, M. Azam, B. Qu, J.-S. Venne,  Autoregressive neural networks with exogenous variables for indoor temperature prediction in buildings,  in: Building Simulation, volume 14, Springer, pp. 165–178.
  * Szegedy et al. [2013] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, R. Fergus,  Intriguing properties of neural networks,  arXiv preprint arXiv:1312.6199 (2013).
  * D’Amour et al. [2020] A. D’Amour, K. Heller, D. Moldovan, B. Adlam, B. Alipanahi, A. Beutel, C. Chen, J. Deaton, J. Eisenstein, M. D. Hoffman, et al.,  Underspecification presents challenges for credibility in modern machine learning,  arXiv preprint arXiv:2011.03395 (2020).
  * Di Natale et al. [2022] L. Di Natale, B. Svetozarevic, P. Heer, C. N. Jones,  Physically consistent neural networks for building thermal modeling: theory and analysis,  Applied Energy 325 (2022) 119806.
  * Geirhos et al. [2020] R. Geirhos, J.-H. Jacobsen, C. Michaelis, R. Zemel, W. Brendel, M. Bethge, F. A. Wichmann,  Shortcut learning in deep neural networks,  Nature Machine Intelligence 2 (2020) 665–673.
  * Beery et al. [2018] S. Beery, G. Van Horn, P. Perona,  Recognition in terra incognita,  in: Proceedings of the European conference on computer vision (ECCV), pp. 456–473.
  * Rosenfeld et al. [2018] A. Rosenfeld, R. Zemel, J. K. Tsotsos,  The elephant in the room,  arXiv preprint arXiv:1808.03305 (2018).
  * Heuer et al. [2016] H. Heuer, C. Monz, A. W. Smeulders,  Generating captions without looking beyond objects,  arXiv preprint arXiv:1610.03708 (2016).
  * Zech et al. [2018] J. R. Zech, M. A. Badgeley, M. Liu, A. B. Costa, J. J. Titano, E. K. Oermann,  Variable generalization performance of a deep learning model to detect pneumonia in chest radiographs: a cross-sectional study,  PLoS medicine 15 (2018) e1002683.
  * Karniadakis et al. [2021] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang,  Physics-informed machine learning,  Nature Reviews Physics 3 (2021) 422–440.
  * Daw et al. [2017] A. Daw, A. Karpatne, W. Watkins, J. Read, V. Kumar,  Physics-guided neural networks (pgnn): An application in lake temperature modeling,  arXiv preprint arXiv:1710.11431 (2017).
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational physics 378 (2019) 686–707.
  * Greydanus et al. [2019] S. Greydanus, M. Dzamba, J. Yosinski,  Hamiltonian neural networks,  Advances in neural information processing systems 32 (2019).
  * Lutter et al. [2019] M. Lutter, C. Ritter, J. Peters,  Deep lagrangian networks: Using physics as model prior for deep learning,  arXiv preprint arXiv:1907.04490 (2019).
  * Cranmer et al. [2020] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, S. Ho,  Lagrangian neural networks,  arXiv preprint arXiv:2003.04630 (2020).
  * Hendriks et al. [2020] J. Hendriks, C. Jidling, A. Wills, T. Schön,  Linearly constrained neural networks,  arXiv preprint arXiv:2002.01600 (2020).
  * LeCun et al. [2010] Y. LeCun, K. Kavukcuoglu, C. Farabet,  Convolutional networks and applications in vision,  in: Proceedings of 2010 IEEE international symposium on circuits and systems, IEEE, pp. 253–256.
  * Karpathy et al. [2015] A. Karpathy, J. Johnson, L. Fei-Fei,  Visualizing and understanding recurrent networks,  arXiv preprint arXiv:1506.02078 (2015).
  * Scarselli et al. [2008] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, G. Monfardini,  The graph neural network model,  IEEE transactions on neural networks 20 (2008) 61–80.
  * Gokhale et al. [2022] G. Gokhale, B. Claessens, C. Develder,  Physics informed neural networks for control oriented thermal modeling of buildings,  Applied Energy 314 (2022) 118852.
  * Nagarathinam et al. [2022] S. Nagarathinam, Y. S. Chati, M. P. Venkat, A. Vasan,  PACMAN: physics-aware control MANager for HVAC,  in: Proceedings of the 9th ACM International Conference on Systems for Energy-Efficient Buildings, Cities, and Transportation, pp. 11–20.
  * Drgoňa et al. [2021] J. Drgoňa, A. R. Tuor, V. Chandan, D. L. Vrabie,  Physics-constrained deep learning of multi-zone building thermal dynamics,  Energy and Buildings 243 (2021) 110992, https://doi.org/10.1016/j.enbuild.2021.110992.
  * Werbos [1990] P. J. Werbos,  Backpropagation through time: what it does and how to do it,  Proceedings of the IEEE 78 (1990) 1550–1560.
  * Empa [2021] Empa, NEST, https://www.empa.ch/web/nest/overview, 2021\. Accessed: 01.10.2022.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala,  PyTorch: An Imperative Style, High-Performance Deep Learning Library,  in: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, R. Garnett (Eds.), Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay,  Scikit-learn: Machine Learning in Python,  Journal of Machine Learning Research 12 (2011) 2825–2830.
  * Seabold and Perktold [2010] S. Seabold, J. Perktold,  statsmodels: Econometric and statistical modeling with python,  in: 9th Python in Science Conference.
  * Faroughi et al. [2022] S. A. Faroughi, N. Pawar, C. Fernandes, S. Das, N. K. Kalantari, S. K. Mahjour,  Physics-Guided, Physics-Informed, and Physics-Encoded Neural Networks in Scientific Computing,  arXiv preprint arXiv:2211.07377 (2022).
  * Nogueira [14 ] F. Nogueira, Bayesian Optimization: Open source constrained global optimization tool for Python, 2014–.

## Appendices

##  Appendix A Proofs of the main theoretical results

###  A.1 Proof of Proposition 1

The proof works by induction on $i$. Based on (20), we can immediately write,
$\forall z,y\in\mathcal{B}$:

$\displaystyle\frac{\partial\bm{T}_{k+j+1}^{z}}{\partial\bm{T}_{k+j}^{y}}$
$\displaystyle=\begin{cases}1-b^{z}-\sum_{y\in\mathcal{N}(z)}c^{zy},&\text{if
}y=z,\\\ c^{zy},&\text{if }y\in\mathcal{N}(z),\\\
0,&\text{otherwise,}\end{cases}$ (40)

where we used the definition of $\bm{\Delta T}$ in (12). By definition, if
(22) and (23) hold, we hence get positive derivatives if $y=z$ or
$y\in\mathcal{N}(z)$ and zeros for any other choice of $y$, satisfying (21)
and completing the base case of the induction.

Let us now assume that:

$\displaystyle\frac{\partial\bm{T}_{k+h}^{x}}{\partial\bm{T}_{k+j}^{y}}\geq
0,\ \forall y,x\in\mathcal{B},\ \forall j<h<i,$ (41)

with equality if and only if $y\notin\mathcal{N}^{(h-j)}(x)$, and show that
the proposition holds for time step $i$. Since we know the temperature in zone
$z$ at time $k+i$ is potentially impacted by the temperature in the entire
building at the previous step, we can decompose the partial derivative of
interest as follows:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j}^{y}}$
$\displaystyle=\sum_{x\in\mathcal{B}}\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+i-1}^{x}}\frac{\partial\bm{T}_{k+i-1}^{x}}{\partial\bm{T}_{k+j}^{y}},$
(42)

for all $y,z\in\mathcal{B}$. Since (20) is time-invariant, we know that:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+i-1}^{x}}$
$\displaystyle=\frac{\partial\bm{T}_{k+j+1}^{z}}{\partial\bm{T}_{k+j}^{x}}\geq
0,$

with equality if and only if $x\notin\mathcal{N}(z)$ by the base case of the
induction (40) if (22) and (23) hold. Similarly, by the induction hypothesis
(41), we know that:

$\displaystyle\frac{\partial\bm{T}_{k+i-1}^{x}}{\partial\bm{T}_{k+j}^{y}}\geq
0,\ \forall y,x\in\mathcal{B},$

with equality if and only if $y\notin\mathcal{N}^{(i-j-1)}(x)$. Putting the
last two equations together, we see that:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j}^{y}}$
$\displaystyle\geq 0,$

with equality only if each term of the sum in Equation (42) is zero. By the
previous arguments, this means $y\notin\mathcal{N}^{(i-j-1)}(x)$ or
$x\notin\mathcal{N}(z)$ for all zones $x$. This is equivalent to say that
there is no path from $y$ to $z$ in $(i-j)$ steps, i.e.
$y\notin\mathcal{N}^{(i-j)}(z)$, which concludes the inductive step.

###  A.2 Proof of Proposition 2

We start by noticing that $\forall y\in\mathcal{B}$, (20) implies:

$\displaystyle\frac{\partial\bm{T}_{k+j+1}^{y}}{\partial\bm{u}_{k+j}^{y}}$
$\displaystyle=\begin{cases}a^{y}_{h},&\text{if }\bm{u}_{k+j}^{y}>0,\\\
a^{y}_{c},&\text{if }\bm{u}_{k+j}^{y}<0,\\\ 0,&\text{otherwise},\end{cases}$
(43)
$\displaystyle\frac{\partial\bm{T}_{k+j+1}^{x}}{\partial\bm{u}_{k+j}^{y}}$
$\displaystyle=0\ \qquad\forall x\in\mathcal{B},x\neq y,$ (44)
$\displaystyle\frac{\partial\bm{T}_{k+j+1}^{y}}{\partial T^{out}_{k+j}}$
$\displaystyle=b^{y},$ (45)

Note that this proves that (24) and (25) for the case $i=j+1$ if
$a^{y}_{h},a^{y}_{c},b^{y}>0$, $\forall y\in\mathcal{B}$. When $i>j+1$,
Proposition 1 implies:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j+1}^{y}}\geq
0,\qquad\forall z,y\in\mathcal{B},\ \forall 0\leq j<i-1,$ (46)

with equality if and only if $y\notin\mathcal{N}^{(i-j-1)}(z)$ if the
conditions in (22) and (23) hold.

Relying on the fact that the temperatures at time $k+i$ are potentially
influenced by the temperatures in the whole building at time $k+j+1$, we have:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{u}_{k+j}^{y}}$
$\displaystyle=\sum_{x\in\mathcal{B}}\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j+1}^{x}}\frac{\partial\bm{T}_{k+j+1}^{x}}{\partial\bm{u}_{k+j}^{y}},$
(47)
$\displaystyle=\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j+1}^{y}}\frac{\partial\bm{T}_{k+j+1}^{y}}{\partial\bm{u}_{k+j}^{y}}\geq
0,$ (48)

where the second equality follows from (44) and the inequality holds as long
as (22) and (23) are respected and $a^{y}_{h},a^{y}_{c}>0$, $\forall
y\in\mathcal{B}$. Furthermore, by Proposition 1, equality is only reached if
$y\notin\mathcal{N}^{(i-j-1)}(z)$.

Similarly, we have:

$\displaystyle\frac{\partial\bm{T}_{k+i}^{z}}{\partial T^{out}_{k+j}}$
$\displaystyle=\sum_{y\in\mathcal{B}}\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j+1}^{y}}\frac{\partial\bm{T}_{k+j+1}^{y}}{\partial
T^{out}_{k+j}},$ (49)
$\displaystyle=\sum_{y\in\mathcal{B}}\frac{\partial\bm{T}_{k+i}^{z}}{\partial\bm{T}_{k+j+1}^{y}}b^{y}>0,$
(50)

where the strict inequality is respected as long as $b^{y}>0$, $\forall
y\in\mathcal{B}$. Indeed, since $z\in\mathcal{N}(z)$ by definition,
Proposition 1 then implies that at least one of the terms in the sum is
strictly positive, while the others are nonnegative.

##  Appendix B Details on the data processing

Once the data had been subsampled and processed as in [39, App. C] and
discarding the $23$% of incomplete measurements, i.e. where at least the
information from one sensor is missing, we were left with over
$80^{\prime}000$ data points. Since the proposed PCNN architectures are based
on NNs in our implementations (see Section 4.2), they are not able to handle
missing values, which prompted us to create a data set of time series without
missing values.

As we aimed to design models that are able to predict the temperature dynamics
over three day-long horizons, we truncated each sequence to a maximum of three
days, and separated the heating and cooling seasons. We allowed the time
series to overlap each hour, i.e. each four steps, to increase the data
efficiency of the approach. Finally, since we implemented a warm-start period
of $3\text{\,}\mathrm{h}$ for all the models, we also made sure the last
$3\text{\,}\mathrm{h}$ of data exist for each time series. To avoid very short
time series, we also ensured they always span at least
$12\text{\,}\mathrm{h}$. Altogether, this allowed us to create more than
$11^{\prime}000$ sequences of data without missing values, which were split in
a training and a validation set with proportions $80\%$-$20\%$, respectively
denoted $\mathcal{D}_{t}$ and $\mathcal{D}_{v}$, and where
$\mathcal{D}_{t}\cap\mathcal{D}_{v}=\emptyset$.

##  Appendix C Solar irradiation preprocessing

To compute the solar irradiation on the windows of a thermal zone $z$ from the
measured irradiation on a horizontal surface $Q^{sun}$, we rely on the
altitude and azimuth angles, respectively $\phi$ and $\theta$, of the sun. The
former captures the elevation of the sun above the horizon while the latter
represents its deviation from the north, in the clockwise direction.

Figure 11: Sketch of the azimuth angles used to compute the solar irradiation
on the windows of a building from the irradiation on a fixed vertical surface.

Firstly, using the altitude of the sun and basic trigonometry, one can easily
show that the measured irradiation on a horizontal surface corresponds to
$Q^{sun}=I\sin\phi$, where $I$ is the global solar irradiation. Similarly, we
know that the irradiation on a vertical surface following the sun, i.e.
tracking its azimuth angle to stay perpendicular to the incoming rays, can be
computed as $I^{vert}=I\cos\phi$. We can hence write the solar irradiation on
a vertical surface following the sun as follows:

$\displaystyle I^{vert}$ $\displaystyle=Q^{sun}\frac{\cos\phi}{\sin\phi}.$
(51)

Since building facades and windows have a fixed orientation in practice and do
not follow the sun azimuth, we again use basic trigonometry to compute the
irradiation on a north-south aligned surface facing east as
$I^{vert}\sin\theta$. Finally, if the facade is not exactly facing east, we
also need to account for its own ”azimuth” $\theta_{0}$, i.e. how much it is
rotated clockwise starting from an east-facing position (Figure 11), which
leads to:

$\displaystyle Q^{win}$ $\displaystyle=I^{vert}\sin(\theta-\theta_{0})$ (52)
$\displaystyle=Q^{sun}\frac{\cos\phi}{\sin\phi}\sin(\theta-\theta_{0}).$ (53)

Once this has be done for each zone $z$, we can populate the required vector
$\bm{Q}^{win}$ used by gray-box architectures in this work.

As one can readily observed, this processing only requires access to the
elevation and azimuth angles of the sun, and to the orientation of the facade
of interest. Furthermore, both solar angles solely depend on the geographical
position of the building, i.e. its latitude and longitude, and the time at
which the measurement was taken. The position and orientation of a building
can easily be found on plans or Google Maps, and we used the Astral Python
library (https://astral.readthedocs.io/en/latest/) to compute the solar angles
corresponding to each time step in our data.

Note that, while this processing works very well for unobstructed facades when
its orientation is known, it cannot be used when for example other buildings
or trees exist in front of the windows and create shading patterns. In that
case, one has to rely on architectures which are able to automatically process
horizontal solar irradiation measurements depending on time information, such
as the LSTMs used in the black-box module of PCNNs. Nonetheless, we can use it
in this paper since UMAR is not obstructed, leading to a very efficient
computation of the true solar irradiation patterns on the windows of each
zone.

##  Appendix D Linear model identification

As is classically done in linear system identification, we firstly used the
least squares method to find the parameters $a^{z}_{h}$, $a^{z}_{c}$, $b^{z}$,
$c^{zy}$, $e^{z}$ best fitting the training data for each thermal zone $z$ and
neighboring zone $y\in\mathcal{N}(z)$, such as in [39, App. A.2]. However,
ensuring none of these parameters is negative, which is necessary to respect
the underlying physics, produced $c^{23}=0$. This is clearly not physically
meaningful, as it would mean there is no heat transfer from Zone $3$ to Zone
$2$. Consequently, we also implemented a BO framework, relying on the
bayes_opt Python library [64]. This allowed us to extensively search for the
best physically consistent parameters for each zone over a five-step
prediction horizon, constraining all the parameters to be positive, for $2300$
iterations starting with $200$ random initial points.

##  Appendix E Impact of the random seed

The mean performance of the five NN-based model architectures, as well as the
corresponding standard deviation, is presented in Table 3. While the LSTM and
S-PCNN models were run on five seeds due to their slightly higher sensitivity,
the other results were obtained over three seeds leading to very consistent
performance. As in the the original PCNN paper [39], this hints at the
robustness of the proposed approach, which does not seem significantly
impacted by the random seed, or at least similarly to classical NN models. On
this case study, the proposed X-PCNN and M-PCNN seem more robust to the choice
of random seed than the S-PCNN. However, as shown in Table 2, the latter
sometimes outperforms M-PCNNs. Nonetheless, overall, X-PCNNs seem to have the
upper hand, always attaining state-of-the-art performance even under different
random seeds.

Model | MAE | MAPE
---|---|---
LSTM | $1.33\pm 0.04$ | $5.7\%\pm 0.2\%$
PiNN | $1.38\pm 0.01$ | $5.9\%\pm 0.1\%$
X-PCNN (Ours) | $\bm{1.18\pm 0.01}$ | $\bm{5.0\%\pm 0.0\%}$
M-PCNN (Ours) | $1.26\pm 0.01$ | $5.4\%\pm 0.0\%$
S-PCNN (Ours) | $1.27\pm 0.04$ | $5.4\%\pm 0.2\%$
Table 3: MAE and MAPE of the methods investigated in this work, average over
the three thermal zone, the three-day long horizon, and more than $750$ time
series.

##  Appendix F X-PCNN gradients

In the case of X-PCNNs, at inference time, we use each single-zone PCNN to
predict the next temperature in the corresponding zone. The new temperatures
in the building are then updated in the data of all the single-zone PCNNs so
they can predict the next step. This is required because the single-zone PCNNs
cannot evolve independently over the prediction horizon since they depend on
temperatures in neighboring zones at teach step. However, overwriting the data
at each step breaks the automatic backpropagation of Python, and we cannot
automatically compute the gradient of the temperature in zone $z$ with respect
to power inputs or temperatures in another zone $z^{\prime}$ without
implementation overhead. We can only retrieve gradients with respect to each
single-zone PCNN’s inputs, i.e. the power $u^{z}$, and the ambient
temperature. Even if we can only compute parts of the gradients automatically,
we still show them in Figure 9 for reference. Note that as we already know
X-PCNNs are physically consistent since they satisfies the criteria of
Corollary 1, these implementation considerations do not put the architecture
in jeopardy.