1

Decentralized and Distributed Temperature Control
via HVAC Systems in Energy Efficient Buildings

arXiv:1702.03308v1 [cs.SY] 10 Feb 2017

Xuan Zhang, Wenbo Shi, Bin Yan, Ali Malkawi and Na Li

Abstract—In this paper, we design real-time decentralized and
distributed control schemes for Heating Ventilation and Air
Conditioning (HVAC) systems in energy efficient buildings. The
control schemes balance user comfort and energy saving, and
are implemented without measuring or predicting exogenous
disturbances. Firstly, we introduce a thermal dynamic model of
building systems and formulate a steady-state resource allocation
problem, which aims to minimize the aggregate deviation between
zone temperatures and their set points, as well as the building
energy consumption, subject to practical operating constraints,
by adjusting zone flow rates. Because this problem is nonconvex,
we propose two methods to (approximately) solve it and to
design the real-time control. In the first method, we present a
convex relaxation approach to solve an approximate version of
the steady-state optimization problem, where the heat transfer
between neighboring zones is ignored. We prove the tightness of
the relaxation and develop a real-time decentralized algorithm to
regulate the zone flow rate. In the second method, we introduce a
mild assumption under which the original optimization problem
becomes convex, and then a real-time distributed algorithm
is developed to regulate the zone flow rate. In both cases,
the thermal dynamics can be driven to equilibria which are
optimal solutions to those associated steady-state optimization
problems. Finally, numerical examples are provided to illustrate
the effectiveness of the designed control schemes.
Index Terms—Temperature control, HVAC systems, convex
relaxation, distributed/decentralized control, gradient algorithms.

I. I NTRODUCTION
It is reported that buildings are responsible for 40% of energy consumption, 70% of electricity consumption, and result
in 30% of greenhouse gas emission [1]. Roughly speaking,
Heating Ventilation and Air Conditioning (HVAC) systems in
buildings account for 40% of the energy use [2]. Therefore,
making HVAC systems more energy efficient is urgent for
improving environmental sustainability.
To date, various control techniques have been developed for
HVAC systems, including gain scheduling, optimal control,
robust control, nonlinear adaptive control, Model Predictive
Control (MPC), intelligent control based on artificial neural
network, fuzzy logic, and genetic algorithm, and so forth [3]–
[7]. Compared with the conventional on/off plus ProportionalSome of the material in this paper appeared in the Proceedings of the 2016
IEEE Global Conference on Signal and Information Processing, December
7-9, 2016, Arlington, VA, USA.
This work was supported by NSF ECCS 1548204, 1608509, NSF CAREER 1553407, and Harvard Center for Green Buildings and Cities. X.
Zhang, W. Shi, B. Yan and A. Malkawi are with Harvard Center for Green
Buildings and Cities, 20 Sumner Road, Cambridge, MA 02138, USA (email:
xuan zhang@g.harvard.edu, {wshi, byan, amalkawi}@gsd.harvard.edu). X.
Zhang and N. Li are with the School of Engineering and Applied Sciences,
Harvard University, 29 Oxford Street, Cambridge, MA 02138, USA (email:
nali@seas.harvard.edu).

Integral-Derivative (PID) control, these techniques are more
robust and energy efficient. However, most of them requires
centralized operation with heavy burdens of sensing, communication and computation, leading to much higher implementation cost than the conventional on/off plus PID control.
On the other hand, smart sensing, communication, computing, and actuation technologies have been stimulating the
emergence of distributed/decentralized control in network systems, including the smart grid [8]–[10], smart buildings [11],
[12], mobile robots [13], and intelligent transportation systems [14]. The advantages of distributed/decentralized control
include: good scalability as the network grows; reduction
in measurement, communication and computation compared
with centralized control; privacy preserving. Thus, applying
distributed/decentralized control to HVAC system control and
optimization becomes an area of active research. Representative work includes, for example, distributed MPC [15]–[17], as
well as mean-field game based distributed control [18]. Among
these popular control mechanisms, distributed MPC is the
most promising one, especially for large buildings. However,
it still requires a large amount of sensing, communication and
computation. And in most cases, it needs good prediction of
future disturbances, i.e., outdoor temperature, sunlight, indoor
occupancy, etc., which may be hard to obtain.
Contribution of this paper. This paper aims to develop
real-time control schemes for HVAC systems in commercial
buildings. Specifically, we aim to design decentralized and
distributed algorithms to guide each thermal zone to properly
adapt their supply air flow rates such that system-wide objectives are achieved under given operating conditions. The proposed control schemes (i) are scalable with respect to building
structures, (ii) satisfy the system operating constraints, and (iii)
ensure system efficiency, reliability and user comfort. Different
from literature [15]–[18], the control schemes designed in this
paper are based on solving steady-state resource allocation
problems via gradient algorithms – they are dynamic feedback
controllers that can be implemented without measuring or
predicting disturbances, which are different from controllers
based on MPC [5]–[7], [15]–[17].
To begin with, we provide the detailed problem setup,
including an introduction of the HVAC system configuration,
a commonly used thermal dynamic model of the building
network, and an optimization problem that aims to minimize a
weighted sum of the aggregate deviation between zone temperatures and their set points and the building energy consumption
subject to practical operating constraints (Section II). Since
the original optimization problem is nonconvex, two different
methods are proposed to (approximately) solve it. Firstly, we
present a convex relaxation approach in which the heat transfer
between neighboring zones is ignored (Section III). We show

2

that this relaxation is always tight under a mild condition, and
develop a decentralized control scheme for real-time zone flow
rate regulation. Also, we extend the proposed method to HVAC
system management in communities/neighborhoods. Secondly,
we introduce a mild assumption under which the original
optimization problem is naturally reformulated as a convex
one (Section IV). Then a distributed algorithm is developed
for real-time zone flow rate regulation. In both scenarios, the
thermal dynamics can be driven to equilibria which are the
optimal solutions of those associated steady-state optimization
problems. Lastly, two numerical examples are provided to
illustrate the effectiveness of the designed control schemes
(Section V), using a building with four adjacent zones.
Notation: ẋ(t) denotes the derivative of a state variable x(t)
d
with respect to time t, i.e., ẋ(t) = dt
x(t). x  ()y denotes
that x is much greater (less) than y. The positive projection
of a function h(y) on a variable x ∈ [0, +∞), (h(y))+
x is:

