arXiv:1608.03090v1 [] 10 Aug 2016

Regression Models for Output Prediction of
Thermal Dynamics in Buildings∗
Georgios C. Chasparis†
Software Competence Center Hagenberg GmbH
Department of Data Analysis Systems
Softwarepark 21
4232 Hagenberg, Austria
Email: georgios.chasparis@scch.at
Thomas Natschlaeger
Software Competence Center Hagenberg GmbH
Department of Data Analysis Systems
Softwarepark 21
4232 Hagenberg, Austria
Email: thomas.natschlaeger@scch.at

Standard (black-box) regression models may not necessarily
suffice for accurate identification and prediction of thermal
dynamics in buildings. This is particularly apparent when
either the flow rate or the inlet temperature of the thermal
medium varies significantly with time. To this end, this paper analytically derives, using physical insight, and investigates linear regression models with nonlinear regressors for
system identification and prediction of thermal dynamics in
buildings. Comparison is performed with standard linear regression models with respect to both a) identification error,
and b) prediction performance within a model-predictivecontrol implementation for climate control in a residential
building. The implementation is performed through the EnergyPlus building simulator and demonstrates that a careful
consideration of the nonlinear effects may provide significant
benefits with respect to the power consumption.
1

Introduction
The increased demand for electricity power and/or fuel
consumption in residential buildings requires an efficient
control design for all heating/cooling equipment. To this end,
recently, there have been several efforts towards a more efficient climate control in residential buildings [2, 3, 4, 5]. Ef-

∗ This work is part of the mpcEnergy project which is supported within
the program Regionale Wettbewerbsfähigkeit OÖ 2007-2013 by the European Fund for Regional Development as well as the State of Upper Austria.
The research reported in this article has been (partly) supported by the Austrian Ministry for Transport, Innovation and Technology, the Federal Ministry of Science, Research and Economy, and the Province of Upper Austria
in the frame of the COMET center SCCH. An earlier version of part of this
paper appeared in [1].
† Address all correspondence related to this paper to this author.

ficiency of power/fuel consumption is closely related to the
ability of the heating/cooling strategy to incorporate predictions of the heat-mass transfer dynamics and exogenous thermal inputs (e.g., outdoor temperature, solar radiation, user
activities, etc.). Naturally this observation leads to modelpredictive control (MPC) implementations (see, e.g., [2, 3]).
The performance though of any such implementation will be
closely related to the performance of the prediction models
for the evolution of both the indoor temperature and the exogenous thermal inputs.
A common approach in the derivation of prediction
models for the indoor temperature evolution relies on the
assumption that the heat-mass transfer dynamics can be
approximated well by linear models. For example, in
the optimal supervisory control formulation introduced in
[2] for controlling residential buildings with an heatingventilation-air-conditioning (HVAC) system, a linear (statespace) model is considered. As pointed out in [2], including
the actual dynamics might be complicated for the derivation
of optimal control strategies. Similar is the assumption in the
dynamic game formulation of HVAC thermally controlled
buildings in [6], where again a linear model is assumed for
the overall system. Also, in references [7, 8], detailed nonlinear representations of the heat-mass transfer dynamics (of
a vapor-compression system in [7] and for a building in [8])
are linearized about an operating point to allow for a more
straightforward controller design.
Given this difficulty in incorporating nonlinear representations in the controller design, several efforts for identification of heat-mass transfer dynamics have adopted linear
transfer functions, such as the standard ARX or ARMAX
black-box model structures (cf., [9]). Examples of such im-

plementations include the MIMO ARMAX model in [10],
the MISO ARMAX model for HVAC systems in [11] and
the ARX model structure considered in [12].
On the other hand, the nonlinear nature of the heatmass transfer dynamics in buildings has been pointed out by
several papers, including the switched linear dynamics for
modeling the intermittent operation of RH systems in [4],
the multiple-time scale analysis considered in [12, 5] and
the physics-based representation of the air-conditioning systems in [7]. In reference [12], a comparison is performed
between standard linear ARX models with two time-scale
transfer models, which according to the authors better represent thermal models. Indeed, the operation of either an
radiant-heating (RH) system or an HVAC system usually
takes place in a faster time-scale than the room dynamics,
an observation that was also utilized for the derivation of a
(nonlinear) model predictive controller in [5] using singularperturbation techniques.
The following questions naturally emerge: What is the
effect of the nonlinearities of an RH or an HVAC system in
the identification error of the indoor temperature? Can more
detailed model representations be utilized to reduce identification error and improve energy efficiency? This paper begins with the observation that thermal dynamics in buildings
are nonlinear in nature due to the operation of the RH and
HVAC systems. We avoid introducing any multiple timescale approximation as in [5] and we adopt a detailed representation of the heat-mass transfer dynamics (using standard
Bond-graph models). In principle, linear transfer models
incorporating linear regressors, such as Output-Error (OE),
ARX or ARMAX structures, may be sufficient for system
identification in buildings. However, we wish to investigate
whether this is a reasonable assumption and whether more
accurate representations of the dynamics are necessary.
To this end, linear regression models with nonlinear regressors are derived analytically using physical insight for
prediction of the indoor temperature in residential buildings
with RH and/or HVAC systems. The derivation accounts for
alternative information structures depending on the sensors
available, which is the main contribution of this paper. The
derived identification models are compared against standard
(black-box) linear regression models trained with simulated
data generated in EnergyPlus (V7-2-0) building simulator
tool developed by the U.S. Department of Energy [13]. Comparison is also performed with respect to the prediction performance of the derived models when used within a standard
MPC for climate control in a residential building. The results indicate that a more detailed design of the thermal-mass
transfer prediction models can be beneficial with respect to
the energy consumption. This paper extends prior work of
the same authors [1], since it extends the derivation of the regression models into a larger family of information structures
(sensors available), while it provides a detailed comparison
with respect to the energy consumption.
In the remainder of the paper, Section 2 provides a description of the overall framework and the objective of this
paper. Section 3 provides a short introduction to a class of
linear regression models (briefly LRM) for system identifica-

tion of heat-mass transfer dynamics. Section 4 analytically
derives using physical insight linear regression models with
nonlinear regressors (briefly NRM). In Section 5, we compare the performance of the derived NRM with LRM with
respect to both a) identification error, and b) prediction performance within a standard MPC for climate control. Finally,
Section 6 presents concluding remarks.
Notation:
• col{x1 , ..., xn }, for some real numbers x1 , ..., xn ∈ R,
denotes the column vector (x1 , ..., xn ) in Rn . Also,
row{x1 , ..., xn } denotes the corresponding row vector.
• diag{x1 , ..., xn } denotes the diagonal matrix in Rn×n
with diagonal entries x1 , ..., xn .
• for any finite set A, |A| denotes the cardinality of A.
.
• = denotes definition.
2

Framework & Objective
When addressing system identification and control of
heat-mass transfer dynamics in buildings, the assumed model
of indoor temperature prediction should be a) accurate
enough to capture the phenomena involved, and b) simple
enough to accommodate a better control design. Given that
prior literature has primarily focused on linear state-space or
transfer models, we wish to address the following:
1. investigate the performance of standard (black-box)
linear (regression) models (LRM) in the identification/prediction of heat-mass transfer dynamics in residential buildings;
2. derive analytically, using physical insight and under alternative information structures, an accurate representation of the dynamics using linear regression models with
nonlinear regressors (NRM);
3. provide a detailed comparison between the derived
NRM with standard LRM and assess the potential energy saving.
To this end, we consider residential buildings divided
into a set of zones I . Each zone is controlled independently
in terms of heating using either a RH and/or an HVAC system. We will use the notation i, j to indicate representative
elements of the set of zones I . The temperature of a zone i is
affected by the temperature of the neighboring zones and/or
the outdoor environment, denoted Ni , i.e., those zones which
are adjacent to i, separated by some form of separator.
In the remainder of this section, we provide some standard background on heat-mass transfer dynamics (based on
which the derivation of the prediction models will be presented). We further discuss the considered assumptions with
respect to the available sensors, and finally some generic
properties of the dynamics.
2.1

