THE POWER OF DIVERSITY: DATA-DRIVEN ROBUST PREDICTIVE CONTROL
FOR ENERGY EFFICIENT BUILDINGS AND DISTRICTS

arXiv:1607.05441v1 [math.OC] 19 Jul 2016

GEORGIOS DARIVIANAKIS, ANGELOS GEORGHIOU, ROY S. SMITH AND JOHN LYGEROS

Abstract. The cooperative energy management of aggregated buildings has recently received a great deal
of interest due to substantial potential energy savings. These gains are mainly obtained in two ways: (i)
Exploiting the load shifting capabilities of the cooperative buildings; (ii) Utilizing the expensive but energy
efficient equipment that is commonly shared by the building community (e.g., heat pumps, batteries and
photovoltaics). Several deterministic and stochastic control schemes that strive to realize these savings, have
been proposed in the literature. A common difficulty with all these methods is integrating knowledge about
the disturbances affecting the system. In this context, the underlying disturbance distributions are often
poorly characterized based on historical data. In this paper, we address this issue by exploiting the historical
data to construct families of distributions which contain these underlying distributions with high confidence.
We then employ tools from data-driven robust optimization to formulate a multistage stochastic optimization
problem which can be approximated by a finite-dimensional linear program. The proposed method is suitable
for tackling large scale systems since its complexity grows polynomially with respect to the system variables.
We demonstrate its efficacy in a numerical study, in which it is shown to outperform, in terms of energy cost
savings and constraint violations, established solution techniques from the literature. We conclude this study
by showing the significant energy gains that are obtained by cooperatively managing a collection of buildings
with heterogeneous characteristics.

1. Introduction
Approximately 20-40% of the total energy consumption in the developed countries is attributed to the
building sector, an amount that often exceeds even the industrial and transportation sectors [1]. Concerns
about the growing environmental impact of building energy consumption, led the EU and the US government
to set the target of a net zero-energy for 50% of their commercial buildings by 2040 [2]. In this context,
active building energy management has attracted considerable attention with substantial efforts to be devoted to developing sophisticated control schemes that are capable of reducing the buildings energy impact
while ensuring comfortable conditions for the building users [3–6]. Nevertheless, the opportunities for large
savings within individual buildings can be limited, and depend on the specific building actuation systems and
construction characteristics [7–9].
Further savings can be envisaged by cooperatively managing the aggregated energy demands of a collection
of buildings in a district via an energy hub. The energy hub is a conceptual entity that provides the interface
between the building community and the power grid by utilizing shared energy generation, conversion and
storage equipment (e.g., heat pumps, batteries and photovoltaics) [10]. In this setting, significant energy
gains can be obtained by exploiting the diversity of the available energy sources, equipment and building
characteristics. However, the main body of the literature separates the optimal control of the energy hub
Date: July 20, 2016.
This research was partially funded by CTI within the SCCER FEEB&D, the Swiss National Science Foundation under the
project IMES and the European Commission under the project Local4Global.
The authors are with the Automatic Control Laboratory, Department of Electrical Engineering and Information Technology,
ETH Zurich, 8092 Zurich, Switzerland. {gdarivia, angelosg, rsmith, lygeros}@control.ee.ethz.ch.
1

2

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

from that of the buildings. A number of papers treat the building energy demands as exogenous signals,
which are typically estimated using building simulation environments such as EnergyPlus [11]. These studies
focus on the control of the devices within the energy hub, employing either deterministic [12,13] or stochastic
[14] formulations.
A large body of literature addresses the building control problems using deterministic schemes [6–9], as they
are scalable and suitable for problems with long prediction horizons. However, deterministic schemes suffer
from frequent constraints violations due to their inability to handle the system disturbances [15]. Stochastic
schemes can potentially address this issue but very often they are either not practically scalable (e.g., see
the discussion in [16] on the application of the scenario approach [17, 18] to building control problems), or
their reliability strongly depends on the exact knowledge of the distributional characteristics of the system
disturbances [19]. Poor estimation of these underlying distributions can lead to significant performance
deterioration [20].
To tackle this issue, distributionally robust methods have recently been introduced in the literature [21].
Instead of considering an exact distribution, these methods account for all distributions in a family which
shares only a few structural parameters, such as moments and/or support information [22,23]. In this context,
the historical data are only partially exploited to obtain a rough estimate of these structural parameters. On
the other hand, purely data-driven approaches which systematically exploit the historical data information to
determine appropriate families of distributions, have also been suggested [24, 25]. These methods provide the
tools to reformulate an originally infinite dimensional problem to a finite dimensional convex semi-definite
optimization problem (SDP). Typically, the number of constraints in this problem depends on the historical
data size, which limits the scalability of the method when dealing with a large data set. Although polynomial algorithms exist for solving SDPs [26], these algorithms are computationally demanding limiting the
applicability of these methods to small problem instances.
Our goal is to develop a data-driven stochastic control scheme that is capable of cooperatively operating the
energy hub and the district buildings. This paper extends the preliminary work in [27], providing a refined
robust approach on handling the system disturbances by systematically exploiting the available historical
realizations of these stochastic processes. In particular:
(1) We propose a data-driven approach that exploits the historical data to train linear models of the
exogenous disturbances, and construct families of distributions that encompass the true disturbance
distributions with high confidence.
(2) We exploit the structure of these families of distributions to formulate a robust multistage stochastic
optimization problem that minimizes the wost-case expected energy costs of the system. The size of
the resulting optimization problem is independent of the historical data size. We approximate this
infinite dimensional problem by a finite-dimensional linear program that scales polynomially with
respect to the prediction horizon length, and more importantly, unlike the SDP approaches cited
above, it can effectively be solved for large scale systems.
(3) We demonstrate in an extensive numerical study the efficacy of the proposed method which is shown
to outperform the optimally tuned deterministic equivalent in terms of energy consumption and
constraint violations. We observe that higher cost benefits are obtained by merging buildings with
dissimilar operation plans, rather than diverse construction characteristics.
The paper is organized as follows. In Section 2, we review in a more compact way, the modeling approach
presented in [27]. The main contributions of this paper are summarized in Sections 3 and 4 where the
developed data-driven distributionally robust methods, and the techniques associated with the derivation of

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

3

Energy Hub
Photovoltaics

Battery

Electricity, dt,elec

pin
t,el

Electrical pout
t,el
grid
Cooling, dt,cool
Heat pump
Chiller

pin
t,gas
Gas grid

Heating, dt,heat
Boiler

District