h(y)
if x > 0
(h(y))+
.
x =
max(0, h(y)) if x = 0
II. P ROBLEM S ETUP
A. HVAC system in buildings
The schematic of a typical HVAC system is illustrated in
Figure 1, which consists of an Air Handling Unit (AHU) for
the whole building and a set of pressure independent Variable
Air Volume (VAV) boxes for each zone [19]. The AHU is
equipped with dampers, a cooling/heating coil, and a Variable
Frequency Drive (VFD) fan: the dampers mix the return air
from the building with the outside air; the cooling/heating coil
cools down/heats up the mixed air; the VFD fan adjusts its own
speed based on the total air flow rate/opening controlled by
VAV boxes to keep the duct pressure at a certain level, and
drives the cooled/heated air to each zone. Each VAV box has a
damper to control the air flow rate supplied to the zone, which
is considered as the only controllable input to the system in
this paper. Also, we will focus on optimizing HVAC operation
when the system is in either cooling or heating mode. In the
cooling mode, the AHU supply air temperature is often set to
12.8◦ C which generally provides the required humidity ratio
to maintain space thermal comfort [20]. In the heating mode,
usually dehumidification is not necessary. Therefore, humidity
control is not considered in this paper.
B. Thermal dynamic model
According to the above configuration, we model a building
as an undirected connected graph (N , E). Here N is the set
of nodes representing zones/rooms, and E ⊆ N × N is the
set of edges. An edge (i, j) ∈ E means that zones i and j
are neighbors. Let N (i) denote the set of neighboring zones
of zone i. The thermal dynamics for the building is described
by a reduced Resistance-Capacitance (RC) model [21] (more
discussion on the model is available in Remark 1):
X Tj − Ti
T o − Ti
Ci Ṫi =
+
+ ca mi (T s − Ti ) + Qi
Ri
Rij
j∈N (i)

i∈N

Fig. 1: Schematic of a typical AHU&VAV system.
resistance of the wall and window separating zone i and
outside, Rij is the thermal resistance of the wall separating
zones i and j, ca is the specific heat of the air, mi is the
flow rate of the supply air, T s is the supply air temperature which is usually a constant [17], and Qi ≥ 0 is the
heat gain from exogenous sources (e.g., user activity, solar
radiation and device operation). If Ri , Rij , Ci are not known
from design specification, they can be obtained via model
identification [22], [23]. In this paper, we consider the flow
rate mi as the only control input to each zone, which is a
common practice in others’ work as well [17], [24], [25].
Proposition 1. When the HVAC system is off, i.e., mi = 0,
system (1) asymptotically converges to an equilibrium point
which is uniquely determined by the inputs T o , Qi . When
the HVAC system is on, the asymptotic convergence property
of system (1) remains and the equilibrium point is uniquely
determined by the inputs T o , mi , Qi .
This proposition can be directly derived by rearranging (1)
in state-space representation, and showing that the system
matrix is Hurwitz. Therefore, the desiderata is to design
mi only for periods when the HVAC system is on, more
specifically, is to design the dynamics of mi to drive (1) to
some desired state.
C. The optimization problem
In reality, each zone has a desired temperature which is the
set point determined by users. The objective of the HVAC
control considered in this paper is to regulate the temperature
to be close to the set point in each zone, and to minimize the
total energy consumption in the building. More specifically,
we consider the following steady-state optimization problem:
X1
wX
min
ri (Zi − Tiset )2 +
ca mi |Zi − T s |
Zi ,mi
2
η
i∈N
i∈N
X
+ws(
m i )3
(2a)
i∈N
o
X
T − Zi
Zj − Zi
s. t.
+
+ ca mi (T s − Zi ) + Qi = 0
Ri
Rij
j∈N (i)

(2b)

(1)

where Ci is the thermal capacitance, Ti is the indoor temperature, T o is the outdoor temperature, Ri is the thermal

Timin
mmin
i

≤ Zi ≤ Timax
≤ mi ≤ mmax
i

(2c)
(2d)

3

X
i∈N

mi ≤ m

(2e)