Background: Heat-mass transfer dynamics
We will use standard modeling techniques of heat-mass
transfer dynamics for each thermal zone i ∈ I , following
a standard Bond-graph approach [14] and the simplified
Fourier’s law for heat conduction and Newton’s law for heat

Fig. 1: Bond-graph approximation of heat-mass transfer dynamics of a thermal zone i ∈ I .
convection (cf., [15]). Such bond-graph approximations of
thermal dynamics are equivalent to a series of electrical
resistance-capacitor (RC) connections. Figure 1 provides a
detailed bond-graph approximation for a thermal zone i ∈ I
with both RH and HVAC heating equipment.

heat transfer dynamics for a separator i j can be written as

Q̇−
s,i j
We use C to denote the thermal capacitance of a separator (e.g., wall) or a thermal zone (e.g., room), and the symbol
R to denote the heat resistance, according to an electricity
equivalent terminology. We also use the symbol Q̇ to denote
heat flow and the symbol T to denote temperature.
The subscript s denotes separator-related parameters;
the subscript r denotes zone-related parameters, and usually
corresponds to a room; the subscript w denotes parameters
related to the RH system; and the subscript a denotes HVACrelated parameters. Finally, a pair i j denotes the interconnection between two neighboring zones i, j ∈ I .
Separators have been modeled as two thermal resistances separated by a capacitance, with the superscript “+”
denoting the “inside” part of the separator, and the superscript “−” denoting the “outside” part of the separator. For
example, R+
s,i j denotes the thermal resistance that it is adjacent to i. Furthermore, Tr,i , i ∈ I , denotes the temperature
of zone i and Ts,i j , j ∈ Ni , denotes the (internal) temperature of the separator i j. The separator exchanges heat with
the neighboring zones i and j due to heat convection, whose
−
resistance has been incorporated within R+
s,i j and Rs,i j . The

1
(Ts,i j − Tr,i )
R+
s,i j
1
= − (Tr, j − Ts,i j )
Rs,i j
+
= Q̇−
s,i j − Q̇s,i j ,

Q̇+
s,i j =

Cs,i j Ṫs,i j

(1)

+
where Q̇−
s,i j and Q̇s,i j denote the heat input and output for
separator i j, respectively.

Regarding the RH system, it is modeled by two thermal
restrictors Rw,i (V̇w,i ) separated by a thermal capacitance Cw,i
(that determines the heat stored in the thermal medium). In
.
particular, Rw,i = Rw,i (V̇w,i ) = 1/ρw cwV̇w,i and Cw,i = cw ρwVw ,
where ρw , cw and Vw are the density, specific heat capacity and volume of the water in the RH system. Due to the
capacitance of the thermal medium, a state variable Tw,i is
introduced that represents the temperature of the water at the
point of heat exchange. The radiator exchanges heat with
the thermal zone due to heat convection of resistance Rc,i .
Furthermore, the inlet water temperature in the RH system is
+
denoted by Tw,i
, while the outlet water temperature, denoted
−
Tw,i , is considered approximately equal to Tw,i (if we assume
uniform temperature distribution within the radiator device).
Thus, the heat transfer dynamics of the RH system are
Q̇w,i =

1
(Tw,i − Tr,i )
Rc,i

Cw,i Ṫw,i =

1
(T + − Tw,i ) − Q̇w,i ,
Rw,i (V̇w,i ) w,i

(2)

vector, respectively. Furthermore, we define
.
A
i (ui ) =

ar,i (V̇a,i ) row{a+
arw,i
rs,i j } j∈Ni
 col{a+

0
s,i j } j∈Ni diag{as,i j } j∈Ni
aw,i
0
awc,i (V̇w,i )

where Q̇w,i is the heat input to zone i through the radiator.
Regarding the HVAC system, we consider the standard
assumption that the fluid is incompressible, thus we are only
interested in the thermal part of the system [16]. A thermal
resistance, Ra,i (V̇a,i ), is introduced to represent a restrictor of
heat flow. It is given by Ra,i = Ra,i (V̇a,i ) = 1/ρa caV̇a,i where ρa
denotes the air density, ca denotes the air specific heat capacity at constant pressure, and V̇a,i is the air volume rate. The
outlet temperature of air in the HVAC system is considered
approximately equal to the zone temperature, while the in+
let air temperature is denoted by Ta,i
. Thus, the heat input
attributed to the HVAC system is
Q̇HVAC,i =

1
(T + − Tr,i ).
Ra,i (V̇a,i ) a,i

Finally, the heat-mass transfer dynamics of the thermal
zone i can be expressed as follows

∑ Q̇+s,i j + Q̇w,i + Q̇HVAC,i + Q̇ext,i .

(4)

j∈Ni

From equations (1), (2), (3) and (4), we can derive in a
straightforward manner the state-space representation of the
overall system, which can be written in a generalized form
as follows
ẋi (t) = Ai (ui (t))xi (t) + Ei (ui (t))di (t),

i ∈ I,

(5)

where

Tr,i
.
xi =  col{Ts,i j } j∈Ni  ,
Tw,i


+
Tw,i
+

Ta,i
.

di = 
 col{Tr, j } j∈N  ,


Q̇ext,i

.
ui =



where
.
ar,i (V̇a,i ) = − ∑ j∈Ni 1/Cr,i R+s,i j − 1/Cr,i Rc,i − 1/Cr,i Ra,i (V̇a,i ),
. 1
.
+
a+
ara,i (V̇a,i ) = 1/Cr,i Ra,i (V̇a,i ),
rs,i j = /Cr,i Rs,i j ,
.
. 1
−
a+ = 1/Cs,i j R+s,i j , a−
s,i j = /Cs,i j Rs,i j ,
(6)
. s,i j+
. 1
−
as,i j = −as,i j − as,i j , aww,i (V̇w,i ) = /Cw,i Rw,i (V̇w,i ),
. 1
. 1
. 1
aw,i = /Cw,i Rc,i , aext,i = /Cr,i , arw,i = /Cr,i Rc,i ,
.
awc,i (V̇w,i ) = −aw,i − aww,i (V̇w,i ).

(3)

The disturbance Q̇ext,i denotes heat inputs from external
sources. It will correspond to a vector of external heat flow
rates, attributed to solar radiation, human presence, etc.

Cr,i Ṫr,i =

.
E
i (ui ) =

0
ara,i (V̇a,i )
0
aext,i
−

0
0
diag{as,i j } j∈Ni 0 
aww,i (V̇w,i )
0
0
0

V̇w,i
V̇a,i


,

i

are the state vector, control-input vector and disturbance

2.2

Information structures
Different information structures will be considered depending on the sensors available. More specifically, for the
remainder of the paper, we consider the following cases:
• Full information (FI): All state, control and disturbance
variables can be measured except for the separator state
variables {Ts,i j } j .
• Medium Information (MI): All state, control and disturbance variables can be measured except for the separator state variables, {Ts,i j } j , and the water state variable
at the point of heat exchange, Tw,i .
• Limited information (LI): All state, control and disturbance variables can be measured except for the separator state variables {Ts,i j } j , the water state variable Tw,i
−
+
and the disturbance variables Tw,i
and Tw,i
.
The specifics can also be visualized in Table 1. The introduced alternative structures intend on demonstrating the
flexibility of the proposed derivation of nonlinear regressors,
since different residential buildings may not necessarily incorporate all sensors of the FI structure.
The case of limited information (LI) corresponds to the
case where only the flow rates are known in the RH and/or
the HVAC system, while the corresponding inlet/outlet temperatures of the thermal medium cannot be measured. Such
assumption may be reasonable for most today’s residential
buildings.
Alternative information structures may be considered,
such as the ones where Q̇ext,i (due to radiation or peoples’
presence) is not known, however such cases can be handled
through the above cases and with small modifications.