Figure 1. Heating, cooling and electricity network of a district.
a tractable approximation to the infinite dimensional stochastic optimization problem, are discussed. We
conclude this paper with an extensive numerical study performed in Section 5. The proofs of the propositions
can be found in the Appendix.
Notation: All random vectors appearing in this paper are defined on an abstract probability space
(Ω, F, P), where E(·) denotes the expectation operator with respect to P. Random vectors are represented in
boldface, while their realizations are denoted by the corresponding symbols in normal font. For given matrices
(A1 , . . . , Am ), we define A := diag(A1 , . . . , Am ) as the block-diagonal matrix with elements (A1 , . . . , Am )
> >
] ∈ Rk with
on its diagonal. Given vectors (v1 , . . . , vm ), vi ∈ Rki , we define [v1 , . . . , vm ] := [v1> , . . . , vm
Pm
k = i=1 ki , as their vector concatenation. We denote by 1 and 0 the vectors with components all one and
zero, respectively. The dimension of the corresponding vectors and the vector concatenations, will be clear
from the context.

2. System modelling
In this section, we describe the energy hub and building dynamics using discrete time, bilinear models
affected by stochastic exogenous disturbances. We assume that these disturbances evolve according to stochastic processes {ξ t }t∈T , where T = {1, . . . , T }, and T is the length of the horizon considered. The vector ξ t
encompasses all the stochastic processes affecting the energy hub and building dynamics, such as the ambient
temperature, solar radiation and internal gains of the buildings.

2.1. Energy hub dynamics
We visualize the energy hub as a conceptual entity that houses and interconnects a number of conversion,
storage and production devices that are shared by the building community. An energy hub essentially provides
the interface between the energy grid and the building community. As depicted in the illustrative example
of Fig. 1, the energy hub is capable of (i) purchasing electricity and gas to meet the electricity, cooling and
heating demand of the building community, and (ii) selling the electricity produced by the photovoltaics or
stored in the battery to maximize profit. In the following, we define the set K to include the energy hub devices,

4

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

and the sets G and O to include the energy streams from the grid and building community, respectively. In
the example of Fig. 1, K = {Battery, Chiller, Heat pump, Boiler, Photovoltaics}, G = {Electrical grid, Gas
grid} and O = {Electricity, Cooling, Heating}.
We model every device i ∈ K using linear dynamics and constraints, as follows:
xt+1,i = Ai xt,i + Bi ut,i + Ci ξ t ,

(1)

Fx,i xt,i + Fu,i ut,i + Fξ,i ξ t ≤ hi .
The vectors xt,i and ut,i denote the hub device internal states and control inputs, respectively. The matrices
Ai , Bi and Ci have appropriate dimensions and can be derived from the device characteristics. Finally, we
assume linear operational constraints captured by the matrices Fx,i , Fu,i , Fξ,i and hi . As discussed in [10,13],
linear approximations for the energy hub devices are reasonable at this level of abstraction. Examples of such
models for the numerical study of Section 5 are given in Appendix 6.
The energy balancing nodes are used to model the interconnection of the energy hub devices, as follows:
Hp pt + Hu ut + Hd dt = 0,

(2)

out
in
out
where the vector pt = [pin
t,1 , pt,1 , . . . , pt,|G| , pt,|G| ] contains the non negative decision variables of the grid
in
out
energy streams, with {pt,i }i∈G and {pt,i }i∈G denoting the energy purchased from, and sold to the grid,
respectively. In a similar way, the vector dt = [dt,1 , . . . , dt,|O| ] concatenates the demands of the building
community, and ut = [ut,1 , . . . , ut,|K| ] captures the decision variables of the energy hub devices. The matrices
Hp , Hu and Hd have proper dimensions and model the power flows affecting the respective balancing node.

2.2. Building dynamics
We model the building dynamics using bilinear state space models, motivated by the resistance-capacitance
models of [28]. The accuracy of these bilinear models was validated against established building simulation
software [11], and real buildings [9]. We denote by B the set of district buildings and for each building i ∈ B
we assume a bilinear model that captures the temperature evolution of the rooms, walls, ceiling and floors,
as follows,


>
xt+1,i = Ai xt,i + Bi + x>
(3)
t,i Ei ut,i + Di + v t,i Ci ξ t ,
where ut,i contains the inputs to these building actuation systems which are not coupled with the disturbances
(e.g., radiators, thermally activated building structures (TABS), air handling unit (AHU), floor heating). The
rest of the control inputs (e.g. , position of the blinds) are captured by v t,i . The system matrices Ai , Bi ,
Ci , Di and Ei depend on the specific building characteristics (e.g., number of rooms, construction material,
window fraction area and actuation units).
The state to input bilinear terms in (3) severely complicate the design of computationally tractable controllers. To address this issue, we replace the bilinear terms x>
b>
t,i Ei ut,i , with the linear terms x
i Ei ut,i,j , where
x
bi is the initial condition for the states of building i. Therefore, the dynamics of the approximated system
are given as:
xt+1,i = Ai xt,i + Bi (b
xi ) ut,i + Ci (v t,i )ξ t + Di ξ t ,

(4a)

The matrix Bi (b
xi ) and the matrix function Ci (v t,i ), are readily constructed by the primitive description of
the dynamics in (3). More sophisticated approaches could also be envisioned (e.g. linearizing about a state
trajectory), but we do not pursue this further here for the shake of simplicity.
We consider linear operational constraints (e.g., radiator and AHU limitations), as follows:
Fx,i xt,i + Fu,i ut,i + Fv,i v t,i + Fξ,i ξ t ≤ hi ,

(4b)

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

5

where the matrices Fx,i , Fu,i , Fv,i , Fξ,i and fi are derived using the BRCM Toolbox [28], and are of appropriate
dimensions. In addition to the operational constraints in (4b), we consider user specified comfort ranges given
as,
lbt,i ≤ xt,i ≤ ubt,i ,

(4c)

where lbt,i and ubt,i are lower and upper bounds. One typically considers bounds only on the room temperatures, so many of the upper/lower bounds can be assumed to be plus/minus infinity, effectively eliminating
the corresponding constraints. These bounds may vary during the day to reflect occupancy patterns (e.g.,
office buildings temperature bounds are often relaxed during the night time since the building is empty). We
will refer to constraint set (4c) as the comfort constraints [29].

2.3. Coupling of buildings to the energy hub
We model the coupling between the energy hub and the buildings with the following set of equality
constraints:
X
dt,j =
η>
∀j ∈ O,
(5)
i,j ut,i ,
i∈B

where η i,j is a (0, 1)-matrix that models whether the building actuation system ut,i is connected to the j-th
output energy source of the hub, dt,j .
To simplify notation, we compactly rewrite Eqns (1), (2), (4) and (5), as follows:
x = B(b
x)u + C(v)ξ + Dξ,
(6)
F π π + F ξ ξ ≤ h,
where the concatenated vector π is defined as π = [p, d, x, u, v]. Note that all the vectors in (6) are considered
over the time horizon T , e.g., p = [p1 , . . . , p|T | ], while the vectors x, u are defined such they include the
state and input variables, respectively, for both the energy hub devices and the buildings.

3. Modelling the uncertainty
To model the evolution of the stochastic process ξ t the exact characterization of its distribution is needed.
However, this information is typically unavailable. Nevertheless, historical data such as past realizations of
ξ t are usually available. A simple approach is to construct an empirical distribution using the historical data
set, and then use it as a proxy for the true distribution. However, if the data set is small, then there might
be several distributions that can describe these data points. As indicated in [20], if one arbitrarily chooses an
element from this family of distributions, e.g., the empirical distribution, then the solution of the resulting
stochastic optimization problem can differ significantly from the solution in which the true distribution is
used. In the following, we adopt a robust perspective to this problem by training linear models of the
exogenous disturbances and constructing families of distributions that describe the disturbance realizations
during past years. We then formulate an optimization problem whose solution addresses all distributions in
the constructed family.

3.1. Dynamics of the stochastic process
We denote by D the set that contains all the sources of disturbances appearing in our problem, e.g.,
solar radiation, ambient temperature, internal gains. For each i ∈ D, we model the disturbance ξ t,i as a

6

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

deterministic forecast ft,i , plus the stochastic error term et,i . To simplify notation, we omit index i in the
subsequent discussion, i.e.,
ξ t = ft + et .

(7a)

Weather forecasts are easily accessible from national weather services (e.g., COSMO-7 of MeteoSwiss [30])
while forecasts of the anticipated internal gains can be obtained from standard weekly profiles of typical
building configurations (e.g., office and residential buildings [29]). The dynamical evolution of et , is captured
by a first order autoregressive system,
et+1 = αt et + wt ,

(7b)

where αt is time-varying constant, and wt ∈ Wt is the stochastic process governing the stochastic evolution
of the error. The autoregressive evolution of et is motivated by systematic errors between forecast and actual
realization. Indeed, the forecast provided by MeteoSwiss can differ from the true realization of the uncertain
parameters due to, among other things, the spatial difference between the local weather station and the
building leading to error correlation over time. Additionally, wt models the noise from imperfect forecasting
and possible measurement noise from the sensing devices.
Equation (7) suggests that instead of constructing the probability distribution associated with ξ t , we can
equivalently construct the probability distribution associated with wt which we denote by Pt . We make the
following structural assumptions for the random variables wt : (i) they are normally distributed, (ii) they are
mutually independent for all time stages t ∈ T and disturbances i ∈ D, and (iii) they are stationary with
respect to different days (i.e., αt and the distribution of wt for the same t but for different days is the same).
Essentially these assumptions imply that first order models adequately capture the correlation over time,
and one can ignore the correlation between different days and disturbances. We validate these hypotheses in
Section 5.2 based on real disturbance data.

3.2. Data-driven uncertainty sets
For each disturbance i ∈ D, the historical data set consists of realization and forecast pairs of the form
{(ξtk , ftk )} with t denoting the time of the k-th day that the data was recorded. We consider N records of
forecast and realization compatible with the sampling time used for the models of Section 2 (typically, daily
records for several years sampled hourly).
We calculate the constants αt describing the dynamics in (7) using least squares fitting. In particular, for
each t ∈ {1, . . . , 24}, we solve,
min

αt ∈R

N
X

2
k
k
(ξt+1
− ft+1
) − αt (ξtk − ftk ) ,

(8)

k=1

and using the optimal solution construct the residual data points as follows:
k
k
wtk := (ξt+1
− ft+1
) − αt (ξtk − ftk ), k = 1, . . . , N.

We denote by St := {wtk }N
k=1 the set of residuals derived from the historical data, and subsequently denote
2
by µ
bt and σ
bt their empirical mean and variance, respectively.
We construct a family of distributions, Pt , that are compatible with the residuals of our historical data
St . The set Pt is sometimes referred as the ambiguity set in the robust optimization literature [21, 24]. We
consider an ambiguity set, Pt , of the following form:
Pt = {N (µt , σt2 ) : µt ≤ µt ≤ µt , σ 2t ≤ σt2 ≤ σ 2t }.

(9a)

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

7

The constants (µt , µt , σ 2t , σ 2t ), are selected to ensure that given St , the true distribution, Pt of wt , is an
element of Pt with high probability. This is achieved by utilizing concepts from statistical hypothesis theory.
In particular, the statistic associated with the chi-square hypothesis test, hχ := (N − 1)b
σt2 /σt2 , follows the
chi-square distribution with (N − 1) degrees of freedom [31]. To ensure that σ 2t ≤ σt2 ≤ σ 2t , with probability
at least 1 − δtχ , we set,
χ
χ
σ 2t = (N − 1) σ
bt2 /q(N
−1) (δt /2),
χ
χ
σ 2t = (N − 1) σ
bt2 /q(N
−1) (1 − δt /2).

(9b)

χ
Here, q(N
−1) (·) denotes the quantile function of the chi-square distribution with (N − 1) degrees of freedom.

We construct the bounds of µt in a similar way. Given that the wt follows a normal distribution and its
p
variance σt2 follows the chi-square distribution, hst := (b
µt − µt )/ σ
bt2 /N , follows a student t-distribution with
N − 1 degrees of freedom, [31]. To ensure that µt ≤ µt ≤ µt , with probability at least 1 − δtst , we set,
p
st st
(δt /2) σ
µt = µ
bt + qN
bt2 /N ,
p
st st
(δt /2) σ
bt − qN
bt2 /N ,
µt = µ

(9c)

st
(·) denotes the quantile function of the Student’s t-distribution with (N − 1) degrees of freedom.
where, qN

The bounds given in (9) explicitly determine the family of distributions, Pt . The following proposition
provides the confidence, 1 − δt , by which the true disturbance distribution, Pt , is a member of Pt .
Proposition 1. Given the sample data, St , let PSt be the |N |-fold product distribution of Pt , then

PSt Pt ∈ Pt ≥ 1 − δt ,
where δt = δtst + δtχ .
Let us now define the multivariate stochastic process w = [w0,1 , . . . , w|T |,|D| ] ∈ W , with W = W0,1 × . . . ×
W|T |,|D| . We characterize the joint distribution P of w, as follows,
n
P = N (µ, Σ) : with µ = (µ1,1 , . . . , µ|T |,|D| )
2
2
and Σ = diag(σ1,1
, . . . , σ|T
|,|D| )o
s.t.P1,1 ∈ P1,1 , . . . , P|T |,|D| ∈ P|T |,|D| .

We close this section, by compactly rewriting the disturbance modelling equations over the horizon T , as
follows:
ξ = H(b
e)w,
where the matrix H(b
e) is readily constructed from (7) and depends on the forecast, and the vector eb =
[e0,1 , . . . , e0,|D| ] containing the current error disturbance measurements.

4. Problem formulation
Our objective is to minimize the worst-case expected cost of energy purchased from the grid by the building
community, over a finite horizon. This must be achieved while satisfying the dynamics and constraints of the
devices and the buildings in the system. This problem can be formulated as a multistage stochastic linear

8

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

program, as follows:
min sup EP c> p



P∈P

s.t. p, u ∈ C, v ∈ R, π = [p, d, x, u, v],
ξ = H(b
e)w,

(10)

x = B(b
x)u + C(v)ξ + Dξ,


inf P F π π + F ξ ξ ≤ h ≥ 1 − ,
P∈P

out
in
out
where ct = [cin
t,1 , ct,1 , . . . , ct,|G| , ct,|G| ] contains the time varying prices
the prices {cin
t,i }i∈G for purchasing energy from the grid are positive

of the grid energy streams. Note that
scalars, while {cout
t,i }i∈G are negative.
Furthermore, we select a chance constraints formulation with worst-case violation probability . The decision
variables v that are coupled with the disturbances ξ, are allowed to take values in the generic finite-dimensional
vector space R. We choose the decision variables p, u to be strictly causal disturbance feedback policies. For
instance, pt is the strictly causal, vector valued function of the energy purchased from the grid at time t ∈ T ,
defined as pt : W1 × . . . × Wt−1 → R. In the following, we denote by C, the infinite-dimensional function
space of strictly causal disturbance feedback policies. Finally, the expected cost of the objective has been
formulated to account for the worst-case multivariate distribution P in the ambiguity set P.

4.1. Constraint relaxation
The proposed control architecture is based on a receding horizon implementation of the system. In this
setting, the mismatch between the building prediction model and real system, can lead to comfort bounds
violations. To address infeasible instances of Problem (10) due to this issue, we relax the comfort constraints
(4c), as follows:
max{lbt,i − xt,i , 0, xt,i − ubt,i } ≤ st,i ,
where we refer st,i as the slack variable. To this end, we rewrite the compactly formulated inequality (6), as
follows:
F π π + F ξ ξ − s ≤ h,

(11)

where s ∈ C denotes the concatenated slack variable vectors over the horizon T . Notice that the slack
variables are only introduced to deal with building model mismatch which can be an issue on a receding
horizon implementation. On the other hand, the chance constraint formulation in Problem (10), is primarily
used to address in a probabilistic fashion the extreme realization of the uncertain parameters w.

4.2. Optimization over linear feedback policies
The optimal solution of Problem (10) remains intractable due to the infinite dimensional structure of its
decision variables. However, a tractable approximation of Problem (10) can be obtained by restricting the
decision variables to the finite dimensional space of affine policies denoted as Caff . A strictly causal affine
policy, e.g. for the variables pt is given by
pt = p0,t +

t−1
X

Ps,t ws ,

(12)

s=1

where p0,t ∈ R, and matrices Ps,t of appropriate dimensions. We refer to (12) as the strictly causal affine
decision rule (ADR) [32]. We refer to the linear policy in (12) which disregards the disturbance history (e.g.,
pt = p0,t ), as open loop policy (OLP). The non-adaptive nature of the OLP policy provides an even more

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

9

conservative controller parametrization with the benefit of a considerably smaller number of optimization
variables.
In this context, the approximated variant of Problem (10) is given as follows:

min sup EP c> p + γ1> s
P∈P

s.t. p, u, s ∈ Caff , v ∈ R, π = [p, d, x, u, v],
ξ = H(b
e)w,

(13)

x = B(b
x)u + C(v)ξ + Dξ,


inf P F π π + F ξ ξ − s ≤ h ≥ 1 − ,

P∈P

where the additional term γ1> s in the objective penalizes the constraint violations in (11), with parameter
γ ∈ R+ . A discussion for appropriate values of γ can be found in [33].

4.3. Chance constraint approximation
The main body of the literature exploits the structure of the ambiguity sets to propose semi-definite [22,25]
and second-order cone [34] reformulations of the distributionally robust chance constraints. Although convex
problems with conic constraints are generally tractable, they are computationally demanding for large systems.
c ⊆ W such that the feasible region of the corresponding robust
An alternative approach is to construct a set W
constraint is a subset of the feasible region of the distributionally robust chance constraint, i.e.,
c,
if F π π + F ξ ξ − s ≤ h, ∀ w ∈ W


then inf P F π π + F ξ ξ − s ≤ h ≥ 1 − .

(14a)
(14b)

P∈P

c described by non-linear constraints. NevSuch methods are discussed in [24] where the authors compute W
ertheless, these approaches typically lead to semi-definite reformulations of the robust constraint (14a). To
tackle this issue, we resort to a more stringent condition, in which we require that at least, 1 − , of the
c . To this end, we construct
probability mass of each distribution, P, in the ambiguity set P, is contained in W
ct,i , as follows:
for every i ∈ D and t ∈ T , the compact convex set W


n
ct,i = wt,i | wt,i ≥ µ − Φ−1 1 − β
σ t,i ,
W
t,i
t,i


o
(15a)
wt,i ≤ µt,i + Φ−1 1 − β t,i σ t,i .
where the constants µt,i , µt,i , and σ t,i are given in (9c) and (9b), respectively. Moreover, Φ−1 (·) denotes the
inverse cumulative normal distribution function, and β t,i and β t,i are positive constants chosen as,
XX


β t,i + β t,i = .

(15b)

t∈T i∈D

c to be,
Finally, we set W
c=W
c1,1 × . . . × W
c|T |,|D| .
W
c defined in (15). Then, the following probabilistic guarantee holds:
Proposition 2. Let W
c ) ≥ 1 − .
inf P(w ∈ W

P∈P

(15c)

10

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

Notice that satisfying (14a) together with the condition in Proposition 2 is more stringent than the actual
chance constraint. In other words,
c ) ≥ 1 − , and (14a) holds,
if inf P(w ∈ W
P∈P


then inf P F π π + F ξ ξ − s ≤ h ≥ 1 − .
P∈P

Imposing this additional condition, we further restrict the feasible region with the benefit of gaining compuc . This, in turn, allows us to address
tational tractability due to the simple hyperrectangular structure of W
the large scale problems examined in this paper.

4.4. Objective function reformulation
We employ the epigraph representation to equivalently rewrite the worst-case expectation in the objective
function of Problem (13). In particular, the linearity of the objective function allows us to replace the expected
value with:


e + γ1> e
sup EP c> p + γ1> s = sup c> p
s = τ,
(16)
P∈P

P∈P

which can equivalently be written as,
e + γ1> e
c> p
s ≤ τ, ∀µ ∈ [µ, µ],

(17)

e, and e
where p
s are derived from (12) by replacing the stochastic variable w by its expected value µ (e.g.,
P
et = p0,t + t−1
p
s=1 Ps,t µt ). Notice that the epigraph representation in (17) is exact with the scalar variable τ
replacing the original worst-case objective.
We conclude this section by providing the final approximation of Problem (10) as follows:
min τ
s.t. p, u, s ∈ Caff , v ∈ R, π = [p, d, x, u, v],
e + γ1> e
c> p
s ≤ τ, ∀µ ∈ [µ, µ],

ξ = H(b
e)w,



c.
∀w ∈ W
x = B(b
x)u + C(v)ξ + Dξ,



F π π + F ξ ξ − s ≤ h,

(18)

Problem (18) still retains its infinite structure involving a continuum space of decision variables and constraints. However, as indicated in [32, 35], Problem (18) can be reformulated into a linear optimization
problem by employing traditional robust optimization techniques that use duality to translate the semiinfinite structure into a finite number of linear constraints. Unlike other robust convex programs which
employ semi-definite programming, this problem is computationally tractable and can be solved in very high
dimensions.

5. Numerical results
In this section, we perform numerical studies to assess the performance of the ADR and OLP control
methodologies. We compare these methods to the classical and commonly used in practice certainty equivalence problem (CEP) in which all random variables of Problem (10) are replaced with their expected value.
To make the comparison of these control strategies more transparent, we do not consider the possibility of
selling energy to the grid. Moreover, we analyze the sensitivity of the considered control methods with respect
to parameters deviations, and investigate the potential gains that can be obtained by cooperatively managing
buildings with heterogeneous construction and operation characteristics.

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

11

Energy Hub
Battery

PV

Electricity
Electricity
grid
Cooling
Heat pump
Chiller

Heating
Boiler

District

Figure 2. Energy hub configuration for numerical studies.

5.1. Problem configuration
We consider districts composed of buildings with roughly the same floor area which are connected through
an energy hub that comprises five devices: chiller, boiler, heat pump, photovoltaics (PV) and battery, as
depicted in Fig. 2. To keep the ratio between demand and supply in the energy hub, relatively constant, we
linearly scale the capacities of the hub devices based on the number of buildings in our system. A detailed
description of the hourly discretized dynamics and constraints governing the energy hub components can be
found in the Appendix 6.
The buildings have heterogeneous construction characteristics, summarized in Table 1. In particular, we
use hourly discretized building models described in [9] characterized by the building type BT ∈ {Swiss Passive
(SP), Swiss Average(SA), Swiss Target (ST)}, and construction type CT = {heavy, light}. Each building
consists of 5 rooms which are characterized by the window fraction area WFA = {30%,50%,80%} and their
corresponding facade orientation. The dynamics and constraints of each building are generated using the
BRCM toolbox [28] to which the individual building construction details along with the specifications of the
control devices (radiators, AHU, TABS, blinds) are provided.
Each of the building models in Table 1, consists of 113 states. Model reduction techniques are used to derive
a simplified, but sufficiently accurate, building model. For instance, using the bilinear balanced truncation

2

No.

Area(m )

1
2
3
4
5

420
420
441
441
374

Building specifications
WFA BT
CT
Input Devices
30%
50%
80%
50%
50%

SP
SP
SA
SA
ST

heavy AHU, blinds, radiator
light AHU, blinds, TABS
light AHU, blinds, TABS
heavy AHU, blinds, radiator
heavy AHU, blinds, radiator

Table 1. Summary of the 5 buildings used in the simulations.

12

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

Time

Cost

Winter Bounds

Summer Bounds

◦

05:00 - 23:00 0.145 CHF/kWh
23:00 - 05:00 0.097 CHF/kWh

[20, 23]◦ C
[15, 30]◦ C

[21, 25] C
[15, 30]◦ C

Table 2. Electricity day/night tariff variations and comfort constraints bounds.

method described in [36], we generate a 59 state model, with a maximum absolute error of less than 0.1◦ C,
during step response simulations for building 1 of Table 1. By evaluating the Hankel matrix of this building
system, we identified the time constants associated with the three most controllable and observable modes as
11.74 days, 4.47 hours and 4.24 minutes. Approximately the same time constants were identified for the other
building systems in Table 1. The reduced order linear models are used only for prediction purposes, while
the original bilinear models, given in (3), are employed to simulate the buildings dynamics in the closed-loop
implementation.
We compare the performance of the ADR, OLP and CEP control designs using the metrics of purchased
grid energy, room constraint violations and solution time. The cost of the energy purchased from the grid
is measured in Swiss Franc (CHF), and the room constraint violations are measured in Kelvin hours (Kh).
We assume time-varying electricity tariffs and comfort constraints bounds, as given in Table 2. The same
comfort bounds are used in every room and building considered in the system. The disturbance forecasts
and realizations are the same as those used by the OptiControl project [37], for the city of Zürich during
the years 2006 and 2007. The data set from year 2006 is used to train the linear disturbance models and
construct the families of distributions given in equations (8) and (9). The 2007 dataset is then used to test
the performance of the controllers developed based on the resulting distribution families. Finally, we choose
χ
st
as soft constraint penalization, γ = 103 , confidence levels δt,i
= δt,i
= 0.01, and constraint violation level
 = 0.01. We select the violation levels for the upper and lower bounds of the i-th disturbance at time t, as
β t,i = β t,i = /(2|T ||D|).

5.2. Disturbance model verification
We consider a disturbance set D which comprises seven sources of uncertainty; ambient and ground temperatures, four sources of solar radiation (North, South, West and East), and building internal gains. The
analysis in Section 3 requires that the residual uncertain parameters wt are normally distributed and independent over time. We verify these assumptions using the Shapiro-Wilk [38] and Pearson [31] hypothesis
tests on historical, weather and occupancy, data for the year 2006.
Time

AT

SRS

SRE

IG

07:00
08:00
09:00
10:00
11:00
12:00
13:00
14:00

0.989
0.996
0.994
0.993
0.994
0.992
0.997
0.995

0.991
0.997
0.995
0.965
0.990
0.941
0.996
0.995

0.972
0.996
0.993
0.992
0.987
0.981
0.988
0.984

0.940
0.945
0.923
0.967
0.920
0.928
0.945
0.965

Table 3. Statistic of the Shapiro-Wilk normality test.

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

13

In Table 3, we report the values of the Shapiro-Wilk test statistic from 7:00 to 14:00 for the ambient
temperature (AT), solar radiation south (SRS) and east (SRE), and internal gains (IG). We choose to present
the results from 7:00 to 14:00 since these are the hours with the greatest variation in mismatch between
forecasts and realizations. Values which are close to 1 indicate that the sample data are compatible with a
normal distribution. In particular, the normality hypothesis is accepted at the significance level of 0.1 for
the atmospheric processes, and 0.01 for the internal gains. These results are in accordance with the study
in [39] where it is shown that normal distributions sufficiently capture the evolution of weather processes.
By contrast, in [40] it is argued that Poisson distributions should be used to generate occupancy profile
trajectories. Although, Poisson distributions can be approximated by Gaussians [41], we stress that the
bounds in (9) are useful even if the underlying data are mildly non-Gaussian, as suggested in [31, §11.3].
Time 07:00 08:00 09:00 10:00 11:00 12:00
07:00
1.00 0.04 0.02 0.01 -0.01 0.01
08:00
0.04 1.00 0.05 0.05 -0.01 0.00
09:00
0.02 0.05 1.00 0.06 0.03 -0.09
10:00
0.01 0.05 0.06 1.00 0.00 -0.09
11:00 -0.01 -0.01 0.03 0.00 1.00 -0.03
12:00
0.01 0.00 -0.09 -0.09 -0.03 1.00
Table 4. Statistic of the Pearson correlation test.

In Table 4, we provide the values of the Pearson test statistic from 7:00 to 12:00 for the SRS and the
results are similar for the other sources of uncertainty. This correlation statistic can range between plus
and minus one with values close to zero denoting probably uncorrelated data. We emphasize that for every
disturbance, the hypothesis of uncorrelated data is accepted at the significance level of 0.05 which provides a
strong statistical evidence for the validity of our uncertain disturbance assumptions.

5.3. Prediction horizon selection
We investigate the effect of the prediction horizon on the performance of the ADR, OLP and CEP control
methodologies. In Fig. 3, we show the respective curves generated by conducting a receding horizon simulation
during the first week of January 2007 for building 1 in Table 1. We observe that the cost of purchased grid
energy associated with the CEP solution method is the least, at the expense of frequent comfort constraint
violations. On the contrary, the stochastic approximations (OLP and ADR), are less cost efficient but lead
to many fewer constraint violations.
We can identify several factors associated with the selection of a suitable prediction horizon:
(i) The dominant time constants of the buildings were roughly identified as 11 days, 4 hours and 4
minutes. Hence, a horizon of T ≥ 7 hours facilitates the controller to excite the modes of the system
that determine its short term evolution. We emphasize that contrary to the ADR and CEP methods,
the long horizons deteriorate the performance of the OLP method due to its inability to adapt on the
growing size of the disturbance uncertainty.
(ii) The six hour gap between day and night comfort bounds, given in Table 2. A prediction horizon of
T ≥ 6, is required for the system to anticipate the comfort bounds of the next day. In this way, it
can utilize the heating capacity of the buildings and the battery storage to exploit the day-night tariff
structure of the electricity prices.

14

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

Energy cost (CHF/Week)

Grid purchased energy
100
95
90
CEP
OLP
ADR

85
80
75
70
65
2

4

6

8

10

12

14

16

18

20

22

24

horizon (hours)
Room constraint violations

Violations (Kh/Week)

12

CEP
OLP
ADR

10
8
6
4
2
0
2

4

6

8

10

12

14

16

18

20

22

24

16

18

20

22

24

horizon (hours)

Time (sec)

Solution time

10

2

10

1

CEP
OLP
ADR

2

4

6

8

10

12

14

Horizon (hours)

Figure 3. Purchased grid energy, room constraint violations and solution time with respect
to the prediction horizon length T during a typical winter week.
(iii) The renewable energy production peaks around midday. A horizon of T ≥ 6, is sufficient for the
system to anticipate the cost-free energy that will be available from the photovoltaic units. In this
way, it can efficiently utilize the battery to fully exploit the potential excess of energy.
(iv) The computation time associated with the ADR and OLP control methods which increases with the
prediction horizon length.
The above analysis suggests that a prediction horizon of T = 8 hours provides a reasonable trade-off between
foresight of the controller and computational tractability. This value is used for the rest of this section.

5.4. Solution method selection
To generalize our observations regarding the comparison of the ADR, OLP and CEP control methodologies.
The simulation experiment is extended to a district composed of the 5 buildings summarized in Table 1. The
system is simulated in a receding-horizon fashion using data realizations of 12 consecutive weeks (restarting
at the beginning of each week) for the winter and summer periods of 2007, starting January 1st and June
29th, respectively.
We calculated the purchased energy cost and room constraint violations for each one of these 12 weeks,
and we report the results in Table 5. The table entries correspond to the (empirical mean, empirical standard

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

Method

Cost (p.u.)

Winter
Violations (Kh/Week)

CEP
OLP
ADR

(1.00, 0.21)
(1.34, 0.19)
(1.07, 0.21)

(4.31, 0.68)
(0.00, 0.00)
(0.03, 0.01)

Method

Cost (p.u.)

Summer
Violations (Kh/Week)

CEP
OLP
ADR

(1.00, 0.17)
(1.07, 0.18)
(1.03, 0.16)

(2.23, 0.14)
(0.03, 0.01)
(0.06, 0.01)

15

Basis (CHF/Week)
327.56

Basis (CHF/Week)
64.19

Table 5. Receding horizon performance results.

deviation) over these 12 weeks scenarios. Notice that we present the energy costs in the per unit (p.u.) system
using as base values the mean costs occurring for the CEP method during the winter and summer period,
respectively. Once again, we observe that the CEP solution method is the most cost effective at the expense
of significant comfort constraint violations. The ADR and OLP control methodologies result in higher cost
of purchased grid energy but results in fewer constraint violations. Notice that the ADR approximation is
considerably more cost efficient than the OLP, while it achieves the same level of constraint violation. This
is attributed to the nature of the ADR method which takes into account the potential adaptation of future
decisions based on the realizations of the disturbance variables.
To better visualize the behavior of the three approximations, the trajectories for the mean room temperature of Building 1 (c.f. Table 1), and the amount of purchased grid energy, are depicted in Figure 4 for
Monday, 1st January 2007. The CEP method operates very near comfort constraints and leads to frequent
violations, while the OLP produces more conservative results by keeping the room temperature well inside
the comfort range at the expense of consumption cost benefits. The poor performance of the CEP method
with respect to the constraint violations can be explained when considering that the approximation assumes
a deterministic evolution of the uncertain parameters. On the other hand, the adaptive nature of the ADR
produces a good compromise between the optimistic decisions made by the CEP and the conservative decisions of the OLP. We remark that all three approximations utilize the load shifting capabilities of the battery
by storing energy during the evening hours when electricity is cheaper, and deploying that energy in the early
morning hours when the building needs to be brought back within the comfort range. As expected the energy
produced from the photovoltaic unit is fully exploited to either reduce the grid energy consumption or charge
the battery.

5.5. Comparison of ADR and tuned-CEP methods
We compare the stochastic ADR method with the deterministic tuned-CEP one. We denote by tuned-CEP,
the CEP method in which the comfort bounds of the building, i ∈ B, at time t ∈ T , have been tightened by
a fixed constant cb ≥ 0 Celcius degrees, as follows,
min{xt,i − lbt,i , ubt,i − xt,i } ≥ cb .
This constant cb is computed as the minimum constraint tightening value for which the ADR and CEP
methods achieve the same amount of comfort bounds violations over the simulation horizon. The tuned-CEP

16

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

Energy cost (CHF/Week)

Grid purchased energy
100
95
90
CEP
OLP
ADR

85
80
75
70
65
2

4

6

8

10

12

14

16

18

20

22

24

horizon (hours)
Room constraint violations

Violations (Kh/Week)

12

CEP
OLP
ADR

10
8
6
4
2
0
2

4

6

8

10

12

14

16

18

20

22

24

16

18

20

22

24

horizon (hours)

Time (sec)

Solution time
CEP
OLP
ADR

10 2

10

1

2

4

6

8

10

12

14

Horizon (hours)

Figure 4. Profile trajectories for the mean room temperature of Building 1, and the total
purchased grid energy by the community during a typical winter day.
can be seen as a practical ad-hoc solution to deal with the increased constraint violations of the CEP method.

Method

Cost (p.u.)

Winter
Tightening cb (◦ C)

CEP
ADR

(1.00, 0.21)
(0.94, 0.19)

[0.38, 0.56]
0

Basis (CHF/Week)
327.56

Table 6. Comparison of ADR and tuned-CEP methods.

We conducted receding horizon simulations for 12 consecutive weeks of the winter 2007, and we report in
Table 6 the empirical mean and standard deviation of the recorded weekly costs. For each of these weeks, we
computed by repetitive trials the respective constraint tightening constant, cb , which ranges from 0.38◦ C and
0.56◦ C. The fact that cb varies considerably from week to week, also suggests that tuning the CEP method is
not a straightforward procedure. Additionally, we observe that the cost of purchased energy associated with
the tuned-CEP solution method is considerably higher than the ADR one. Therefore, the inherent ability of
the ADR method to deal with the uncertainty allows it to achieve an efficient trade-off between constraint
violations and energy cost without the extra effort of tuning.

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

17

5.6. District energy benefit
Finally, we investigate the potential energy cost savings that can be obtained by considering diverse,
heterogeneous buildings in a district. We define the heterogeneity on buildings by means of operational and
construction characteristics. Specifically, we distinguish two operation types for the buildings, commercial
(COM) and residential (RES). The day-night room temperature bounds defining these types are summarized
in Table 7. We select building (Bd) 1 and 3 of Table 1 to be classified as commercial, while building 2 and 4
as residential. Notice that these buildings also differ in their construction characteristics.

Time
06:00 - 09:00
09:00 - 19:00
19:00 - 23:00
23:00 - 06:00

COM Bounds
Lower Upper
15◦ C 30◦ C
21◦ C 25◦ C
15◦ C 30◦ C
15◦ C 30◦ C

RES Bounds
Lower Upper
21◦ C 25◦ C
15◦ C 30◦ C
21◦ C 25◦ C
15◦ C 30◦ C

Table 7. Comfort constraints bounds for commercial and residential buildings.

We conduct five receding horizon simulations during the first week of January 2007. In each of these
experiments, we consider a different combination of residential and commercial buildings from Table 1, and
we report the cost of purchased grid energy for the CEP, ADR and OLP methods in Table 8. Notice that in
each simulation experiment, we appropriately scale the hub devices with respect to the number of buildings
that are connected to the energy hub.

RES (Bd 2)
RES (Bd 4)
COM (Bd 1)
COM (Bd 3)
RES+RES (Bd 2+4)
COM+COM (Bd 1+3)
COM+RES (Bd 1+2)

ADR
1.13
1.15
1.28
1.31
2.27
2.58
1.51

OLP
1.39
1.41
1.47
1.49
2.80
2.96
1.90

CEP Basis (CHF/Week)
1.00
56.37
1.02
1.23
1.26
2.02
2.49
1.37

Table 8. Cost of purchased energy for various district configurations during a typical winter week.

We observe that the cost benefits obtained by merging buildings with different construction characteristics
but similar operation plans can be limited. On the contrary, significant gains are obtained by aggregating
buildings with dissimilar operation plans. There are two reasons for this: (i) The energy shifting mainly occurs
among buildings with complementary demand profiles. In this occasion, the efficient, but of limited capacity,
devices of the hub (e.g., heat pump) are fully utilized during the course of day; (ii) The energy produced by
the photovoltaics units is better exploited through the storage when the peak demands of buildings do not
coincide over time. This flexibility on the demand profiles results to a considerably less purchased energy
from the grid.

18

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

6. Conclusion
This paper presents a unified data-driven control framework for the problem of cooperatively managing
the aggregated energy demands of buildings in a district. It exploits the available historical data to train
linear models with additive uncertainties that effectively capture the evolution of the stochastic processes in
the system. The underlying distributions of these additive uncertainties are shown to belong to Gaussian
families of distributions which are constructed off-line using the historical data information. We exploit
the simple structure of these sets to approximate the resulting robust stochastic optimization problem by a
finite dimensional linear program. This program is tractable and its complexity scales polynomially. This is
particularly important when addressing large scale problems such as cooperative building energy management.
An extensive simulation study based on realistic data demonstrates the efficacy of the proposed method.
Among other things, our study demonstrates that additional cost benefits can be obtained by aggregating
buildings with heterogeneity in the demand profiles.
As future work, we note that our problem has a decoupled structure. Indeed, the linear structure of
the objective and and the weakly coupled structure of the constraints can be exploited in a number of
distributed and decentralized optimization algorithms. Additionally, the problem’s decoupled structure can be
exploited by recently developed optimization algorithms such as the alternating direction method of multipliers
(ADMM), for a fast numerical solution of the linear optimization problem.

Acknowledgments
The authors would like to thank the Building Science and Technology Laboratory (EMPA) for providing the
building characteristics and the occupancy data, and MeteoSwiss for making available the weather forecasts
and realizations. The authors would also like to thank Annika Eichler, Marc Hohmann and Ben Flamm for
fruitful discussions on the topic.

Appendix
Proof of Proposition 1
We construct the family of distributions, Pt as the intersection of (i) the family, Ptχ , of Gaussian distributions with unknown mean and bounds on variance, and (ii) the family, Ptst , of Gaussian distributions with
unknown variance and bounds on the mean, as follows:
Pt = {Pt such that (Pt ∈ Ptχ ) ∧ (Pt ∈ Ptst )},
where,
Ptχ = {Pt is Gaussian and condition (9b) is satisfied},
Ptst = {Pt is Gaussian and condition (9c) is satisfied}.
Given the data, St , equations (9b) and (9c) also provide the confidence levels associated with the respective
ambiguity sets,
PSt (Pt ∈ Ptχ ) ≥ 1 − δtχ ,
PSt (Pt ∈ Ptst | Pt ∈ Ptχ ) ≥ 1 − δtst .
Notice that the bounds in (9c) are derived under the assumption that the sample variance is following the
chi-square distribution. Therefore, the confidence level associated with the construction of the distribution
family, Ptst , is given by the conditional probability.

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

19

We use Bayes rule to evaluate the significance level, δt , associated with the ambiguity set Pt , as follows:
PSt (Pt ∈ Pt ) = PSt (Pt ∈ Ptχ ∧ Pt ∈ Ptst )
= PSt (Pt ∈ Ptst | Pt ∈ Ptχ ) · PSt (Pt ∈ Ptχ )
≥ (1 − δtst )(1 − δtχ )
= 1 − δtst − δtχ + δtst δtχ
≥ 1 − δtst − δtχ
= 1 − δt .
with the significance level, δt , given as, δt = δtst + δtχ .



Proof of Proposition 2
c
The uncertain vector w comprises |T ||D| uncorrelated components. Therefore, we seek a bounded set W
that is composed by 2|T ||D| inequalities that upper and lower bound every component i = 1, . . . , |T ||D|, of
c can be written as,
w. Denoting by wt,i and wt,i the upper and lower bound, respectively, the set W

c = w ∈ R|T | |D| : w ≤ wt,i ≤ wt,i , ∀t ∈ T , i ∈ D .
W
t,i
We require that,
c) ≥ 1 −  .
inf P(w ∈ W

P∈P

We exploit the conservative Bonferonni approximation [20] to decouple the joint chance constraints into a set
of individual ones, as follows:
XX

c) =
P(w ∈
/W
P(wt,i ≥ wt,i ) + P(wt,i ≤ wt,i )
i∈D t∈T 

XX
β t,i + β t,i = ,
=
i∈D t∈T

with β t,i , β t,i ≥ 0. Notice that in this case the Bonferonni approximation is exact since we are considering
mutually exclusive events (i.e. P(wt,i ≥ wt,i ∧ wt,i ≤ wt,i ) = 0).
2
Assume that wt,i ∼ N (µt,i , σt,i
), then the chance constraints,

P(wt,i ≤ wt,i ) ≥ 1 − β t,i ,
P(wt,i ≥ wt,i ) ≥ 1 − β t,i ,
are equivalently reformulated (see [23]), as follows,

wi ≥ µi − Φ−1 1 − β t,i · σi ,

wi ≤ µi + Φ−1 1 − β t,i · σi .
In this setting, the robustified individual chance constraints,
inf P(wt,i ≤ wt,i ) ≥ 1 − β t,i ,

P∈P

inf P(wt,i ≥ wt,i ) ≥ 1 − β t,i ,

P∈P

are equivalently be reformulated as follows,
)