where i ∈ N in (2b)-(2d), ri , w are positive weight coefficients, η is a given constant named Coefficient of Performance
(COP, i.e., the ratio of the produced cold/heat to the consumed
energy) of the cooling/heating coil, s is a given coefficient
regarding the fan power consumption, Tiset is the temperature
set point satisfying Timin < Tiset < Timax , [Timin , Timax ] is
] is the range of mi in
, mmax
the user comfort range, [mmin
i
i
min
which mi
is usually close to zero to guarantee a minimal
ventilation level [16], and m is the upper bound of the total
flow rate in the building. Note that (i) to avoid confusion
between steady-state values and temperature dynamics, we use
Zi to denote the steady-state temperature value whereas Ti is
the temperature in the dynamic model (1), and (ii) T o and
Qi are exogenous disturbances. We assume that problem (2)
is feasible and satisfies Slater’s condition [26]. Moreover, we
have four important remarks for problem (2).
• In (2a), the term regarding the building energy consumption
contains two parts, i.e., the energy
consumption of the coolP
ing/heating coil given by η1 i∈N ca mi |Zi − T s | (weighted
by w), and the power consumption of the supply fan which is
approximately proportional
to the cube of the total supply air
P
flow [17], i.e., s( i∈N mi )3 (weighted by w).
• Parameters ri , w are determined by users: if preferring more
comfort, they can increase ri and decrease w, and vice versa.
• In the cooling mode, T s  Zi (or Ti ), ∀i; in the heating
mode, T s  Zi (or Ti ), ∀i. This is usually true in practice.
For example [17], in the cooling mode, T s = 12.8◦ C while
Ti ≥ 21.5◦ C. Note that once the mode is determined, the sign
of Zi − T s (or Ti − T s ) is determined.
• In (2e), we do not impose any lower bound constraint
for
flow rate as theP
lower bound usually equals
P the total
min
. In addition, m < i∈N mmax
holds, otherwise,
i
i∈N mi
this constraint would become redundant.
To conclude, the goal is to design the regulating rule for mi
in a decentralized or distributed way so that system (1) reaches
a steady state which is an optimal solution to problem (2).
Remark 1. Though system (1) is a 1st-order RC model, using
higher order RC models does not affect the formulation of (2)
since it is a steady-state optimization problem. For example,
for the 2nd-order model in [17], [21]
X Tij − Ti
T o − Ti
+
+ ca mi (T s − Ti ) + Qi
Ci Ṫi =
Ri
Rij
j∈N (i)

Cij Ṫij =

Ti − Tij
Tj − Tij
+
Rij
Rij

where Tij is the temperature of the wall separating zones i
and j, and Cij is the thermal capacitance of the wall, the
corresponding steady-state Equation (2b) is given by
X Zj − Zi
T o − Zi
+
+ ca mi (T s − Zi ) + Qi = 0
Ri
2Rij
j∈N (i)

which results from the steady-state equation Zij = (Zi +Zj )/2
(Zij is the steady-state temperature of the wall), i.e., the
steady-state equation of the higher order model can be reduced
to formulation (2b) by eliminating states of solids in the
building envelope. Because our control design procedures proposed later are based on solving the steady-state optimization

problem (2), using higher order models will not affect them.
III. M ETHOD I: A N A PPROXIMATE S OLUTION P ROCEDURE
A. Approximation and convex relaxation
Problem (2) is nonconvex because the quadratic equality
constraint (2b) is nonconvex in mi , Zi . In reality, due to
that Zi , Zj (or Ti , Tj ) of neighboring zones are often very
close to each other and Rij is not small, the total heat gain
from neighboring zones is (sometimes much) less dominant
compared with the heat gain from the outside plus the indoor
heat gain P
in every zone. Therefore, in this section,
P we ignore
Zj −Zi
Tj −Ti
in
(1)
as
well
as
the term
j∈N (i) Rij
j∈N (i) Rij
in (2b). In addition, we approximate the power consumption
P
of the supply fan by one of its upper bound sφ i∈N 12 m2i
where φ is a given constant satisfying φ ≥ m, because using
a Cauchy-Schwarz inequality, we have
P
P
X1
sφ( i∈N mi )2
s( i∈N mi )3
2
m ≥
≥
sφ
2 i
2|N |
2|N |
i∈N

where |N | is the number of zones. By doing this, the objective function of the optimization problem becomes separable,
which will enable us to design a decentralized controller later.
To summarize, we obtain an approximate model to (1) as
T o − Ti
+ ca mi (T s − Ti ) + Qi
(3)
Ri
together with an approximate problem to (2) given by

X 1
w
1
ri (Zi − Tiset )2 + ca mi |Zi − T s | + wsφm2i
min
Zi ,mi
2
η
2
i∈N
(4a)
T o − Zi
s
s. t.
+ ca mi (T − Zi ) + Qi = 0
(4b)
Ri
min
max
Ti
≤ Zi ≤ Ti
(4c)
Ci Ṫi =

mmin
≤ mi ≤ mmax
i
i
X
mi ≤ m

(4d)
(4e)

i∈N

where i ∈ N in (3), (4b)-(4d). Remark that the approximate
model (3) still keeps the global asymptomatic stability as
in Proposition 1. We assume that problem (4) is feasible
and satisfies Slater’s condition, moreover, we require ri >
wc2a
sφη 2 , i ∈ N so that (4a) can be strictly convex in Zi , mi .
Although (4) is still nonconvex due to the term mi Zi in (4b),
we can actually convexify it and show that theo convexification
is indeed exact. From (4b), we have mi =
o

∆

fi (Zi ) =

T −Zi
Ri

+ Qi

ca (Zi − T s )

T −Zi
Ri

+Qi

ca (Zi −T s ) .

Define

= mi ≥ mmin
> 0, i ∈ N
i

where the physical meaning of fi (Zi ) is the approximate air
flow rate for zone i to stay at temperature Zi . Then
fi0 (Zi )
fi00 (Zi ) =

=

T s −T o
Ri

− Qi

ca (Zi − T s )2
s
o
−2( T R−T
− Qi )
i
ca (Zi − T s )3

, i∈N
> 0, i ∈ N

4

hold, i.e., fi (Zi ) is convex in Zi , because
T s  T o holds in
s
o T s −T o
the cooling mode, and T  T , Ri > Qi hold (as Qi is
less dominant [16], [18], especially in winter) in the heating
mode. Therefore, we can relax (4) as a convex optimization
problem given by

X 1
w
1
set 2
s
2
min
ri (Zi − Ti ) + ca mi |Zi − T | + wsφmi
Zi ,mi
2
η
2
i∈N
(5a)
T o −Zi
Ri

+ Qi

≤ mi , i ∈ N
ca (Zi − T s )
Timin ≤ Zi ≤ Timax , i ∈ N

(5b)

mmin
≤ mi ≤ mmax
, i∈N
i
i
X
mi ≤ m.

(5d)

s. t.

i∈N

(5c)
(5e)

Remark 2. In the cooling mode, fi0 (Zi ) < 0, ∀i holds while
in the heating mode, fi0 (Zi ) > 0, ∀i holds, i.e., fi (Zi ) is
monotonic in Zi , ∀i.

B. Tightness of the convex relaxation
Before investigating the tightness of the convex relaxation,
we make an assumption on the system parameters.
Assumption 1. In the cooling mode, Tiset < T o , ∀i holds; in
the heating mode, Tiset > T o , ∀i holds. Moreover, the zone
temperature set point satisfies fi (Tiset ) ≥ mmin
, ∀i.
i
This assumption holds in usual practice. Firstly, when the
outdoor temperature is higher (in hot days)/lower (in cold
days) than users’ expectation (quantified as Tiset ), the HVAC
system should be turned on. So we usually have Tiset < T o , ∀i
for cooling and Tiset > T o , ∀i for heating. Secondly, note that
mmin
≈ 0, ∀i holds as stated earlier, and fi (Zi ) stands for
i
the approximate flow rate for zone i to stay at temperature
Zi . The inequality fi (Tiset ) ≥ mmin
means that in order to
i
make the temperature of zone i be the set point, the flow rate
needs to be no smaller than its minimum. Otherwise, we will
have fi (Tiset ) < mmin
≈ 0 which can lead to Tiset ≥ T o
i
in the cooling mode, and Tiset ≤ T o in the heating mode, a
contradiction to the first part of this assumption. To sum up,
users can choose their set points to satisfy this assumption in
reality, which is completely decentralized.
The Lagrangian of problem (5) is

X 1
w
1
L=
ri (Zi − Tiset )2 + ca mi |Zi − T s | + wsφm2i
2
η
2
i∈N
!
o
T −Zi
X

X
+ Qi
Ri
+
+
ζi
−
m
+
λ
m
−
m
i
i
ca (Zi − T s )
i∈N
i∈N
X
X
+
νi+ (Zi − Timax ) +
νi− (Timin − Zi )
i∈N

+

X

i∈N

µ+
i (mi

−

mmax
)
i

+

i∈N
−
+
ζi , νi+ , νi− , µ+
i , µi , λ

X

min
µ−
− mi )
i (mi

i∈N

where
are the Lagrange multipliers
(dual variables) for constraints (5b)-(5e). Since problem (5) is
convex, feasible and satisfies Slater’s condition, the Karush-

Kuhn-Tucker (KKT) conditions are necessary and sufficient
conditions for optimality [26], given by
s

o

T −T
− Qi
∂L
w
= ri (Zi − Tiset ) ± ca mi + ζi Ri
∂Zi
η
ca (Zi − T s )2
+
+νi − νi− = 0
(6a)
∂L
w
−
+
= wsφmi ± ca (Zi − T s ) − ζi + µ+
i − µi + λ = 0
∂mi
η
(6b)
ζi (fi (Zi ) − mi ) = 0, ζi ≥ 0, fi (Zi ) − mi ≤ 0
(6c)

νi+ (Zi − Timax ) = 0, νi+ ≥ 0, Zi − Timax ≤ 0
νi− (Timin − Zi ) = 0,
max
µ+
) = 0,
i (mi − mi
−
min
µi (mi − mi ) = 0,

νi−
µ+
i
µ−
i

≥

0, Timin

− Zi ≤ 0

≥ 0, mi − mmax
≤0
i
≥

0, mmin
i

(6d)
(6e)
(6f)

− mi ≤ 0

(6g)

mi − m ≤ 0

(6h)

!
λ+

X

mi − m

= 0, λ+ ≥ 0,

i∈N

X
i∈N

where i ∈ N in (6a)-(6g), and “±” takes “+” in the cooling
mode and “−” in the heating mode in (6a)-(6b).
Remark 3. The convex relaxation from problem (4) to (5)
is tight if and only if any solution of (6) satisfies ζi > 0 or
ζi = 0, fi (Zi ) = mi , ∀i.
As for the tightness, we have the following theorem.
Theorem 1. Under Assumption 1, the convex relaxation from
problem (4) to (5) is always tight.
P roof : Let us assume that there exists an i in the solution
of (6) such that ζi = 0, fi (Zi ) < mi holds. Then we can
obtain ri (Zi − Tiset ) + wη ca mi = νi− − νi+ in the cooling
mode and ri (Zi − Tiset ) − wη ca mi = νi− − νi+ in the heating
+
w
s
mode from (6a), and µ−
i = wsφmi ± η ca (Zi − T ) + µi +
λ+ > 0 from (6b). When νi− = νi+ = 0 holds, we have
Zi < Tiset in the cooling mode and Zi > Tiset in the heating
mode. When νi− > νi+ = 0 holds which only happens in
the cooling mode (otherwise in the heating mode, we end up
with Zi = Timin > Tiset , a contradiction), we have Zi =
Timin < Tiset . When νi+ > νi− = 0 holds which only happens
in the heating mode, we have Zi = Timax > Tiset . These facts
lead to fi (Tiset ) < fi (Zi ) < mi = mmin
by Remark 2, a
i
contradiction to Assumption 1. Based on Remark 3, the convex
relaxation from problem (4) to (5) is always tight.

C. A decentralized algorithm
Theorem 1 indicates that solving (4) is equivalent to solving (5), while (5) can be solved in either a centralized or distributed/decentralized way. Any centralized algorithm requires
to measure the outdoor temperature T o and the indoor heat
gain Qi in every zone. Because these exogenous disturbances
can fluctuate frequently and are not easy to obtain, the cost
of centralized algorithms would be expensive. In this section,
we develop a real-time decentralized algorithm that does not
need measurement of these exogenous disturbances.
Since (5) is convex, feasible and satisfies Slater’s condition,
we design the following algorithm to solve (5) based on a

5

standard primal-dual gradient method [27]:
!
T s −T o
− Qi
∂L
set
Żi = − kZi
= kZi ri (Ti − Zi ) − ζi Ri
∂Zi
ca (Zi − T s )2
!
w
∓ ca mi − νi+ + νi−
(7a)
η
!
∂L
w
ṁi = − kmi
= kmi (−wsφmi ∓ ca (Zi − T s ) + ζi
∂mi
η

ζ̇i
ν̇i+
ν̇i−
µ̇+
i
µ̇−
i
λ̇+

−
+
− µ+
i + µi − λ )
!+
!+
T o −Zi
+ Qi
∂L
Ri
= k ζi
− mi
=kζi
∂ζi
ca (Zi − T s )
ζi
ζi
!+
∂L
=kν +
= kν + (Zi − Timax )+
νi+
i
i
∂νi+ +
νi
!+
∂L
= kν − (Timin − Zi )+
=kν −
νi−
i
i
∂νi− −
νi
!+
∂L
=kµ+
= kµ+ (mi − mmax
)+
i
µ+
i
i
i
∂µ+
+
i
µi
!+
∂L
= kµ− (mmin
− m i )+
=kµ−
i
µ−
i
i
i
∂µ−
i
µ−
i
!+
!+
X
∂L
+
=kλ+
=
k
mi − m
λ
∂λ+ +
+
λ

i∈N

(7b)
(7c)

(7d)

(7e)

(7f)

(7g)

(7h)

Fig. 2: Information exchange of the decentralized controller.
This theorem requires T o , Qi to be either constant, step
change, or slow-varying, which holds in practice as these
disturbances vary at a time-scale of minutes. Remark that our
controller operates in real-time, i.e., at a time-scale of seconds.
In (7a) and (7c), the disturbances T o , Qi appear. Motivated
by [29], to make the algorithm implementable without measuring these terms, we substitute (3) into (7a) and (7c) to get

− ζi

λ

where i ∈ N in (7a)-(7g), kZi , kmi , kζi , kν + , kν − , kµ+ , kµ− ,
i
i
i
i
kλ+ are positive scalars representing the controller gains, and
“∓” takes “−” in the cooling mode and “+” in the heating
mode in (7a)-(7b). Note that Ti has its own dynamics given
by (3) and thus can not be designed, which is why we replace
Ti with Zi initially, i.e., Zi , i ∈ N are ancillary state variables.
According to [27], [28], it is true that (7) asymptotically
converges to an equilibrium point which is the unique optimal
solution of (5), since the optimization problem is convex
and the objective function is strictly convex in the decision
variables. Now using mi in (7b) as the input to system (3),
we can naturally obtain a real-time decentralized controller to
regulate (3) to a steady state which solves (4).
Theorem 2. Given constant/step change/slow-varying T o , Qi ,
under Assumption 1, the equilibrium point of the overall
system (3) and (7) is unique, asymptotically stable and Ti , mi
of the equilibrium point is the optimal solution of (4).
P roof : For constant/step change/slow-varying disturbances
T o , Qi , since (5) is convex and the objective function is strictly
convex in Zi , mi , its optimal solution is unique [26]. Therefore, the equilibrium point of (7) is unique, asymptotically
stable [27], [28] and is the optimal solution of (5)/(4) by
Theorem 1. So the equilibrium point of (3) and (7) is also
asymptotically stable, due to the cascade nature, i.e., (7)→(3).
Under the fact that the equilibrium point of (3) is uniquely determined by the inputs T o , mi , Qi where mi is given by (7b),
we have Ti = Zi (Ti is the state given by (3)) when the overall
system reaches steady state, which completes the proof.

w
ca mi − νi+ + νi−
η
!
− Ci Ṫi − ca mi (Ti − T s )

Żi =kZi ri (Tiset − Zi ) ∓

ζ̇i =kζi

T s −Ti
Ri

(8a)

ca (Zi − T s )2
Ti −Zi
Ri

+ Ci Ṫi + ca mi (Ti − T s )
ca (Zi − T s )

!+
− mi

(8b)
ζi

in which the derivative action Ṫi can be implemented by using
a differentiator with some form of filtering [30].
Implementation. The proposed control scheme (7b), (7d)(7h) and (8) is completely decentralized as shown in Figure 2
and can be implemented as follow. Given T s , Ri , Ci , ri , w,
φ, s, η, mmin
, mmax
, each zone in the building collects Tiset ,
i
i
min
max
[Ti , Ti
] from users, locally measures its indoor temperature Ti , receives the feedback signal λ+ from the supply
fan/duct, and then uses these information to update Zi , mi , ζi ,
−
νi+ , νi− , µ+
i , µi based on (7b), (7d)-(7g) and (8). On the
other hand, given P
m, the supply fan/duct locally measures
the total flow rate i∈N mi which is proportional to the fan
speed [24], and uses (7h) to update λ+ , then broadcasts it. Here
T s , Ri , Ci , s, η, mmin
, mmax
, m are building parameters and
i
i
set
min
max
Ti , [Ti , Ti
], ri , w, φ are parameters specified by users.
D. An extension
In this section, we provide an alternative scenario which
fits the above solution procedure. Consider a community or
a neighborhood consisting of a number of separated houses.
Each house is equipped with an HVAC system. Since these
houses are separated (i.e., Rij = ∞), we can naturally use
system (3) to model their temperature dynamics. The objective
is to regulate the temperature in each house to be close to its set
point, meanwhile, to minimize the total energy consumption
of the community. As a result, we end up with a steady-

6

(a)

(b)

(c)

Fig. 3: m1 , m2 as functions of Z1 , Z2 . (a): m1 . (b): m2 . (c): the constraint set/domain for Z1 , Z2 .
state optimization problem given by (4). Note that (i) each
entire house is modeled by a RC model (we do not specify
zones inside the house) and has a unique temperature set point
determined by users; (ii) mmin
= 0 holds in (4d) because we
i
now allow houses to turn off their HVAC systems; (iii) m
in (4e) is an indicator of the upper bound of the total energy
consumption in the community; (iv) the supply air temperature
can be different in houses, i.e., T s can be replaced by Tis
in (3) and (4a)-(4b) (the mixed mode case, i.e., some of HVAC
systems are in the cooling mode while the others are in the
heating mode, is naturally included), and similar for w, η, s, φ.
Correspondingly, we modify Assumption 1 a bit as follow.
Assumption 2. For any house i, when the HVAC system is
on, the house temperature set point satisfies fi (Tiset ) > 0, ∀i.
Otherwise, the HVAC in house i is turned off and all terms and
constraints relating to house i in problem (4)/(5) are excluded.
As before, users can choose their set points to satisfy this
assumption, which is decentralized. Under this assumption and
following the proof of Theorem 1, we can still ensure the
tightness of the convex relaxation procedure in Section III.A.
Moreover, we can still derive a decentralized algorithm similar
to (7)-(8) for HVAC system management in the community.
What needs to be emphasized is that now (7h) is updated in
a control center of the community, which receives mi from
each house and sends λ+ to them (when mi = 0, the control
center can know that house i has turned off its HVAC system).
IV. M ETHOD II: T HE G ENERAL C ASE
In this section, we focus on solving problem (2) directly.
The key step is to reformulate (2) as an optimization problem
without using the decision variable mi by substituting (2b)
into (2a) and (2d)-(2e) to eliminate mi . Before showing the
details, let us first consider a simple scenario, i.e., a building
with two adjacent zones, to gain an insight into this method.
A. An illustrating example
Let the constraint set of problem (2) be
30 − Z1
Z2 − Z1
+
+ 1.012m1 (12.8 − Z1 ) + 0.1 = 0
15
18
30 − Z2
Z1 − Z2
+
+ 1.012m2 (12.8 − Z2 ) + 0.2 = 0
16
18
0.01 ≤ m1 ≤ 0.5, 0.01 ≤ m2 ≤ 0.5, m1 + m2 ≤ 0.7.

Fig. 4: The constraint sets for Z1 , Z2 when the upper bound
of m1 + m2 changes.

where we do not add user comfort constraint (2c) (the reason
will become clear later). The domain for m1 , m2 , Z1 , Z2
is shown in Figure 3(a)-(b), which are nonconvex surfaces.
However, by rewriting m1 , m2 as functions of Z1 , Z2
m1 =

−Z1
+ Z218
+ 0.1
,
1.012(Z1 − 12.8)

30−Z1
15

m2 =

−Z2
+ Z118
+ 0.2
1.012(Z2 − 12.8)

30−Z2
16

and substituting them into all three inequality constraints, we
can see that the constraint set for Z1 , Z2 is actually convex as
shown in Figure 3(c), which is also the two-dimensional view
of Figure 3(a)-(b). The reason is as follows. Since the cooling
mode is considered (i.e., Z1 , Z2  T s = 12.8), the constraints
0.01 ≤ m1 ≤ 0.5 and 0.01 ≤ m2 ≤ 0.5 are equivalent to
30 − Z1
Z2 − Z1
+
+ 0.1
15
18
≤0.5 × 1.012(Z1 − 12.8)
30 − Z2
Z1 − Z2
0.01 × 1.012(Z2 − 12.8) ≤
+
+ 0.2
16
18
≤0.5 × 1.012(Z2 − 12.8)
0.01 × 1.012(Z1 − 12.8) ≤

which are linear, and therefore, convex. So the convexity of the
constraint set only depends on the convexity of the function
m1 + m2 . Then we calculate the Hessian matrix of
m1 + m2 =

−Z1
+ Z218
+ 0.1
+
1.012(Z1 − 12.8)

30−Z1
15

−Z2
+ Z118
+ 0.2
1.012(Z2 − 12.8)

30−Z2
16

7

with respect to Z1 , Z2 . Any point in the domain in Figure 3(c)
has a positive semi-definite Hessian matrix of m1 + m2 .
Next, we modify the upper bound of m1 + m2 to be 0.5
and 0.3, and draw the constraint sets of Z1 , Z2 , as illustrated
in Figure 4. It can be seen that these constraint sets are also
convex, because any point in the domains displayed in Figure 4
has a positive semi-definite Hessian matrix of m1 + m2 . Note
that including user comfort constraint (2c) will not affect the
above analysis since it does not affect the convexity of the
constraint set of Z1 , Z2 but only makes the set smaller, i.e.,
the above analysis is actually less conservative than that with
constraint (2c). These numerical results indicate that once we
exclude mi from problem (2), convexity can be obtained.
B. Solution procedure
Now we provide the solution procedure for (2). Define
P
o
Z −Z
i
X
X T R−Z
+ j∈N (i) jRij i + Qi
∆
i
h(Z) =
mi =
> 0.
ca (Zi − T s )
i∈N

i∈N

where Z is a collection of Zi , i ∈ N . By substituting (2b)
into (2a) and (2d)-(2e) to eliminate mi , problem (2) becomes

This fact is demonstrated in the previous example: even though
Z1 , Z2 deviate about 8◦ C from each other, the constraint set
remains convex as shown in Figures 3(c)-4.
Under Assumption 3 (and the sign of Zi −T s is determined),
problem (9) is convex since the objective function (9a) is
strictly convex in Zi as convexity is preserved from h(Z) > 0
to (h(Z))3 [26]. Motivated by a standard primal-dual gradient
method [27], we design the following algorithm to solve (9)
(assuming the cooling mode here, i.e., Zi > T s , the heating
mode case is similar):
 1

X 1
max
Żi =kZi ri (Tiset − Zi ) + µ+
+
+
m
c
a
i
i
Ri
Rij
j∈N (i)

 1

X µ+
X 1
j
−
− µ−
+
+ mmin
ca − νi+
i
i
Rij
Ri
Rij
j∈N (i)
j∈N (i)
P
T s −Z
T s −T o
X µ−
+ j∈N (i) Rij j − Qi
Ri
j
−
+ νi −
+
Rij
ca (Zi − T s )2
j∈N (i)
!
!
X
1
w
2
+
+
(3wsh (Z) + λ ) +
Rij ca (Zj − T s )
ηRi
j∈N (i)

(10a)
X1
w X T o − Zi
ri (Zi − Tiset )2 +
+ Qi
min
Zi
2
η
Ri
+ws(h(Z))3

mmin
≤
i

s. t. Timin ≤ Zi ≤ Timax
P
Z −Z
T o −Zi
+ j∈N (i) jRij i + Qi
Ri
ca (Zi − T s )
h(Z) ≤ m.

(9a)
(9b)

≤ mmax
i

(9d)

is positive semi-definite.
Assumption 3 is not conservative and it usually holds in
2
h(Z)
is usually diagonal dominant due to
practice, because ∂ ∂Z
2
that Zi , Zj of neighboring zones do not deviate too much from
each other (this can also be guaranteed by properly adjusting
user comfort range), i.e.,

j∈N (i)

(10c)

i

µ̇+
i =kµ+

(9c)

Assumption 3. The Hessian matrix of h(Z), i.e.,
 2

∂ h(Z)
∂ 2 h(Z)
2
2
∂ h(Z)  ∂Zi
∂Zi ∂Zj 
=
∂ 2 h(Z)
∂ 2 h(Z)
∂Z 2
∂Zj ∂Zi
∂Zj2
s
P
s
o
T
−Z
−2( T R−T
+ j∈N (i) Rij j − Qi )
∂ 2 h(Z)
i
=
>0
∂Zi2
ca (Zi − T s )3
1
1
∂ 2 h(Z)
=−
−
<0
∂Zi ∂Zj
Rij ca (Zi − T s )2
Rij ca (Zj − T s )2

∂ 2 h(Z)
≈
∂Zi ∂Zj

ν̇i− =kν − (Timin − Zi )+
ν−
i

where i ∈ N in (9b)-(9c). Constraint (9c) can actually become
linear once the cooling/heating mode (the sign of Zi − T s )
is determined (also, the absolute value sign in (9a) can be
removed). The convexity of the constraint set only depends on
constraint (9d). Now we claim that the HVAC system satisfies
the following assumption under normal operating conditions.

X
∂ 2 h(Z)
−
2
∂Zi

(10b)

i

i∈N

i∈N

ν̇i+ =kν + (Zi − Timax )+
ν+

s
o
−2( T R−T
i

− Qi )

ca (Zi − T s )3

> 0.

i

i

X Zj − Zi
T o − Zi
+
+ Qi
Ri
Rij
j∈N (i)
!+

− mmax
ca (Zi − T s )
i

(10d)
µ+
i

min
µ̇−
ca (Zi − T s ) −
i =kµ− mi
i

T o − Zi
Ri

!+
X Zj − Zi
−
− Qi
(10e)
Rij
j∈N (i)
µ−
i
!+
P
o
Zj −Zi
i
X T R−Z
+
j∈N (i) Rij + Qi
i
+
−m
λ̇ =kλ+
ca (Zi − T s )
i∈N
λ+
(10f)
−
+
where νi+ , νi− , µ+
i , µi , λ are the Lagrange multipliers (dual
variables) for constraints (9b)-(9d), kZi , kν + , kν − , kµ+ , kµ− ,
i
i
i
i
kλ+ are positive scalars representing the controller gains, and
i ∈ N in (10a)-(10e). In addition, we adopt the following low
pass dynamics for each mi as the control input to system (1):
!
P
Z −Z
T o −Zi
+ j∈N (i) jRij i + Qi
Ri
− mi
(11)
ṁi = km
ca (Zi − T s )

where km > 0 is the homogeneous controller gain. The reason
for using the low pass dynamics is that it can attenuate high
frequency noises to help improve system performance.
Theorem 3. Given constant/step change/slow-varying T o , Qi
(remark that they vary at a time-scale of minutes), under
Assumption 3, any trajectory of system (1), (10)-(11) asymptotically converges to a unique equilibrium point at which

8

problem (9)/(2) is solved and Ti = Zi , ∀i (this is because of
the strict convexity of the function (9a) in Zi and the cascade
nature of the overall system, i.e., (10)→(11)→(1)).
In (10)-(11), the disturbances T o , Qi appear. Motivated
by [29], to make the algorithm implementable without measuring these terms, we substitute (1) into (10)-(11) to get
w
− νi+ + νi−
ηRi
 1

X 1
X µ+
j
+
max
+ µi
+
+ mi ca −
Ri
Rij
Rij

Żi = kZi ri (Tiset − Zi ) +

j∈N (i)

j∈N (i)

 1

X 1
X µ−
j
min
− µ−
+
+
m
c
+
a
i
i
Ri
Rij
Rij
j∈N (i)

− (3wsh2 (Z) + λ+ )

j∈N (i)
T s −Ti
Ri

P

− Ci Ṫi +

j∈N (i)

1
)+
Rij ca (Zj − T s )