Tr,i

{Ts,i j } j∈Ni

Tw,i

V̇w,i

V̇a,i

+
Tw,i

+
Ta,i

{Tr, j } j∈Ni

Q̇ext,i

(FI)

X

X

X

X

X

X

X

X

(MI)

×

X

×

X

X

X

X

X

X

(LI)

×

X

×

×

X

X

×

×

X

X

Table 1: Information structures with measured variables.
2.3

The RH dynamics can be written as follows

System separation under FI

Under the full information structure (FI), a natural separation of the system dynamics can be introduced between
the zone and RH dynamics. Such separation will be used at
several points throughout the paper, and it is schematically
presented in Figure 2.

ẋw,i (t) = Aw,i (uw,i (t))xw,i (t) + Ew,i (uw,i (t))dw,i (t),
where
.
xw,i = Tw,i ,

dw,i
uw,i

dr,i
Tw,i

RH


. 
Aw,i (ui ) = awc,i (V̇w,i ) ,

ur,i
Fig. 2: System separation architecture under FI.

ẋr,i (t) = Ar,i (ur,i (t))xr,i (t) + Er,i (ur,i (t))dr,i (t)

(7)

.
where ur,i = V̇a,i is the control-input vector and

Tw,i
+

Ta,i
. 

dr,i = 
 col{Tr, j } j∈N  ,
i
Q̇ext,i


.
xr,i =

Tr,i
col{Ts,i j } j∈Ni


,

are the state vector and disturbance vector of the zone dynamics, respectively. Furthermore, we define
.
Ar,i (ui ) =



row{a+
rs,i j } j∈Ni

ar,i (V̇a,i )
,
col{a+
s,i j } j∈Ni diag{as,i j } j∈Ni




0
aext,i
. arw,i ara,i (V̇a,i )
Er,i =
.
0
0
diag{a−
s,i j } j∈Ni 0

.
dw,i =



+
Tw,i
Tr,i




. 
Ew,i = aww,i (V̇w,i ) aw,i .

3

In particular, the overall (heat-mass transfer) dynamics
for a single zone can be written as follows



.
uw,i = V̇w,i ,

are the state vector, control-input vector and disturbance vector of the RH system, respectively. Furthermore, we define
the following system matrices

Tr,i

Zone

(8)

Linear Regression Models (LRM)
We would like to focus on models for system identification which (to great extent) retain the physical interpretation of the parameters. To this end, an output-error model
structure (cf., [9]) is considered, where a one-to-one correspondence between the identification parameters and the
discrete-time model parameters can be established. In the
remainder of this section, we provide some background on
the Output-Error (OE) model structure, while we also discuss the implicit assumptions made when a linear structure
of the dynamics is assumed.

3.1

Background: Output-Error (OE) model
If we assume that the relation between input u and the
undisturbed output of the system w can be written as a linear difference equation, and that the disturbance consists of
white measurement noise v, then we obtain:
w(k) + f1 w(k − 1) + ... + fn f w(k − n f )
= b1 u(k − 1) + ... + bnb u(k − nb )

(9)

and the output is y(k) = w(k) + v(k), for some positive integers n f ≥ nb . The parameter vector to be determined is:
.
θ = (b1 , . . . , bnb , f1 , . . . , fn f ). Since w(k) is not observed, it
should be constructed from previous inputs and it should
carry an index θ.
The natural predictor (resulting from a maximum a pos-

teriori predictor1 ) is: ŷ(k|θ) = w(k, θ), and it is constructed
from past inputs only. If we define the vector:

.
ϕ(k, θ) = u(k − 1), . . . , u(k − nb ),

−w(k − 1, θ), . . . , −w(k − n f , θ) ,

may lead to significant identification errors. To see this, let
us consider the full information structure (FI). In a forthcoming section, we will show that the following approximation
of the RH part of the dynamics holds for sufficiently small
sampling period ε > 0:

(10)

then, the predictor can be written more compactly as:
ŷ(k|θ) = ϕ(k, θ)T θ, leading to a linear regression model
(LRM), where ϕ(k, θ) is called the regression vector. To simplify notation, in several cases we will write ϕ(k) instead of
ϕ(k, θ). Note that in the above expression, the w(k − j, θ),
j = 1, 2, ..., n f , are not observed, but using the above maximum a-posteriori predictor, they can be computed using
previous inputs as follows: w(k − j, θ) = ŷ(k − j|θ), j =
1, 2, ..., n f .
Such output-error model predictors will be used
throughout the paper. However, note that, depending on
the physics of the identified system, the regression vector
ϕ defined in (10) may be nonlinear with respect to the input/output lagged variables. In such cases, the resulting predictor will be referred to as a linear regression model with
nonlinear regressors (NRM).
3.2

Tw,i (k + 1) ≈
+
(1 + εawc,i (V̇w,i (k)))Tw,i (k) + εaww,i (V̇w,i (k))Tw,i
(k)+
εaw,i Tr,i (k),
plus higher-order terms of ε.
According to this (finite-response) approximation, a linear regression vector ϕ may describe well the evolution of
Tw,i as long as either the flow rate V̇w,i or the temperatures
+
Tw,i and Tw,i
are not varying significantly with time. However, variations in the flow rate V̇w,i and in the water temperature Tw,i may be large with time. Similar are the conclusions
when investigating the effect of the air flow rate V̇a,i in the
evolution of the zone temperature. The exact effect of these
nonlinear effects cannot be a-priori determined.

4

Linear Regression Models with Nonlinear Regressors
(NRM)
In this section, we will explore the formulation of OE
model structures when the regressors in ϕ may be nonlinear
functions of lagged input/output variables. Such investigation will be performed under the FI, MI and LI structure,
extending prior work of the same authors [1] to a larger set
of possible information structures (beyond the FI structure).
For the derivation, we will be using the following notation: Let q−1 {·} denote the one-step delay operator, i.e.,
q−1 {x(k)} = x(k − 1). Note that the delay operator is linear.
Let us also define the following operators:
.
Ps,i j = [(1 + εas,i j )q−1 ]
.
Pwc,i (V̇w,i ) = [(1 + εawc,i (V̇w,i ))q−1 ]
.
Pr,i (V̇a,i ) = [(1 + εar,i (V̇a,i ))q−1 ]

Discussion
Output-error model structures do not consider any disturbance terms in the process dynamics, instead they only
consider measurement noise. Hence, such model structures
are rather appropriate for providing a clear picture of the impact of the assumed (physical) model into the prediction error. On the other hand, when considering structures with process noise, such as ARMAX model structures (cf., [9, Section 4.2]), the process dynamics are perturbed by artificial
terms, which are not easily motivated by the physical insight.
Since the goal of this paper is to evaluate the impact of the
assumed (physical) dynamics, we consider OE models more
appropriate for evaluating predictions. This, however, does
not imply that the performance of an OE model structure is
necessarily better compared to an ARMAX model. The goal
of this paper is not to provide such a comparison.
Furthermore, identifiability (cf. [9, Section 4.6]) of the
considered OE model structures will be guaranteed by the
(inherent) controllability of the heat-mass transfer dynamics.

where ε > 0 defines the sampling period. Define also
.
Qs,i j = [1 − Ps,i j ],
.
Qwc,i (V̇w,i ) = [1 − Pwc,i (V̇w,i )],
.
Qr,i (V̇a,i ) = [1 − Pr,i (V̇a,i )].

3.3

For any operator P{·}, P−1 {·} will denote its inverse operator, i.e.,