wt,i ≥ µt,i − Φ−1 1 − β t,i · σt,i
∀µt,i ≤ µt,i ≤ µt,i ,

wt,i ≤ µt,i + Φ−1 1 − β t,i · σt,i
∀σ t,i ≤ σt,i ≤ σ t,i .

20

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

Exploiting the linearity of these constraints, we can explicitly compute the upper and lower bounds, as follows,

wt,i = µt,i − Φ−1 1 − β t,i · σ t,i ,

wt,i = µt,i + Φ−1 1 − β t,i · σ t,i .
c in (15)
This analysis provides the constructive proof for the structure of W



Dynamics and constraints of the energy hub devices
We choose to present the devices characteristics sized per building (Bd). In particular, we model the chiller,
boiler and heat pump using a coefficient of performance [12], which gives rise to the following constraints:
in
uout
t,chiller = 0.7ut,chiller ,
in
uout
t,boiler = 0.9ut,boiler ,
out
in
ut,HP = 3ut,HP ,
out
in
out
where ut,i = [uin
t,i , ut,i ] with ut,i , ut,i , denoting the input and output power flows on the i-th device,
respectively. The capacities of the conversion units are given as follows,

0 ≤ uout
t,chiller ≤ 20 kW/Bd,
out
0 ≤ ut,boiler ≤ 25 kW/Bd
0 ≤ uout
t,HP ≤ 5 kW/Bd,
We consider a South oriented photovoltaic array with maximum output of 4.10 kW/Bd. The photovoltaic
dynamics are generated by linearizing the non-linear model of [42] with respect to the ambient temperature
and the solar radiation,
0 ≤ uout
t,PV ≤ 0.1280 − 0.0019ξ t,AT + 3.7ξ t,SRS ,
for all t ∈ T . Note that the units of ξ t,SRS are measured in kW/m2 and can typically take values ξ t,SRS ∈ [0, 1].
We consider a lead-acid battery [43], with a 5kW/Bd capacity, giving rise to the following linear dynamical
system:
!
!
!
0.51 0.22
0.61
−0.83
xt+1 =
xt +
uin
uout
t +
t ,
0.47 0.78
0.25
−0.39
where the states and control are constrained by:
out
0 ≤ uin
≤ 8,
t , ut
1 ≤ xt,1 + xt,2 ≤ 5, xt ≥ 0,
0.62xt,1 + 0.27xt,2 − uout
≥ 0,
t
in
0.84xt,1 + 0.37xt,2 + ut ≤ 2.58,
0.73xt,1 + 0.73xt,2 + uin
t ≤ 3.66.