T s −Zj +Tj −Ti
Rij

ca (Zi −

− ca mi (Ti − T s ) !!

T s )2
(12a)

X Zj − Zi − Tj + Ti
Ti − Zi
+ Ci Ṫi +
Ri
Rij
j∈N (i)
!+

µ̇+
i = kµ+
i

+ ca mi (Ti − T s ) − mmax
ca (Zi − T s )
i

(12b)
µ+
i

min
µ̇−
ca (Zi − T s ) −
i = kµ− mi
i

−

Ti − Zi
− Ci Ṫi
Ri
!+

X Zj − Zi − Tj + Ti
− ca mi (Ti − T s )
Rij

j∈N (i)

m)+
λ+

+

λ̇ = kλ+ (h(Z) −
ṁi = km
Ti −Zi
Ri

(12c)
µ−
i

(12d)

P
j∈N (i)

Zj −Zi −Tj +Ti
Rij

+ ca mi (Ti − T s ) !

ca (Zi − T s )
(12e)
h(Z) =

X
i∈N

Ti −Zi
Ri

+ Ci Ṫi +

0.7
To

35
30
25
20

Q1
Q2

0.6

Q3

0.5

Q4

0.4
0.3

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

0.2

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

Fig. 6: Profiles of the disturbances (Q1 = Q2 ).
−
λ+ from the supply fan/duct and Tj , Zj , µ+
j , µj from its
neighboring zones, and then uses these information to up−
date Zi , νi+ , νi− , µ+
i , µi , mi based on (10b)-(10c), (12a)-(12c)
and (12e). On the other hand, given
P m, the supply fan/duct
locally measures the total flow rate i∈N mi which is proportional to the fan speed [24], and uses (12d) to update λ+ , and
then broadcasts 3wsh2 (Z)+λ+ . Note that the supply fan/duct
can compute λ+ and h(Z) via

λ̇+ = kλ+ (h(Z) − m)+
λ+
X
1 X
ṁi +
mi .
h(Z) =
km
i∈N

(13a)
(13b)

i∈N

Also, T s , Ri , Rij , Ci , s, η, mmin
, mmax
, m are building pai
i
set
min
max
rameters and Ti , [Ti , Ti
], ri , w are given by users.

− mi +

+ Ci Ṫi +

40

Power(kW)

X

Indoor heat gains

Outdoor temperature
Temperature(c−degree)

j∈N (i)

Fig. 5: Information exchange of the distributed controller.

P
j∈N (i)

Zj −Zi −Tj +Ti
Rij

+ ca mi (Ti − T s )

ca (Zi − T s )
(12f)
in which the derivative action Ṫi can be implemented by using
a differentiator with some form of filtering [30].
Implementation. The proposed controller (10b)-(10c)
and (12) is completely distributed as shown in Figure 5 and can
be implemented as follow. Given T s , Ri , Rij , Ci , ri , w, s, η,
mmin
, mmax
, each zone in the building collects Tiset ,
i
i
min
max
[Ti , Ti
] from users, locally measures its indoor temperature Ti , receives the feedback signals 3wsh2 (Z) +

V. N UMERICAL I NVESTIGATIONS
In this section, we present two numerical examples for
scenarios described in Sections III-IV respectively, using a
house with four adjacent zones as illustrated in Figure 1.
Only the cooling case is presented in the following, while the
heating case is similar under the proposed control schemes.
The parameters of the simulations are from [31], [32]: all
Ci = 20kJ/◦ C, all Ri = 15◦ C/kW, all Rij = 23◦ C/kW, ca =
1.012kJ/kg/◦ C, T s = 12.8◦ C (cooling), s = 2kW/(kg/s)3 , η =
2.9, all comfort ranges are within ±1.5◦ C of their set points,
all [mmin
, mmax
] = [0.01, 0.45]kg/s, m = 0.5kg/s, φ = 1kg/s,
i
i
all ri = 0.1p.u., w = 1p.u., all kZi = 0.067p.u., all kmi =
kζi = kν + = kν − = kµ+ = kµ− = kλ+ = km = 1p.u. (i.e.,
i
i
i
i
per unit), and the disturbance injection is shown in Figure 6.
The simulation result of the first scenario is illustrated in
Figures 7-9, in which we changed w = 1p.u. to w = 0.1p.u.
at 12h. The curves labelled with “app” indicate the case of
using (7b), (7d)-(7h) and (8) for the approximate model (3)
while the curves labelled without “app” are for the accurate

9

1

Zone 2

Tset1

22.5

22

21.5

21

0

2

4

6

appT2
Tset

2

22.5

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

Zone 3

Tset3

22.5
22

0

2

4

6

Temperature(c−degree)

Temperature(c−degree)

T3

23

1

appm2
m2

0.4

appm

3

0.3

m

3

0.2

appm

0.1

m4

0

4

0

2

4

6

8

Zone 4
24

appT3

23.5

1

m

23

22

10
12
14
Time(hour)

16

18

20

22

24

Ancillary variables

zeta1

T4

23

appsumm
summ

appzeta1

4
appT4

appzeta2

3

Tset4

zeta2
appzeta3

2

22

zeta3
1

appzeta4

21

zeta4
0

20

8 10 12 14 16 18 20 22 24
Time(hour)

appm

0.5

T2

23.5

8 10 12 14 16 18 20 22 24
Time(hour)

24

21.5

Flow rates
0.6

24

1

Temperature(c−degree)

Temperature(c−degree)

Zone 1

T

Flow rate(kg/s)

appT
23

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

Fig. 7: Temperatures in each zone under controller (7b), (7d)(7h) and (8): curves labeled with “app” indicate the case of
using the approximate model (3).

0

2

4

6

8

10
12
14
Time(hour)

16

18

20

22

24

applam
lam

Fig. 9: Air flow rates in each zone and ancillary variables
under controller (7b), (7d)-(7h) and (8): curves labeled with
“app” indicate the case of using the approximate model (3).

Fig. 8: Heat gain components in each zone by the accurate
model (1) under controller (7b), (7d)-(7h) and (8).

lows: all comfort ranges are within ±1.8◦ C of their set points,
[mmin
, mmax
] = [0.01, 0.15]kg/s, m = 0.5kg/s before 16h
4
4
while m decreases to 0.4kg/s after 16h, w = 0 before 8h and
w = 1p.u. after 8h (the other parameters remain the same).
The simulation result is illustrated in Figures 10-11. We can
see that the curves of the ancillary states Zi almost coincide
with the temperature curves Ti as expected. Before 8h, the
temperatures are exactly the same as their set points since
there is no consumption reduction purpose and the total flow
rate dose not saturate during this period. After 8h, deviations
appear due to the consideration of energy saving. The total
flow rate reaches its maximum around 15h, and from 16h to
18h, thus, the zone temperatures further deviate from their set
points (due to the lack of capacity). All these two scenarios
inspire us that tuning the weight coefficient w (or equivalently
ri , as the optimal solution of (2)/(4) depends on the ratio
ri /w, i ∈ N ) can balance user comfort and HVAC system
energy consumption – there always exists a tradeoff between
user comfort and energy saving.

model (1). It can be seen that the difference between these
two cases is very small due to that the heat transfer between
neighboring zones is negligible compared with the heat gain
from the other (outside and indoor) sources, as shown in
Figure 8 and explained in Section III as well. The temperature
deviations with respect to their set points before decreasing
the weight coefficient w are larger than those thereafter since
starting from 12h, user comfort becomes more important
while energy consumption saving becomes less. The total
flow rate reaches its maximum from 12h to 16h, while does
not saturate in other periods. Interestingly, there are always
deviations from the real temperatures to their set points (unless
w = 0), because there is a tradeoff between user comfort
and energy saving when using the proposed controller. All
ζi are always positive, indicating that the convex relaxation
procedure proposed in Section III always ensures tightness.
In the second scenario, we modify the parameters as fol-

VI. C ONCLUSION AND F UTURE W ORK
This paper presents decentralized/distributed control frameworks on real-time temperature regulation for HVAC systems
in energy efficient buildings. The proposed controllers adjust
air flow rates in each zone, which balances user comfort
and energy saving. Moreover, they can automatically adapt
to changes of disturbances such as the outdoor temperature
and indoor heat gains, without measuring or predicting those
values. Also, the implementation of the controllers are simple.
Future work includes: investigating, e.g., the H2 and H∞
performances of the controlled systems; considering the fact of
reheating by each VAV box, and natural ventilation scheduled
by users (these will be characterized as extra exogenous
inputs to the thermal dynamics); studying the interaction
between the controlled building HVAC systems and power
grids/microgrids; last but not least, extending the control
design frameworks to other heating and cooling systems, for
instance, ground source heat pump systems.

Zone 2
1.5

1

1
Power(kW)

Power(kW)

Zone 1
1.5

0.5

0

0

0

2

4

6

−0.5
8 10 12 14 16 18 20 22 24
0 2 4 6
Time(hour)
Outside and indoor heat gain
Heat
gain
from
neighbors
Zone 3

1.5

1.5

1

1
Power(kW)

Power(kW)

−0.5

0.5

0.5

0

−0.5

8 10 12 14 16 18 20 22 24
Time(hour)
Zone 4

0.5

0

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

−0.5

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

10

Zone 2
24

1
1

23

Tset

1

22.5
22
21.5
21
20.5

0

2

4

6

T

2

Z

Temperature(c−degree)

Temperature(c−degree)

Zone 1

T

23.5

2

Tset

2

23

[10]

22

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

Zone 3

Z3

23

Tset3

22.5

22

21.5

0

2

4

6

[11]

Zone 4
24

T3
Temperature(c−degree)

Temperature(c−degree)

23.5

[9]

22.5

21.5

8 10 12 14 16 18 20 22 24
Time(hour)

[8]

Z

23.5

Z4

23

[12]

Tset4

22

[13]

21

20

8 10 12 14 16 18 20 22 24
Time(hour)

T4

0

2

4

6

8 10 12 14 16 18 20 22 24
Time(hour)

[14]

Fig. 10: Temperatures in each zone under (10b)-(10c) and (12).
[15]
Flow rates

m

Flow rate(kg/s)

m2
m3

0.4

m4

2

[17]

summ

0.2
0

Power(kW)

[16]

1

0.6

0

2

4

6

8

10

12

14

16

18

20

22

24

[18]

1.5
Total consumption
Coil consumption
Fan consumption

1
0.5
0
0.5

0

0

2

4

6

8

10

12

14

16

[19]
[20]
18

20

22

24
mu+4

0

2

4

6

8

10

12

14

16

18

20

22

[21]
24

4
lambda

2
0

0

2

4

6

8

10

12
14
Time(hour)

16

18

20

22

24

[22]
[23]

Fig. 11: Air flow rates, the total consumption and ancillary
variables under (10b)-(10c) and (12).

[24]

R EFERENCES

[25]

[1] UNEP Sustainable Buildings & Climate Initiative, “Buildings and climate change: Summary for decision-makers,” Paris CEDEX 09, France:
Sustainable United Nations, 2009.
[2] U.S. Energy Information Administration Office of Integrated and International Energy Analysis, U.S. Department of Energy, “Annual energy
outlook 2012,” DOE/EIA-0383, 2012.
[3] D. S. Naidu and C. G. Rieger, “Advanced control strategies for heating,
ventilation, air-conditioning, and refrigeration systems – an overview:
Part I: Hard control,” HVAC&R Research, vol. 17, no. 1, pp. 2–21,
2011.
[4] ——, “Advanced control strategies for heating, ventilation, airconditioning, and refrigeration systems – an overview: Part II: Soft and
fusion control,” HVAC&R Research, vol. 17, no. 2, pp. 144–158, 2011.
[5] F. Oldewurtel, A. Parisio, C. N. Jones, D. Gyalistras, M. Gwerder,
V. Stauch, B. Lehmann, and M. Morari, “Use of model predictive control
and weather forecasts for energy efficient building climate control,”
Energy and Buildings, vol. 45, pp. 15–27, 2012.
[6] 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, vol. 100, no. 1,
pp. 240–253, 2012.
[7] Y. Ma, A. Kelman, A. Daly, and F. Borrelli, “Predictive control for
energy efficient buildings with thermal storage: modeling, stimulation,

[26]
[27]
[28]
[29]
[30]
[31]

[32]

and experiments,” IEEE Control System Magazine, vol. 32, no. 1, pp.
44–64, 2012.
H. Farhangi, “The path of the smart grid,” IEEE Power Energy Magazine, vol. 8, no. 1, pp. 18–28, 2010.
X. Zhang and A. Papachristodoulou, “A real-time control framework for
smart power networks: Design methodology and stability,” Automatica,
vol. 58, pp. 43–50, 2015.
X. Zhang, R. Kang, M. McCulloch, and A. Papachristodoulou, “Realtime active and reactive power regulation in power systems with tapchanging transformers and controllable loads,” Sustainable Energy,
Grids and Networks, vol. 5, pp. 27–38, 2016.
W. Shi, N. Li, X. Xie, C.-C. Chu, and R. Gadh, “Optimal residential
demand response in distribution network,” IEEE Journal on Selected
Areas in Communications, vol. 32, no. 7, pp. 1441–1450, 2014.
W. Shi, N. Li, C. C. Chu, and R. Gadh, “Real-time energy management
in microgrids,” IEEE Transactions on Smart Grid, vol. 8, no. 1, pp.
228–238, 2017.
A. Jadbabaie, J. Lin, and A. S. Morse, “Coordination of groups of mobile
autonomous agents using nearest neighbor rules,” IEEE Transactions on
Automatic Control, vol. 48, no. 6, pp. 988–1001, 2003.
F.-Y. Wang, “Parallel control and management for intelligent transportation systems: Concepts, architectures, and applications,” IEEE
Transactions on Intelligent Transportation Systems, vol. 11, no. 3, pp.
630–638, 2010.
P. D. Moroşan, R. Bourdais, D. Dumur, and J. Buisson, “Building
temperature regulation using a distributed model predictive control,”
Energy and Buildings, vol. 42, no. 9, pp. 1445–1452, 2010.
Y. Ma, G. Anderson, and F. Borrelli, “A distributed predictive control
approach to building temperature regulation,” in Proc. of 2011 American
Control Conference, 2011, pp. 2089–2094.
H. Hao, J. Lian, K. Kalsi, and J. Stoustrup, “Distributed flexibility characterization and resource allocation for multi-zone commercial buildings
in the smart grid,” in Proc. of 54th IEEE Conference on Decision and
Control, 2015, pp. 3161–3168.
K. Deng, P. Barooah, and P. G. Mehta, “Mean-field control for energy
efficient buildings,” in Proc. of 2012 American Control Conference,
2012, pp. 3044–3049.
S. C. Sugarman, HVAC Fundamentals. CRC Press, 2005.
G. Wang, Z. Wang, K. Xu, and M. Liu, “Air handling unit supply
air temperature optimization during economizer cycles,” in Proc. of
International Conference for Enhanced Building Operations, 2007.
Y. Lin, T. Middelkoop, and P. Barooah, “Issues in identification of
control-oriented thermal models of zones in multi-zone buildings,” in
Proc. of 51st IEEE Conference on Decision and Control, 2012, pp.
6932–6937.
M. J. Jiménez, H. Madsen, and K. K. Andersen, “Identification of the
main thermal characteristics of building components using MATLAB,”
Building and Environment, vol. 43, no. 2, pp. 170–180, 2008.
P. Bacher and H. Madsen, “Identifying suitable models for the heat
dynamics of buildings,” Building and Environment, vol. 47, no. 7, pp.
1511–1522, 2011.
H. Hao, Y. Lin, A. S. Kowli, P. Barooah, and S. Meyn, “Ancillary
service to the grid through control of fans in commercial building HVAC
systems,” IEEE Transactions on Smart Grid, vol. 5, no. 4, pp. 2066–
2074, 2014.
Y. Lin, P. Barooah, S. Meyn, and T. Middelkoop, “Experimental evaluation of frequency regulation from commercial building HVAC system,”
IEEE Transactions on Smart Grid, vol. 6, no. 2, pp. 776–783, 2015.
S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge
University Press, 2004.
D. Feijer and F. Paganini, “Stability of primal-dual gradient dynamics
and applications to network optimization,” Automatica, vol. 46, no. 12,
pp. 1974–1981, 2010.
A. Cherukuri, E. Mallada, and J. Cortés, “Asymptotic convergence of
constrained primal-dual dynamics,” System and Control Letters, vol. 87,
pp. 10–15, 2016.
X. Zhang, A. Papachristodoulou, and N. Li, “Distributed optimal steadystate control using reverse- and forward-engineering,” in Proc. of 54th
IEEE Conference on Decision and Control, 2015, pp. 5257–5264.
K. H. Ang, G. Chong, and Y. Li, “PID control system analysis, design,
and technology,” IEEE Transactions on Control Systems Technology,
vol. 13, no. 4, pp. 559–576, 2005.
M. Deru, K. Field, D. Studer, K. Benne, B. Griffith, P. Torcellini, B. Liu,
M. Halverson, D. Winiarski, M. Rosenberg, M. Yazdanian, J. Huang,
and D. Crawley, “Department of energy commercial reference building
models of the national building stock,” in Technical Report NREL/TP5500-46861, 2011.
X. Li, J. Wen, and E. W. Bai, “Developing a whole building cooling
energy forecasting model for on-line operation optimization using proactive system identification,” Applied Energy, vol. 164, pp. 69–88, 2016.