Linear approximation & implicit assumptions
We wish to explore the utility of an OE model structure
using a linear regression vector ϕ into identifying the nonlinear system dynamics of Equation (5). Note that the original
dynamics (5) are bilinear in nature due to multiplications of
the flow rates (V̇w,i , V̇a,i ) with state or disturbance variables.
An investigation of the original dynamics (5) reveals that
when an output-error model is considered that admits a linear regression vector ϕ, we implicitly admit assumptions that

1 i.e., the one that maximizes the probability density function of the output given observations up to the previous time instant.

P−1 {P{x(k)}} = x(k).
Property 4.1. For each zone i ∈ I , the family of operators
{Qs,i j } j∈Ni are pairwise commutative, i.e.,
Qs,i j Qs,i j0 {x(k)} = Qs,i j0 Qs,i j {x(k)}
for any j, j0 ∈ I with j 6= j0 .

Proof. For each zone i ∈ I , and for any two neighboring
zones j, j0 ∈ Ni , we have the following:

T̂w,i (k|θw,i ) ≈ ϕw,i (k)T θw,i plus higher-order terms of ε, where
θw,i is a vector of unknown parameters and

Qs,i j Qs,i j0 {x(k)}


= [1 − Ps,i j ]  1 − Ps,i j0 {x(k)}
= [1 − Ps,i j ] x(k) − (1 + εas,i j0 )x(k − 1)
= x(k) − (1 + εas,i j0 )x(k − 1)−
(1 + εas,i j )x(k − 1) + (1 + εas,i j )(1 + εas,i j0 )x(k − 2).


T̂w,i (k − 1|θw,i )
.  V̇w,i (k − 1)T̂w,i (k − 1|θw,i ) 
.
ϕw,i (k) = 
+
 V̇w,i (k − 1)Tw,i
(k − 1) 
Tr,i (k − 1)

It is straightforward to check that the same expression also
occurs if we expand Qs,i j0 Qs,i j , due to the fact that (1+εas,i j )
commutes with (1 + εas,i j0 ). Thus, the conclusion follows. •



Proof. By a Taylor-series expansion of the RH dynamics, the
finite-step response of the water temperature Tw,i (k) is
Tw,i (k) ≈
(1 + εawc,i (V̇w,i (k − 1)))Tw,i (k − 1)+
+
εaww,i (V̇w,i (k − 1))Tw,i
(k − 1) + εaw,i Tr,i (k − 1),

Another two properties that will be handy in several
cases are the following:
Property 4.2. For each zone i ∈ I , we have

plus higher-order terms of ε. Equation (12) directly results
from the definition of the delay operator Pwc,i (V̇w,i ). From
Equation (12), a regression vector may be derived of the desired form. •

[Qs,i j Qwc,i (V̇w,i (k))] = [Qwc,i (V̇w,i (k))Qs,i j ]−


(1 + εas,i j )ε (1 − q−1 )awc,i (V̇w,i (k)) q−2 .
Proof. See Appendix A. •
Property 4.3. For some finite index set A and signal x(k),

h

∏ Qs,i j
j∈A

i
{x(k)} = x(k) +

|A|

∑ αm x(k − m),

(13)

(11)

m=1

for some constants α1 , ..., α|A| .
Proof. The conclusion follows in a straightforward manner
by induction and the definition of the delay operator Qs,i j . •

4.1.2 Zone dynamics
A similar approach to Proposition 4.1 for the derivation
of a nonlinear regression vector can also be employed for the
zone dynamics.
Proposition 4.2 (NRM for Zone under FI). For
sufficiently small sampling period ε > 0, the maximum
a-posteriori predictor of Tr,i (k) can be approximated by
T̂r,i (k|θr,i ) ≈ ϕr,i (k)T θr,i plus higher-order terms of ε, where
θr,i is a vector of unknown parameters and

|Ni |+1
col{T̂r,i (k − m|θr,i )}m=1


n
o


|Ni |+1

 col col{Tr, j (k − m)}m=2

j∈Ni 

Ni |+1 

.  col{V̇a,i (k − m)T̂r,i (k − m|θr,i )}|m=1
.
ϕr,i (k) = 


 col{V̇a,i (k − m)T + (k − m)}|Ni |+1 


a,i
m=1


|Ni |+1


col{T̂w,i (k − m|θw,i )}m=1


|Ni |+1
col{Q̇ext,i (k − m)}m=1


4.1

Full-Information (FI) structure
Using the natural decomposition of the dynamics under
the FI structure described in Section 2.3, in the following
subsections we provide a derivation (using physical insight)
of nonlinear regression vectors under an OE model structure
for the two subsystems of the overall dynamics.
4.1.1 RH dynamics
By isolating the RH dynamics, we first formulate a
(finite-response) prediction of the water temperature Tw,i , as
the following proposition describes.
Proposition 4.1 (NRM for RH under FI). For
sufficiently small sampling period ε > 0, the RH dynamics (8)
can be approximated by:
Tw,i (k) ≈ Pwc,i (V̇w,i (k − 1)){Tw,i (k)}+

+
εaww,i (V̇w,i (k − 1))Tw,i
(k − 1) + εaw,i Tr,i (k − 1), (12)

plus higher-order terms of ε. Furthermore, the maximum
a-posteriori predictor of Tw,i (k) can be approximated by

Proof. See Appendix B. •
4.2

Medium-Information (MI) structure
The case of MI structure is slightly different from the
FI structure, since the water temperature, Tw,i , at the point
of heat exchange cannot be directly measured, therefore the
dynamics of the RH and the zone cannot be separated as was
demonstrated in Figure 2. The following proposition provides the corresponding NRM predictor for the MI structure.
Proposition 4.3 (NRM for Zone under MI). For sufficiently small sampling period ε > 0, the maximum aposteriori predictor of Tr,i (k) can be approximated by

T̂r,i (k|θr,i ) ≈ ϕr,i (k)T θr,i plus higher-order terms of ε, where
θr,i is a vector of unknown parameters and
.
ϕr,i (k) =


|Ni |+2
col{T̂r,i (k − m|θr,i )}m=1
n
o
|Ni |+2
col col{Tr, j (k − m)}m=2











j∈Ni


N
+2
|
|
i


col{V̇a,i (k − m)T̂r,i (k − m|θr,i )}m=1




|Ni |+2
+


col{
V̇
(k
−
m)T
(k
−
m)}
a,i
a,i
m=1




N
+2
|
|
i
+

 col{V̇a,i (k − m)V̇w,i (k − m)Ta,i
(k − m)}m=2


(14)

N
+1
| i|
.

col{
V̇
(k
−
m
−
1)
T̂
(k
−
m|θ
)}
w,i
r,i
r,i
m=1




|Ni |+2


col{V̇w,i (k − m)T̂r,i (k − m|θr,i )}m=2


|Ni |+2 

 col{V̇w,i (k − m)V̇a,i (k − m)T̂r,i (k − m|θr,i )}m=2 


|Ni |+2


+
col{V̇w,i (k − m)Tw,i
(k − m)}m=2




|Ni |+2


col{Q̇ext,i (k − m)}m=1


|Ni |+2
col{V̇w,i (k − m)Q̇ext,i (k − m)}m=2
Proof. The proof follows similar reasoning with the proof of
Proposition 4.2 and it is presented in Appendix C. •
4.3

Limited-Information (LI) structure
Similarly to the derivation of the NRM for the FI and
MI structures, we may also derive the corresponding NRM
for the case of the LI structure. Note that the difference of
such information structure (compared to MI) is the fact that
+
+
the inlet temperatures of the thermal mediums (Tw,i
and Ta,i
)
cannot be measured. In such case, the only valid assumption
is the fact that the inlet temperature of the thermal mediums
remain constant with time. In this case, Proposition 4.3 con+
+
tinue to hold, when Tw,i
and Ta,i
are replaced with unity in
(14) (i.e., they become part of the parameters in θr,i ).
5