To this end, we provide the equations describing the interconnection of the energy hub with the building
community. In particular, the electricity balancing constraint is given by,
X
X
pt,elect +
uout
uin
t,i =
t,i + dt,elect ,
i∈E+

i∈E−

where E+ = {Photovoltaics, Battery} and E− = {Heat Pump, Chiller, Boiler, Battery}. Similarly, the
heating and cooling energy balancing constraints are given by,
X
out
uout
t,i = dt,heat , and ut,chiller = dt,cool ,
i∈H+

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

21

where H+ = {Heat Pump, Boiler}. Finally, the demand for electricity, heating and cooling of the building
community is given by
X
dt,elect =
ut,i,AHU ,
i∈B
X
dt,cool =
ut,i,TABS ,
i∈B
X

dt,heat =
ut,i,radiator + ut,i,TABS ,
i∈B

where TABS are typically used both for heating and cooling of the buildings.

References
[1] L. Perez-Lombard, J. Ortiz, and C. Pout. A review on buildings energy consumption information. Energy Buildings,
40(3):394–398, 2008.
[2] I. Sartori, A. Napolitano, and K. Voss. Net zero energy buildings: A consistent definition framework. Energy Buildings,
48(1):220–232, 2012.
[3] L. Deori, L. Giulioni, and M. Prandini. Optimal building climate control: A solution based on nested dynamic programming
and randomized optimization. In Proc. IEEE Conf. Decision Control, pages 4905–4910, Los Angeles, California, USA, 2014.
[4] A. Aswani, N. Master, J. Taneja, D. Culler, and C. Tomlin. Reducing transient and steady state electricity consumption in
HVAC using learning-based model-predictive control. Proc. IEEE, 100(1):240–253, 2012.
[5] M. Maasoumy, M. Razmara, M. Shahbakhti, and A. Vincentelli. Handling model uncertainty in model predictive control
for energy efficient buildings. Energy Buildings, 77(1):377–392, 2014.
[6] F. Oldewurtel, D. Sturzenegger, and M. Morari. Importance of occupancy information for building climate control. Appl.
Energy, 101(1):521–532, 2013.
[7] F. Oldewurtel, A. Parisio, C. N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, and M. Morari. Use of model
predictive control and weather forecasts for energy efficient building climate control. Energy Buildings, 45(1):15–27, 2012.
[8] J. Siroky, F. Oldewurtel, J. Cigler, and S. Privara. Experimental analysis of model predictive control for an energy efficient
building heating system. Appl. Energy, 88(9):3079–3087, 2011.
[9] D. Sturzenegger, D. Gyalistras, M. Morari, and R. S. Smith. Model predictive climate control of a Swiss office building:
Implementation, results and cost-benefit analysis. IEEE Trans. Control Syst. Technol., 24(1):1–12, 2016.
[10] M. Geidl, G. Koeppel, P. Favre-Perrod, B. Klockl, G. Andersson, and K. Frohlich. Energy hubs for the future. IEEE Power
Energy Mag., 5(1):24–30, 2007.
[11] D. Crawley, L. Lawrie, C. Pedersen, and F. Winkelmann. Energy plus: Energy simulation program. ASHRAE J., 42(4):49–
56, 2000.
[12] R. Evins, K. Orehounig, V. Dorer, and J. Carmeliet. New formulations of the energy hub model to address operational
constraints. Energy, 73(1):387–398, 2014.
[13] E. Fabrizio, M. Filippi, and J. Virgone. An hourly modelling framework for the assessment of energy sources exploitation
and energy converters selection and sizing in buildings. Energy Buildings, 41(10):1037–1050, 2009.
[14] A. Parisio, C. Del Vecchio, and A. Vaccaro. A robust optimization approach to energy hub management. Int. J. Elec. Power
Energy Sys., 42(1):98–104, 2012.
[15] D. Bertsimas, D. Brown, and C. Caramanis. Theory and applications of robust optimization. SIAM Rev., 53(3):464–501,
2011.
[16] X. Zhang, S. Grammatico, G. Schildbach, P. Goulart, and J. Lygeros. On the sample size of randomized mpc for chanceconstrained systems with application to building climate control. In Proc. Eur. Control Conf., pages 478–483, Strasbourg,
France, 2014.
[17] Jitka Dupačová, Nicole Gröwe-Kuska, and Werner Römisch. Scenario reduction in stochastic programming. Math. Program.,
95(3):493–511, 2003.
[18] G. Calafiore and M. Campi. The scenario approach to robust control design. IEEE Trans. Autom. Control, 51(5):742–753,
2006.
[19] A. Shapiro. Analysis of stochastic dual dynamic programming method. Eur. J. Oper. Res., 209(1):63–72, 2011.
[20] A. Nemirovski and A. Shapiro. Convex approximations of chance constrained programs. SIAM J. Optim., 17(4):969–996,
2006.
[21] W. Wiesemann, D. Kuhn, and M. Sim. Distributionally robust convex optimization. Oper. Res., 62(6):1358–1376, 2014.