Performance Evaluation
In this section, we provide a comparison between standard linear regression models and the nonlinear ones derived
in Section 4. The comparison will be performed both with
respect to a) identification error, and b) prediction performance within a standard MPC implementation for climate
control in a residential building. The second part is considered the most representative of the performance of an identification model since even small prediction errors might lead
to significant performance degradation.
5.1

Simulation platform
In order to compare the proposed nonlinear regression
vectors with standard linear regressors, we used the EnergyPlus (V7-2-0) building simulator developed by the U.S.
Department of Energy [13]. The BCVTB simulation tool
has also been used for allowing data collection and also climate control developed in MATLAB to be implemented during run-time. A typical (three-storey) residential building in

Linz, Austria, was accurately modeled and simulated with
the EnergyPlus environment to allow for collecting data from
a realistic residential environment.
5.2

Data generation
The data collection for system identification was performed under normal operating conditions of the heating system during the winter months (October - April) under the
weather conditions of Linz, Austria. To replicate normal operating conditions, a standard hysteresis controller was employed, according to which the water flow V̇w,i (k) is updated
regularly at time instances tk = kTsam where Tsam = 1/12h.
The control law is as follows:
.
V̇w,i (k) =


V̇w,max






0

if pi (k) > 0 and
(Tr,i (k) < Tset − ∆T or
(Tr,i (k) ≥ Tset − ∆T and Tr,i (k − 1) ≤ Tr,i (k))),
else

where ∆T determines a small temperature range about the
desired (set) temperature, Tset , in which the control flow
maintains its value and V̇w,max denotes the maximum water
flow set to 0.0787kg/sec. Furthermore, pi (k) ∈ {0, 1} indicates whether people are present in thermal zone i. In other
words, the hysteresis controller applies only if someone is
present in the thermal zone i, which can be easily determined
through a motion detector. Furthermore, the inlet water tem+
perature Tw,i
is determined by the following heating curve:

+
Tw,i

.
=

(
ρ0 + ρ1 · (Tset − Tout )ζ , if Tset > Tout
ρ0 ,
else.

(15)

where ρ0 , ρ1 and ζ are positive constants. For the simulation,
we set ∆T = 0.1, ρ0 = 29.30, ρ1 = 0.80 and ζ = 0.97. The
set temperature was set equal to Tset = 21oC.
There exists a natural ventilation system that operates
autonomously with an intermittent flow pattern V̇a,i . Thus,
the only control parameter is V̇w,i . All parameters mentioned
in Table 1 can be measured, except for {Ts,i j } j , which allows
for evaluating all considered structures (FI, MI and LI).
5.3

Recursive system identification
To utilize the linear and nonlinear regression vectors for
system identification of heat-mass transfer dynamics of a residential building, we use an OE model structure as described
in Section 4, while we employ a regularized recursive leastsquares implementation (cf., [17, Section 12.3]) for training its parameters. The reason for implementing a recursive identification procedure is primarily due to the fact that
the size of the available data is quite large, which makes the
use of standard least-squares approaches practically infeasible. Besides, predictions for the zone temperature will be

needed continuously during run-time, demanding for more
efficient computational schemes. Furthermore, a recursive
least squares implementation allows for an adaptive response
to more recent data, thus capturing more recent effects.

the form (17) the following condition is satisfied:2
#2

"

E ∆Gi0 (q)yi (k) + ∑ ∆Gil (q)zl (k)
l∈L

5.4

Identification experiments
The experiments performed in this study will involve a)
a standard linear regression model, with b) the derived NRM
of Proposition 4.3 (which corresponds to the MI structure).
Our intention is to evaluate the benefit of the more accurate
physical representation of the models derived in Section 4.
In both regression vectors (linear or nonlinear) we consider
|Ni | = 1, i.e., we implicitly assume that there is a single
neighboring zone. This can be justified by the fact that the
building can be considered as a single thermal zone, since
the same heating schedule is employed in all rooms.
In particular, the linear regression model (LRM) implemented corresponds to the following structure:

|Ni |+2
col{
T̂
(k
−
m|θ
)}
r,i
r,i m=1
 n

o


|Ni |+2
 col col{Tr, j (k − m)}m=1


j∈Ni 


|Ni |+2


col{V̇a,i (k − m)}m=1


. 

|Ni |+2
+
ϕr,i (k) = 
,
col{Ta,i (k − m)}m=1




|Ni |+2


col{V̇w,i (k − m)}m=1




N
+2
|
|
 col{T + (k − m|θr,i )} i

w,i
m=1


|Ni |+2
col{Q̇ext,i (k − m)}m=1


(16)

i.e., it only considers delayed samples of the state and the
available input and disturbance parameters. The NRM model
implemented corresponds to the one of Equation (14).

5.5

Persistence of Excitation
Following the analysis of Section 4, for any thermal
.
zone i ∈ I , the output of the zone yi (k) = Tr,i (k) can briefly
be expressed in the following form:
(17)

l∈L

for some linear transfer functions Gi0 (q), Gil (q), l ∈ L and
some index set L ⊂ N. The terms zil (k) represent time sequences of input, disturbance or output signals, or products
of such terms. For example, if we consider the MI model
presented in (14), the terms zil may correspond to a disturbance, such as Tr, j , j ∈ Ni , or to a product of a control input
and the output, such as V̇w,i (k)Tr,i (k).
A model of the generic form of equation (17) is uniquely
determined by its corresponding transfer functions Gi0 (q),
Gil (q), l ∈ L. Thus, in order for the data to be informative,
the data should be such that, for any two different models of

(18)

where ∆Gi0 (q), ∆Gil (q) denote the corresponding differences between the transfer functions of the two different
.
models. Moreover, E = limN→∞ 1/N ∑Nk=1 E[·] and E[·] denotes the expectation operator. Since the models are assumed
different, either ∆Gi0 (q) or ∆Gil (q), for some l ∈ L, should
be nonzero. Thus, if ∆Gil (q) 6= 0, for some l ∈ L, the persistence of excitation 3 of zil (k) will guarantee (18).
Since data are collected through a closed-loop experiment, special attention is needed for terms zil (k) which involve the product of the water flow rate V̇w,i (k) with the output yi (k) of the system, such as the term V̇w,i (k)yi (k) in (14).
Note, however, that this term is nonlinear in both yi (k) and
the presence/occupancy indicator, since the hysteresis controller a) applies only when someone is present in the zone
and b) is a nonlinear function of the output yi (k). Thus, it
is sufficient for either the input signal or the occupancy indicator to be persistently exciting in order for the experiment
to be informative. In Figure 3, we provide the positive spectrum generated through the discrete Fast Fourier Transform
(DFFT) of the disturbance and control signals.
For either the simplified model of (16) or the more detailed model of (14), note that any transfer function (from
an input/disturbance to the zone temperature) requires the
identification of at most 2( Ni + 2) parameters. Thus, in
order for an experiment to be informative enough, it is sufficient that the considered inputs/disturbances are persistently exciting of order 2( Ni + 2) or higher. In the case of
Ni = 1 (considered in this experiment), it suffices that any
input/disturbance signal is persistently exciting of order 6 (or
5 in case one of the frequencies is at zero, as it is the case for
all signals in Figure 3). As demonstrated in Figure 3, this
condition is satisfied by the input/disturbance signals, since
for all of them the positive spectrum is non-zero in at least 3
distinct frequencies.
5.6

yi (k) = Gi0 (q)yi (k) + ∑ Gil (q)zil (k),

6= 0,

Identification error comparison
In Figure 4, we demonstrate the resulting identification error under the LRM of (16) and the NRM of (14) for
Ni = 1. Note that the NRM achieves a smaller identification error of about 10%. This reduced error is observed later
in time (due to the larger training time required by the larger
number of terms used in the nonlinear regression vector).
The data used for training correspond to data collected between October and April from the simulated building, however the same data have been reused several times for better
fitting. Thus, the time axis in Figure 4 does not correspond
to the actual simulation time, but it corresponds to the accumulated time index of the reused data.

2 cf.,
3 cf.,

discussion in [9, Section 13.4].
[9, Definition 13.2]

Presence

Running Average Squared Error (oC2 )
0.4

Power

400

0.3

200
0

LRM
NRM

0

0.1

0.2

0.3
0.4
Water Flow

0.5

0.6

0.2
0.1

0

1

0.5

Power

400

Power
Power

1
0.8
0.6
0.4
0.2
0

1
0.8
0.6
0.4
0.2
0

Power

2

0
·105

0
·104

0

0.1
0.2
0.3
Water Inlet Temperature

1

0.5

1.5

2

2.5

Outdoor Temperature

5 · 10−2

·107

0.1

0.4

0.15

3
·10−2

0.2

min

·105

Nhor n

α∑

o
2
p̂i (k) T̂r,i (k) − Tset (k) /Nhor +

k=0
Nhor −1 n

Solar Radiation

∑

k=0

s.t.

1
0.5
0

5 · 10−2

0.1

0.15

0.2

Frequency (1/h)

Fig. 3: Positive Spectrum of Input and Disturbance Signals.

5.7

3.5

weather conditions may dramatically change within a few
days). This may potentially result in a significant degradation of the prediction performance when the prediction has
not taken into account the effect of the nonlinearities present
in the dynamics. The effect of these nonlinearities may not
be known a-priori and a thorough investigation is necessary
to understand their potential economic impact.
To this end, we designed a standard MPC for the RH
system of the main living area of the residential building.
The goal is to evaluate the performance of the previously designed (nonlinear) regression vectors for the zone temperature prediction compared to standard (black-box) linear regression vectors. The structure of the MPC is rather simple
and addresses the following optimization problem.



o
−
+
βTsam Tw,i
(k) − T̂w,i
(k) + γTsamV̇w,i (k)
(19a)

1.5

0

3

Fig. 4: Identification performance for the zone temperature.

200
0

2
1.5
2.5
time (hours)