22

Darivianakis, Georghiou, Smith, Lygeros: Data-Driven Robust MPC for Energy Efficient Buildings and Districts
Article submitted to IEEE Transactions on Control Systems Technology

[22] B. Van Parys, D. Kuhn, P. Goulart, and M. Morari. Distributionally robust control of constrained stochastic systems. IEEE
Trans. Autom. Control, 61(2):430–442, 2016.
[23] G. Lanckriet, L. El Ghaoui, C. Bhattacharyya, and M. Jordan. A robust minimax approach to classification. J. Mach.
Learn. Res., 3(1):555–582, 2003.
[24] D. Bertsimas, V. Gupta, and N. Kallus. Data-driven robust optimization. arXiv:1401.0212, 2013.
[25] E. Delage and Y. Ye. Distributionally robust optimization under moment uncertainty with application to data-driven
problems. Oper. Res., 58(3):595–612, 2010.
[26] H. Wolkowicz, R. Saigal, and L. Vandenberghe. Handbook of semidefinite programming: Theory, Algorithms, and Applications. Springer Science & Business Media, 2012.
[27] G. Darivianakis, A. Georghiou, R. S. Smith, and J. Lygeros. A stochastic optimization approach to cooperative building
energy management via an energy hub. In Proc. IEEE Conf. Decision Control, pages 7814–7819, Osaka, Japan, 2015.
[28] D. Sturzenegger, D. Gyalistras, V. Semeraro, M. Morari, and R. S. Smith. BRCM Matlab Toolbox: Model generation for
model predictive building control. In Proc. Amer. Control Conf., pages 1063–1069, Portland, USA, 2014.
[29] Indoor environmental input parameters for design and assessment of energy performance of buildings addressing indoor air
quality, thermal environment, lighting and acoustics. European committee for standardization, 2007.
[30] J. Steppeler, G. Doms, U. Schättler, H.W. Bitzer, A. Gassmann, U. Damrath, and G. Gregoric. Meso-gamma scale forecasts
using the nonhydrostatic model LM. Meteor. Atmos. Phys., 82(1-4):75–96, 2003.
[31] E. Lehmann and J. Romano. Testing Statistical Hypotheses. Springer, 2006.
[32] A. Georghiou, W. Wiesemann, and D. Kuhn. Generalized decision rule approximations for stochastic programming via
liftings. Math. Program., 152(1):301–338, 2014.
[33] E. Kerrigan and J. Maciejowski. Soft constraints and exact penalty functions in model predictive control. In Proc. UKACC
Int. Conf. Control, Cambridge, UK, 2000.
[34] G. Calafiore and L. El Ghaoui. On distributionally robust chance-constrained linear programs. J. Optim. Theory App.,
130(1):1–22, 2006.
[35] A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski. Adjustable robust solutions of uncertain linear programs. Math.
Program., 99(2):351–376, 2004.
[36] P. Benner and T. Damm. Lyapunov equations, energy functionals, and model order reduction of bilinear and stochastic
systems. SIAM J. Control Optim., 49(2):686–711, 2011.
[37] D. Gyalistras, A. Fischlin, M. Morari, C. Jones, F. Oldewurtel, A. Parisio, F. Ullmann, C. Sagerschnig, and A. Gruner.
Final report: Use of weather and occupancy forecasts for optimal building climate control (OptiControl). Technical report,
ETH Zurich, 2010.
[38] N. Razali and Y. Wah. Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests. J.
Stat. Model. Analyt., 2(1):21–33, 2011.
[39] St. Jewson and R. Caballero. Seasonality in the statistics of surface air temperature and the pricing of weather derivatives.
Meteorol. App., 10(04):367–376, 2003.
[40] J. Page, D. Robinson, N. Morel, and J. L. Scartezzini. A generalised stochastic model for the simulation of occupant presence.
Energy Buildings, 40(2):83–98, 2008.
[41] David L Wallace. Asymptotic approximations to distributions. Ann. Math. Stat., 29(3):635–654, 1958.
[42] H. Fakham, P. Degobert, and B. François. Control system and power management for a PV based generation unit including
batteries. In Int. Aegean Conf. Elect. Mach. Power Electron., pages 141–146, Bodrum, Turkey, 2007.
[43] E. I. Vrettos and S. A. Papathanassiou. Operating policy and optimal sizing of a high penetration RES-BESS system for
small isolated grids. IEEE Trans. Energy Convers., 26(3):744–756, 2011.