Prediction performance comparison
Small identification error may not necessarily guarantee
a good prediction performance, especially when predictions
are required for several hours ahead. From the one hand,
training for system identification is based upon the one-step
ahead innovations, which does not necessarily guarantee a
similar accurate prediction for several hours ahead. On the
other hand, training might have been performed under operating conditions that are (possibly) different from the operating conditions under which predictions are requested (e.g.,

var.

T

(19b)

−
T̂w,i
(k) ≈ T̂w,i (k) ≈ ϕw,i (k) · θw,i
+
Tw,i (k) ∈ {40oC, 45oC},

(19c)

T̂r,i (k) ≈ ϕr,i (k) θr,i

V̇w,i (k) ∈ {0, V̇w,max },

(19d)
(19e)

k = 0, 1, 2, ..., Nhor − 1,

Note that the first part of the objective function (19a)
corresponds to a comfort cost. It measures the average
squared difference of the zone temperature from the desired
(or set) temperature entered by the user at time k. The set
temperature was set equal to 21oC throughout the optimization horizon. The variable p̂i (k) ∈ {0, 1} holds our estimates
on whether people are present in zone i at time instance k.
The second part of the objective function (19a) corresponds to the heating cost, while the third part corresponds
to the pump-electricity cost. The nonnegative parameters
β, γ were previously identified for the heating system of
the simulated building and take values: β = 0.3333kW /oCh,

γ = 0.5278 · 103 kW sec/hm3 . The non-negative constant α is
introduced to allow for adjusting the importance of the comfort cost compared to the energy cost. A large value has been
assigned equal to 106 to enforce high comfort.
The p̂i (k), k = 1, 2, ...,, as well as the outdoor temperature, Tout , and the solar gain, Q̇ext (k), are assumed given
(i.e., predicted with perfect accuracy). This assumption is
essential in order to evaluate precisely the impact of our temperature predictions (19b) in the performance of the derived
optimal controller.
The sampling period was set to Tsam = 1/12h, the optimization period was set to Topt = 1h, and the optimization
horizon was set to Thor = 5h. This implies that Nhor = 5 ·12 =
60. Furthermore, the control variables are the inlet water
temperature which assumes two values (40oC and 45oC) and
the water flow which assumes only two values, the minimum
0kg/sec and the maximum V̇w,max = 0.0787kg/sec.
For the prediction model of (19b), we used either the
LRM of (16) or the NRM of (14), both trained offline using data collected during normal operation of a hysteresis
controller (as demonstrated in detail in Section 5.2). For
the prediction model of the outlet water temperature used in
the computation of the cost function, we used the prediction
model derived in (13). As in the case of Section 5.6 where the
LRM and the NRM were compared, we evaluated the value
of the cost function under these two alternatives. The performances of the two models with respect to the comfort cost
(which corresponds to the 1st term of the objective (19a))
are presented in Figure 5. The performances of the two models with respect to the energy spent (i.e., heating and pump
electricity cost), which correspond to the 2nd and 3rd term
of the objective (19a), are presented in Figure 6. Note that
the NRM achieves lower comfort cost using less amount of
energy which is an indication of small prediction errors.

3
2.5
2
1.5
1
0.5
0

Comfort Cost |Tr,i (k) − Tset | (oC)
LRM
NRM

0

200

400

600 800
Time (h)

1,000 1,200

Fig. 5: Running-average comfort cost under LRM and NRM.

5.8

Discussion
Note that the NRM achieves an improvement to both
prediction accuracy and cost performance. The fact that the
LRM may not provide a similarly good performance as the

Heating and Pump Electricity Cost (kW h/Tsam )
0.2

0.1

0

LRM
NRM
0

200

400

600 800
Time (h)

1,000 1,200

Fig. 6: Running-average heating and pump electricity cost
under LRM and NRM.
NRM could be attributed to the fact that the operating conditions for testing the MPC controller might not be identical with the training conditions. Note, for example, that
the training phase was performed under the heating-curve
pattern of (15), while the MPC employs only two alternative water inlet temperatures. This variation in the inlet water temperature may have exhibited different patterns in the
nonlinear parts of the dynamics, thus affecting the prediction
performance.
The simulation study presented in this paper considered
a regression vector generated with Ni = 1, which corresponds to the case of a single neighboring zone (e.g., the
outdoor temperature). The proposed regression models are
also applicable for higher-order models (i.e., when more than
one neighboring zones are present), however it might be necessary that some of the input signals are designed to include
larger number of frequencies. Note, however, that the input
signals that can be designed (e.g., medium flow rates and inlet temperatures) appear as products with other disturbance
variables in the derived regression models (e.g., equation
(14)). Thus, persistence of excitation for such products of
non-correlated signals is an easier task. On the other hand,
some of the disturbance signals (e.g., the solar radiation or
the outdoor temperature) cannot be designed. Of course, the
prediction model will operate at conditions that are similar to
which it was identified. However, the impact of such disturbance signals when they are not informative enough needs to
be further investigated.

6

Conclusions
We derived analytically regression models admitting
nonlinear regressors specifically tailored for system identification and prediction for thermal dynamics in buildings. The
proposed models were compared with standard (black-box)
linear regression models derived through the OE structure,
and an improvement was observed with respect to both the
resulting prediction performance as well as the energy cost.
The greatest advantage of the proposed identification scheme
relies on the fact that it provides a richer and more accurate
representation of the underlying physical phenomena, contrary to standard black-box identification schemes.

References
[1] Chasparis, G., and Natschläger, T., 2014. “Nonlinear system identification of thermal dynamics in buildings”. In Control Conference (ECC), 2014 European,
pp. 1649–1654.
[2] Nghiem, T., and Pappas, G., 2011. “Receding-horizon
supervisory control of green buildings”. In Proc. of the
2011 American Control Conference, pp. 4416–4421.
[3] Oldewurtel, F., Parisio, A., Jones, C., Morari, M., Gyalistras, D., Gwerder, M., Stauch, V., Lehmann, B., and
Morari, M., 2012. “Use of model predictive control and
weather forecasts for energy efficient building climate
control”. Energy and Buildings, 45, pp. 15–27.
[4] Nghiem, T., Pappas, G., and Mangharam, R., 2013.
“Event-based green scheduling of radiant systems in
buildings”. In Proc. of 2013 American Control Conference (ACC), pp. 455–460.
[5] Touretzky, C., and Baldea, M., 2013. “Model reduction
and nonlinear MPC for energy management in buildings”. In Proc. of 2013 American Control Conference
(ACC), pp. 455–460.
[6] Coogan, S., Ratliff, L., Calderone, D., Tomlin, C., and
Sastry, S., 2013. “Energy management via pricing in
LQ dynamic games”. In Proc. of 2013 American Control Conference (ACC), pp. 443–448.
[7] Rasmussen, B., Alleyne, A., and Musser, A., 2005.
“Model-driven system identification of transcritical vapor compression systems”. IEEE Transactions on Control Systems Technology, 13(3), pp. 444–451.
[8] Maasoumy, M., Razmara, M., Shahbakhti, M., and Vincentelli, A. S., 2014. “Handling model uncertainty in
model predictive control for energy efficient buildings”.
Energy and Buildings, 77, pp. 377 – 392.
[9] Ljung, L., 1999. System Identification: Theory for the
User, 2nd ed. Prentice Hall Ptr, Upper Saddle River,
NJ.
[10] Yiu, J.-M., and Wang, S., 2007. “Multiple ARMAX
modeling scheme for forecasting air conditioning system performance”. Energy Conversion and Management, 48, pp. 2276–2285.
[11] Scotton, F., Huang, L., Ahmadi, S., and Wahlberg, B.,
2013. “Physics-based modeling and identification for
HVAC systems”. In Proc. of 2013 European Control
Conference (ECC), pp. 1404–1409.
[12] Malisani, P., Chaplais, F., Petit, N., and Feldmann, D.,
2010. “Thermal building model identification using
time-scaled identification models”. In Proc. of 49th
IEEE Conference on Decision and Control, pp. 308–
315.
[13] EnergyPlus. EnergyPlus Energy Simulation software
Version: 7-2-0.
[14] Karnopp, D., Margolis, D., and Rosenberg, R., 2012.
System Dynamics: Modeling, Simulation and Control
of Mechatronic Systems, 5th ed. John Wiley & Sons,
Inc, Hoboken, NJ.
[15] Thirumaleshwar, M., 2009. Fundamentals of Heat &
Mass Transfer. Dorling Kindersley (India) Pvt. Ltd,
New Dehli, India.

[16] Karnopp, D., 1978. “Pseudo bond graphs for thermal
energy transport”. ASME Journal of Dynamic Systems,
Measurement and Control, 100, pp. 165–169.
[17] Sayed, A., 2003. Fundamentals of Adaptive Filtering.
John Wiley & Sons, Inc., New Jersey.

A

Proof of Property 4.2
We can write:
[Qs,i j Qwc,i
 (V̇w,i (k))] {x(k)}
= Qs,i j Qwc,i (V̇w,i (k)) {x(k)}
= Qs,i j x(k) − (1 + εawc,i (V̇w,i (k)))x(k − 1)
= x(k) − (1 + εawc,i (V̇w,i (k)))x(k − 1)−
(1 + εas,i j )x(k − 1) + (1 + εas,i j )x(k − 2)+
(1 + εas,i j )εawc,i (V̇w,i (k − 1))x(k − 2).

Similarly, we have:
[Qwc,i (V̇w,i (k))Qs,i j ] {x(k)}
= x(k) − (1 + εawc,i (V̇w,i (k)))x(k − 1)−
(1 + εas,i j )x(k − 1) + (1 + εas,i j )x(k − 2)+
(1 + εas,i j )εawc,i (V̇w,i (k))x(k − 2).
Thus, we have:
[Qs,i j Qwc,i (V̇w,i (k))] {x(k)} =
[Qwc,i (V̇w,i (k))Qs,i j ] {x(k)} −

(1 + εas,i j )ε awc,i (V̇w,i (k)) − awc,i (V̇w,i (k − 1)) x(k − 2).
Finally, we may write:
[Qs,i j Qwc,i (V̇w,i (k))] =
[Qwc,i (V̇w,i (k))Qs,i j ]−


(1 + εas,i j )ε (1 − q−1 )awc,i (V̇w,i (k)) q−2 .

(20)

B

Proof of Proposition 4.2
A Taylor-series expansion of the continuous-time zone
dynamics leads to the following approximation:
Tr,i (k) ≈
Pr,i (V̇a,i (k − 1)){Tr,i (k)} + ε

∑ a+rs,i j Ts,i j (k − 1)+

j∈Ni

+
εarw,i Tw,i (k − 1) + εara,i (V̇a,i (k − 1))Ta,i
(k − 1)+
εaext,i Q̇ext,i (k − 1),
(21)

plus higher-order terms of ε. Furthermore, note that the finite
response of the separator-dynamics can be approximated by:
n
o
−
+
Ts,i j (k) ≈ Q−1
(22)
s,i j εas,i j Tr,i (k − 1) + εas,i j Tr, j (k − 1)

plus higher-order terms of ε. Given that the family of operators Qs,i j , j ∈ Ni are pairwise commutative by Property 4.1,
if we replace Ts,i j into (21) and we multiply both sides of the
above expression by the composite operator Qs,i j1 · · · Qs,i j|N | ,
i
briefly denoted by ∏ j∈Ni Qs,i j {·}, we have:
Tr,i
h (k) ≈
1−

∏ Qs,i j

j∈Ni

h

Qs,i j

∏

j∈Ni
2

ε
ε

i
Pr,i (V̇a,i (k − 1)){Tr,i (k)} +

∑

a+
rs,i j

∑

a+
rs,i j

j∈Ni
2

j∈Ni

εarw,i

h

∏

j∈Ni

h
ε

in
o
Tr,i (k) +

∏

j∈Ni

εaext,i

h

in
o
∏ Qs,i` a+s,i j Tr,i (k − 2) +

h `6= j

in
o
−
Q
a
T
(k
−
2)
+
r,
j
s,i`
∏
s,i j

`6= j

in
o
Qs,i j Tw,i (k − 1) +

in
o
+
Qs,i j ara,i (V̇a,i (k − 1))Ta,i
(k − 1) +

h

∏ Qs,i j

j∈Ni

i
Q̇ext,i (k − 1) ,

plus higher-order terms of ε. The conclusion follows directly
by using Property 4.3 and expanding the terms of the above
expression.


εaext,i Qwc,i (V̇w,i (k − 2)) Q̇ext,i (k − 1)
plus higher-order terms of ε. Applying to both sides of the
above expression the composite operator Qs,i j1 · · · Qs,i j|N | ,
i
briefly denoted by ∏ j∈Ni Qs,i j {·}, and making use of Property 4.1, we have:
n
o
Tr,i (k) ≈

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

j∈Ni

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

j∈Ni
2

ε


Pr,i (V̇a,i (k − 1)){Tr,i (k)} +

∑ a+rs,i j ∏

Qs,i` [Qs,i j Qwc,i (V̇w,i (k − 2))]Q−1
s,i j

∑ a+rs,i j ∏

−1
Qs,i` [Qs,i j Qwc,i (V̇w,i (k − 2))]Qs,i
j

`∈Ni \ j
n j∈Ni
o
+
as,i j Tr,i (k − 2) +

ε2

`∈Ni \ j
n j∈Ni
o
−
as,i j Tr, j (k − 2) +
n
o
+
ε2 arw,i ∏ Qs,i j aww,i (V̇w,i (k − 2))Tw,i
(k − 2) +
j∈Ni

2

ε arw,i aw,i

∏ Qs,i j {Tr,i (k − 2)} +

j∈Ni

ε

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

j∈Ni

n
o
+
ara,i (V̇a,i (k − 1))Ta,i
(k − 1) +
εaext,i

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

j∈Ni

C

Proof of Proposition 4.3
Under the MI structure, the finite-response approximation of (21) continues to hold for sufficiently small step-size
sampling interval ε. Note that the separator dynamics can
still be approximated according to Equation (22). Furthermore, a Taylor-series expansion of the water dynamics leads
to the following finite-response approximation:
Tw,i (k) ≈ Q−1
wc,i (V̇w,i (k − 1))
n
o
+
εaww,i (V̇w,i (k − 1))Tw,i
(k − 1) + εaw,i Tr,i (k − 1)

plus higher-order terms of ε. According to Property 4.2, we
may perform the substitution:
[Qs,i j Qwc,i (V̇w,i (k − 2))] =
[Qwc,i (V̇w,i (k − 2))Qs,i j ]−


(1 + εas,i j )ε (1 − q−1 )awc,i (V̇w,i (k − 2)) q−2 .
Thus, we may write:

∏

plus higher-order terms of ε. If we replace Ts,i j and Tw,i into
(21) and we apply to both sides of the resulting expression
the operator Qwc,i (V̇w,i (k − 2)), we get
n
o
Qwc,i (V̇w,i (k − 2)) Tr,i (k) ≈

Qwc,i (V̇w,i (k − 2)) Pr,i (V̇a,i (k − 1)){T
n r,i (k)} + o
−1
+
2
ε ∑ ars,i j Qwc,i (V̇w,i (k − 2))Qs,i j a+
s,i j Tr,i (k − 2) +
j∈Ni

ε

2

∑

j∈Ni

n
o
−1
−
a+
Q
(
V̇
(k
−
2))Q
a
T
(k
−
2)
+
r,
j
wc,i
w,i
rs,i j
s,i j
s,i j

+
ε2 arw,i aww,i (V̇w,i (k − 2))Tw,i
(k − 2)+
2
ε arw,i aw,i Tr,i (k − 2)+
n

o
+
εQwc,i (V̇w,i (k − 2)) ara,i (V̇a,i (k − 1))Ta,i
(k − 1) +


Q̇ext,i (k − 1)

j∈Ni
2

ε

n
o
Qs,i j Qwc,i (V̇w,i (k − 2))Qr,i (V̇a,i (k − 1)) Tr,i (k) ≈

∑ a+rs,i j ∏

`∈Ni \ j

j∈Ni

Qs,i` Qwc,i (V̇w,i (k − 2))

n
o
a+
s,i j Tr,i (k − 2) +

ε2

∑ a+rs,i j ∏

`∈Ni \ j

j∈Ni

Qs,i` Qwc,i (V̇w,i (k − 2))

n
o
a−
s,i j Tr, j (k − 2) +
n
o
+
ε2 arw,i ∏ Qs,i j aww,i (V̇w,i (k − 2))Tw,i
(k − 2) +
j∈Ni

ε2 arw,i aw,i

∏ Qs,i j {Tr,i (k − 2)} +

j∈Ni

ε

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

j∈Ni

n
o
+
ara,i (V̇a,i (k − 1))Ta,i
(k − 1) +
εaext,i

∏

j∈Ni

|Ni |+2


Qs,i j Qwc,i (V̇w,i (k − 2)) Q̇ext,i (k − 1) , (23)

Qwc,i (V̇w,i (k − 2))Qr,i (V̇a,i (k − 1)){Tr,i (k)} =
Tr,i (k) − 2Tr,i (k − 1) + Tr,i (k − 2)−
εar,i (V̇a,i (k − 1))Tr,i (k − 1)−
εawc,i (V̇w,i (k − 2))Tr,i (k − 1)+
εawc,i (V̇w,i (k − 2))Tr,i (k − 2)+
εar,i (V̇a,i (k − 2))Tr,i (k − 2)+
ε2 awc,i (V̇w,i (k − 2))ar,i (V̇a,i (k − 2))Tr,i (k − 2).

∏

n
o
Qs,i` Qwc,i (V̇w,i (k − 2)) a+
T
(k
−
2)
r,i
s,i j

`∈Ni \ j

≈

Note that

α5,mV̇w,i (k − m)V̇a,i (k − m)Tr,i (k − m),

m=2

plus higher-order terms of ε. In the following, we will approximate the terms of the above expression, by neglecting
terms of order of ε3 or higher.

|Ni |+2

∑

m=2

∏

j∈Ni

=

|Ni |+2

∑

(4)

+
α1,mV̇w,i (k − m)Tw,i
(k − m),

∏ Qs,i j

j∈Ni

The following also hold:
=

|Ni |+2

∑

m=2

n
o
Qwc,i (V̇w,i (k − 2)) a+
T
(k
−
2)
r,i
s,i j

+
=a+
s,i j Tr,i (k − 2) − (1 + εawc,i (V̇w,i (k − 2)))as,i j Tr,i (k − 3)

n
o
Tr,i+ (k − 2)

(5)

α1,m Tr,i+ (k − m),

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

+
≈a+
s,i j Tr,i (k − 2) − as,i j Tr,i (k − 3)

n
o
+
Qwc,i (V̇w,i (k − 2)) ara,i (V̇a,i (k − 1)Ta,i
(k − 1)

(2)

α1,m Tr,i (k − m),

n
o
+
Qs,i j aww,i (V̇w,i (k − 2))Tw,i
(k − 2)

m=2

(24)

(1)

∑

j∈Ni

=

|Ni |+2

∑

m=1

|Ni |+2

+
=ara,i (V̇a,i (k − 1))Ta,i
(k − 1)−

∑

(5)

+
α1,mV̇a,i (k − m)Ta,i
(k − m)+
(5)

+
α2,mV̇w,i (k − m)V̇a,i (k − m)Ta,i
(k − m),

+
(1 + εawc,i (V̇w,i (k − 2)))ara,i (V̇a,i (k − 2))Ta,i
(k − 2)

m=2


Qwc,i (V̇w,i (k − 2)) Q̇ext,i (k − 1)

∏ Qs,i j Qwc,i (V̇w,i (k − 2))

j∈Ni

=Q̇ext,i (k − 1) − (1 + εawc,i (V̇w,i (k − 2)))Q̇ext,i (k − 2)
where the first expression ignores terms of ε since, in (23),
it multiplies an expression of order of ε2 . Using the above
approximations and Property 4.3, the terms of (23) can be
approximated as follows:
n
o
Tr,i (k)

∏ Qs,i j Qwc,i (V̇w,i (k − 2))Qr,i (V̇a,i (k − 1))

j∈Ni

=Tr,i (k) +
|Ni |+2

∑

m=1

|Ni |+1

∑

m=1

|Ni |+2

∑

m=2

|Ni |+2

∑

m=1

(1)

α1,m Tr,i (k − m)+

(1)

α2,mV̇a,i (k − m)Tr,i (k − m)+
(1)

α3,mV̇w,i (k − 1 − m)Tr,i (k − m)+
(1)

α4,mV̇w,i (k − m)Tr,i (k − m)+

=

|Ni |+2

∑

m=1

|Ni |+2

∑

m=2

n
o
+
ara,i (V̇a,i (k − 1))Ta,i
(k − 1)


Q̇ext,i (k − 1)

(6)

α1,m Q̇ext,i (k − m)+
(6)

α2,mV̇w,i (k − m)Q̇ext,i (k − m).
(×)

for some constant parameters α×,m ∈ R.
It is straightforward to check that using the above approximations, Equation (23) can be written as a linear regression with a nonlinear regression vector of the form depicted
in Equation (14).

