# A Comfort-Aware Energy Efficient HVAC System Based on the Subspace Identification Method

Authors: O. Tsakiridis; D. Sklavounos; E. Zervas; J. Stonham
Journal: Journal of Energy
 (2016)

---

## Abstract

A proactive heating method is presented aiming at reducing the energy consumption in a HVAC system while maintaining the thermal comfort of the occupants. The proposed technique fuses time predictions for the zones’ temperatures, based on a deterministic subspace identification method, and zones’ occupancy predictions, based on a mobility model, in a decision scheme that is capable of regulating the balance between the total energy consumed and the total discomfort cost. Simulation results for various occupation-mobility models demonstrate the efficiency of the proposed technique.


---

## Body

## 1. Introduction

As the control and limitation of the energy consumption remain a field with an exceptional technological and economical interest, areas that have been considered as high energy consumers comprise very challenging research issues. According to the US Energy Information Administration from 2013 through 2040 the electricity consumption in the commercial and residential sectors will be increasing by 0.5% and 0.8% per year [1]. It is well established and widely accepted through research studies that the main energy consumers in the commercial and residential buildings are the Heat Ventilation and Air Conditioning (HVAC) systems, as well as the lighting systems.Buildings contributed a 41% (or 40 quadrillion btu) to the total US energy consumption in 2014. On an average, about 43% of the energy consumption in a commercial and residential building is due to HVAC systems [2]. Therefore, due to the high rate of energy consumption, the necessity of the demand-driven control in the HVAC systems has become inevitable. In modern buildings several sophisticated systems have been applied aiming at providing this type of control in the HVAC systems. The state-of-the-art technology of the HVAC control considers the occupancy of the zones a very important parameter, playing a key role in methods aiming at reducing energy consumption. The detection and prediction of the zones’ occupancy are a very challenging research field and several techniques have been proposed based on historic statistical data as well as on probabilistic models.Another equally important factor taken into account in advanced HVAC control systems is the thermal comfort of the occupants. The objective of modern HVAC control systems is to reduce energy consumption without compromising the comfort of the occupants. Wireless sensor networks, equipped with temperature, humidity, and occupancy detection sensor nodes, are nowadays the basic platform to build automated HVAC control systems. A number of methods aiming to maintain the thermal comfort while saving energy have been proposed. Some of them are described in Section2. Towards this direction, a novel technique is proposed in this paper, which aims at balancing the comfort and energy costs in a multizone system. The decisions on heating the zones or not may be taken either centrally or in a distributed manner by wireless sensor nodes scattered in the multizone system. In any case temperature and zone occupancy information must be exchanged between a node, responsible for a zone, and its neighboring nodes. The decision process itself relies on two kinds of predictions: (a) temperature-time predictions for the zones and (b) the zones’ occupancy profile. The emphasis in this paper is on the zones’ temperature predictions and to this end a deterministic subspace identification method is used for modelling the thermal dynamics of each zone. That is, each zone is modelled by a simple state-space model capable of producing accurate predictions based on the surrounding temperatures, the heating power of the zone, and the current state that summarizes the temperature history of the zone. For the zones’ occupancy predictions we consider a semi-Markov model, where occupants (moving as a swarm) stay in a zone for a random period of time and then move to adjacent zones with given probabilities. What is needed by the decision process is the distribution of the first entrance time to unoccupied zones. Aiming at a proactive action, the proposed method periodically computes the risk of activating the heater or not and decides in favor of the action that produces the smaller risk. The computation of the risks relies on the relative weights of the energy and discomfort costs so that the balance between the total energy consumed and the total discomfort cost may be regulated.The structure of the paper is as follows: In Section2, prior and related work of the field of the demand-driven HVAC systems that utilize occupancy and prediction methods is presented. In Section 3, the proposed comfort-aware energy efficient mechanism is described, along with the subspace identification method, used for obtaining temperature-time predictions, and a discussion on the distribution of the first entrance time to unoccupied zones. Section 4 contains simulation results and assesses the effectiveness of the proposed algorithm. Finally, Section 5 summarizes the paper and proposes research directions for future work.


## 2. Prior Work

The necessity of demand-driven HVAC systems, for energy efficient solutions, has orientated the researchers towards occupancy-based activated systems. In multizone spaces the reactive and proactive activation of the heating/cooling of the zones contributes to significant energy savings and improves the thermal comfort for the occupants. Several research works with valuable results of energy savings have utilized the occupancy detection and prediction, in order to control the HVAC system appropriately. For the occupancy detection several types of sensors are used (CO2, motion, etc.), while, for the prediction, combinational systems are usually applied, utilizing mathematical predictive models (e.g., Markov chains) alongside with detected real data.The authors in [3] proposed an automatic thermostat control system which is based on an occupancy prediction scheme that predicts the destination and arrival time of the occupants in the air-conditioned areas, in order to provide a comfortable environment. For the occupants’ mobility prediction the cell tower information of the mobile telephony system is utilized and the arrival time prediction is based on historical patterns and route classification. For the destination prediction of locations very close to each other (intracell), the time-aided order-Markov predictor is used. The authors in [4] proposed a closed-loop system for optimally controlling HVAC systems in buildings, based on actual occupancy levels named the Power-Efficient Occupancy-Based Energy Management (POEM) System. In order to accurately detect occupants’ transition, they deployed a wireless network comprising two parts: the occupancy estimation system (OPTNet) consisting of 22 camera nodes and a passive infrared (PIR) sensors system (BONet). By fusing the sensing data from the WSN (OPTNet and BONet) with the output of an occupancy transition model in a particle filter, more accurate estimation of the current occupancy in each room is achieved. Then, according to the current occupancy in each room and the predicted one from the transition model, a control schedule of the HVAC system takes over the preheat of the areas to the target temperature. In [5] the authors used real world data gathered from a wireless network of 16 smart cameras called Smart Camera Occupancy Position Estimation System (SCOPES) and in this way, they developed occupancy models. Three types of Markov chain (MC) occupancy models were tested and these are the single MC, the closest distance, and finally the blended (BMC). The authors concluded that BMC is the most efficient and thus they embodied it as the occupancy prediction method in their proposed “OBSERVE” algorithm, which is a temperature control strategy for HVAC systems.The authors in [6] have developed an integrated system called SENTINEL that is a control system for HVAC systems utilizing occupancy information. For the occupancy detection and localisation, the system utilizes the existing Wi-Fi network and the clients’ smart phones. The occupancy localization mechanism is based on the access point (AP) communication with a client’s smart phone, so that if an occupant’s phone sends packets to an AP then he is located within the range of the AP. By classifying the building areas into two main categories, namely, the personal and the shared spaces, the proposed mechanism activates the HVAC system when an owner of a personal space (e.g., office) has been detected within an area where his/her office is located or when occupants are detected within shared spaces. In [7] an integrated heating and cooling control system of a building is presented aiming to reduce energy consumption. The occupancy behavior prediction as well as the weather forecast, as inputs to a virtual (software based) building model, determines the control of the HVAC system. The occupancy detection technique utilizes Gaussian Mixture Models (GMM) for the categorization of selected features, yielding the highest information gain according to the different number of occupants. This categorization was used for observation to a Hidden Markov Model for the estimation of the number of occupants. A semi-Markov model was developed based on patterns comprised by sensory data of CO2, acoustics, motion, and light changes, to estimate the duration of occupants in the space. The work in [8] proposes a model predictive control (MPC) technique aiming to reduce energy consumption in a HVAC system while maintaining comfortable environment for the occupants. The occupancy predictive model is based on the two-state Markov chain, with the states modelling the occupied and occupied condition of the areas. The authors in [9] propose a feedback control algorithm for a variable air volume (VAV) HVAC system for full actuated (zones consisting of one room) and underactuated (zones consisting of more than one room) zones. The proposed algorithm is called MOBSua (Measured Occupancy-Based Setback for underactuated zones) and it utilizes real-time occupancy data, through a WSN, for optimum energy efficiency and thermal comfort of the occupants. Moreover, the algorithm can be applied on conventional control systems with no need of occupancy information and it is scalable to arbitrary sized buildings.


## 3. Comfort-Aware Energy Efficient Mechanism

We consider a multizone HVAC system consisting ofZ zones z
i, i
=
1,2
,
…
,
Z. Each zone is equipped with a wireless sensor node capable of conveying the zone’s temperature and occupancy information to its neighbors or to a central processing unit. The occupancy information may be very simple, that is, binary information indicating the presence or absence of individuals in the zone, or more advanced like the number of occupants in the zone. The objective is to minimize the total energy and discomfort cost defined as(1)
E
=
E
e
n
e
r
g
y
+
E
d
i
s
c
o
m
f
o
r
t
=
∑
i
=
1
Z
∫
-
∞
∞
W
i
I
H
i
=
on
⁡
d
t
+
∑
i
=
1
Z
∫
-
∞
∞
C
i
I
T
i
<
T
c
d
t
.Note that only the heating problem is addressed in this paper. Thus, the case of spending energy for cooling the zones and keeping the temperature T
i within a certain comfort zone is not treated for simplicity. In (1) I
(
A
) denotes the indicator function which takes the value 1 or 0 depending on whether condition A is true or false. W
i is the power of a heater that covers zone z
i and H
i is the state of the heater; that is, the heater is on or off. Thus, the first term of (1) is the total energy consumed by the multizone system. For the discomfort cost we define the cost per unit time C
i and the target comfort threshold T
c. As long as the zone’s temperature T
i is higher than T
c the occupants do not feel discomfort. The parameter C
i may depend on the zone, the number of the occupants (e.g., the total system’s discomfort cost is proportional to the number of people experiencing discomfort), and the difference of the zone’s temperature T
i and the comfort threshold T
c. It is worth noting that there is a tradeoff between energy savings and discomfort cost. We may preheat all zones to the comfort level thus rendering the discomfort cost equal to zero. This policy is inefficient due to the large amount of energy consumed. In the other extreme we could act reactively by heating a zone only upon detecting occupants in it. In this case the energy cost would be as small as possible but the discomfort cost could increase dramatically. In the sequel we will describe a method that acts proactively by taking periodically decisions on whether to heat a zone or not. If the entrance to a zone is delayed then we postpone heating until the next decision epoch. If, on the contrary, it is anticipated that a zone will be occupied before the zone’s temperature reaches the comfort level, then we preheat the zone.The decision of whether or not the heater of an unoccupied zonez
i should be turned on depends on (a) the current state of the zone, (b) the relative value of the energy cost (W
i) and discomfort cost (C
i) per unit time, (c) the temperatures of the surrounding zones, and (d) an estimate of the time that the zone will become occupied. This decision may be taken either centrally from the central processing unit that gathers information by all zones or in a distributed manner by each zone’s node after collecting the relative information. There is also the possibility of dividing the decision process functionality between the central unit and the wireless sensor network nodes. For example, the prediction of the first entrance times to unoccupied zones may be taken by the central unit that is aware of the location of people in the multizone system, and then the predicted values may be communicated back to the zones’ nodes for the final decision. The specific implementation of the decision process is irrelevant to this paper and it will not be further analyzed.We assume that time is discretizedt
=
m
T
s, T
s is the sampling period, and that each node takes its decisions periodically every D sampling periods. Note that there is no need for the nodes to be synchronized to each other. Let Y
i
(
t
) denote the random variable that models the remaining time from the current time instant t until entrance to the unoccupied zone z
i. It is obvious that the random variables Y
k
(
t
), for the various unoccupied zones z
k, do not follow the same distribution. We further assume that the node is capable of making predictions for the time it takes to exceed the comfort threshold T
c. To this end, let S
i
o
n
(
t
) denote the time period that is required to exceed T
c if the heater is turned on immediately. Similarly, we define S
i
o
f
f
(
t
) as the predicted time to reach the comfort threshold if the heater remains off for a period D and then, at the next decision epoch, it is turned on. The relation of the aforementioned parameters is shown in Figure 1.Figure 1
Predicted time intervals related to the decision process.Suppose now that at timet the node of the unoccupied zone z
i has to take a decision whether to turn the zone’s heater on or remain in the off state. At the decision epoch t the risk to turn the heater on is(2)
R
i
o
n
t
=
W
i
DP
S
i
o
f
f
t
≤
Y
i
t
+
W
i
Y
i
t
-
S
i
o
n
t
P
S
i
o
n
t
<
Y
i
t
≤
S
i
o
f
f
t
.That is, if S
i
o
f
f
(
t
)
≤
Y
i
(
t
), there is plenty of time to reach the comfort threshold even if we start heating at the next decision epoch and therefore we consume unnecessarily W
i
D units of energy (this is the case depicted in Figure 1). Similarly, if S
i
o
n
(
t
)
<
Y
i
(
t
)
≤
S
i
o
f
f
(
t
) we waste W
i
(
Y
i
(
t
)
-
S
i
o
n
(
t
)
) units of energy. Note that there is no risk of turning on the heater if the entrance to the zone happens sooner than it is predicted. In a similar fashion we calculate the risk of keeping the heater off. In this case(3)
R
i
o
f
f
t
=
C
i
DP
Y
i
t
≤
S
i
o
n
t
+
C
i
S
i
o
f
f
t
-
Y
i
t
P
S
i
o
n
t
<
Y
i
t
≤
S
i
o
f
f
t
.Using (2) and (3) the decision criterion takes the form(4)
R
i
o
n
t
≶
o
f
f
o
n
R
i
o
f
f
t
.Consider for the moment that there is no randomness on Y
i
(
t
), that is, strict time scheduled occupation of zones. Typical examples are the classrooms in a university campus. In this case, criterion (4) takes a very simple form since one of the involved probabilities is equal to one and the rest assume the value of zero. Thus, if Y
i
(
t
)
≤
S
i
o
n
(
t
), then the heater is turned on, if S
i
o
f
f
(
t
)
<
Y
i
(
t
), the heater remains off, and if S
i
o
n
(
t
)
<
Y
i
(
t
)
≤
S
i
o
f
f
(
t
) the heater will be turned on/off depending on the values of the relative energy and discomfort cost. The balance between these two costs is regulated by the weights W
i and C
i.IfY
i
(
t
) is random, then criterion (4) needs to be slightly modified as the value of Y
i
(
t
) is unknown. In this case we replace Y
i
(
t
) with its conditional expected value given that S
i
o
n
(
t
)
<
Y
i
(
t
)
≤
S
i
o
f
f
(
t
). Therefore, in criterion (4), we substitute Y
i
(
t
) by E
[
Y
i
(
t
)
∣
S
i
o
n
(
t
)
<
Y
i
(
t
)
≤
S
i
o
f
f
(
t
)
].From the previous discussion, it becomes clear that the decision process needs the estimatesS
i
o
n
(
t
), S
i
o
f
f
(
t
) and the distribution of the process Y
i
(
t
). The former is the subject of the next subsection, whereas a discussion on the distribution of Y
i
(
t
) is left for Section 3.2.


### 3.1. Temperature-Time Predictions Based on SID

It is conceivable that in order to make predictions about a zone’s temperature evolution we need a model able to capture the thermal dynamics of the tagged zone. To this end we resort to the deterministic subspace identification method. Each zone is treated separately and is represented by a discrete time, linear, time invariant, state-space model. That is,(5)
x
k
+
1
=
A
x
k
+
B
u
k
+
w
k
,
(6)
y
k
=
C
x
k
+
D
u
k
+
v
k
,where y
k is the zone temperature at time k and is in general an l-dimensional vector if l temperature sensors scattered in the zone have been deployed. Assuming a uniform zone temperature we use only one sensor node and thus y
k is scalar. v
k represents the measurement noise due to imperfections of the sensor. u
k is a vector of input measurements and its dimensionality varies from zone to zone. This vector contains the measurable output of sources that influence zone’s temperature. Such sources are the power of a heater located in the zone, the air temperature at various positions outside the zone, and so forth. The vector x
k is the state vector of the process at the discrete time k. The states have a conceptual relevance only and they are not assigned physical interpretation. Of course a similarity transform can convert the states to physical meaningful ones. The unmeasurable vector signal w
k represents noise due to the presence of humans in the zone, lamps’ switching on and off, and so forth. For the rest of the paper we neglect the effects of state and output noise and thus we deal with a deterministic identification problem, which is the computation of the matrices A, B, C, and D from the given input-output data. The no noise assumption makes sense since we need temperature-time predictions for the empty zones only which are “free” from occupants and other disturbances.In order to estimate the state-space matricesA, B, C, and D of each zone, we need a training phase during which all the nodes in the structure report their relative measurements to a central computing system. These measurements consist of the temperature of the zone and the current power of a heater, if the sensor node is in charge of a heated zone, or only the temperature if the sensor node monitors an area immaterial to the decision process but relevant to other zones, that is, the exterior of a building. The training phase should be performed only once but under controlled conditions, for example, no extra sources of heating and closed windows. Moreover, the length of this period should be quite large in order to exhibit several variations in the input signals to excite the modes of the system. In the simulation part we considered a training period of 24 h, that is, 86400 samples with sampling period 1 sec. The sampling period of 1 sec is too small for all practical reasons. However, as the dominant factor of performance degradation is the estimation of the residual arrival times we consider an ideal temperature prediction model (high accuracy sensors and temperature noiseless zones) and we access only the uncertainty of the occupation model.After receiving all the measurements, the central system organizes them to input-output data for each zone. For example, if a zone is namedA, neighbors zones are named B, C, and D and the exterior to the building zone is named E, then the temperature of zone A is the output of the system to be identified, whereas the temperatures of zones B, C, D, and E as well as the power profile of the heater of zone A are the input data for the system. A deterministic subspace identification process, described in the following paragraphs, is then run for each zone. The relevant matrices A, B, C, and D of the state-space model of each zone are communicated back to the sensor node responsible for the monitoring of the zone.Following the notation and the derivation in [10] we define the block Hankel matrices(7)
U
p
=
U
0
∣
i
-
1
=
u
0
u
1
⋯
u
j
-
1
u
1
u
2
⋯
u
j
⋮
⋮
⋮
⋮
u
i
-
1
u
i
⋯
u
i
+
j
-
2
,
U
f
=
U
i
∣
2
i
-
1
=
u
i
u
i
+
1
⋯
u
i
+
j
-
1
u
i
+
1
u
i
+
2
⋯
u
i
+
j
⋮
⋮
⋮
⋮
u
2
i
-
1
u
2
i
⋯
u
2
i
+
j
-
2which represent the “past” and the “future” input with respect to the present time instant i. Stacking U
p on top of U
f we have the block Hankel matrix(8)
U
0
∣
2
i
-
1
=
U
0
∣
i
-
1
U
i
∣
2
i
-
1
=
U
p
U
fwhich can also be partitioned as(9)
U
0
∣
2
i
-
1
=
U
0
∣
i
U
i
+
1
∣
2
i
-
1
=
U
p
+
U
f
-by moving the “present” time i one step ahead to time i
+
1. Similarly, we define the output block Hankel matrices Y
0
∣
2
i
-
1, Y
p, Y
f, Y
p
+, Y
f
- and the input-output data matrices(10)
W
p
=
W
0
∣
i
-
1
=
U
0
∣
i
-
1
Y
0
∣
i
-
1
=
U
p
Y
p
,
W
p
+
=
U
p
+
Y
p
+
.The state sequence X
m is defined as(11)
X
m
=
x
m
x
m
+
1
⋯
x
m
+
j
-
1and the “past” and “future” state sequences are X
p
=
X
0 and X
f
=
X
i, respectively. With the state-space model (13), (15) we associate the extended observability matrix Γ
i and the reversed extended controllability matrix Δ
i, where(12)
Γ
i
=
C
C
A
C
A
2
⋮
C
A
i
-
1
,
Δ
i
=
A
i
-
1
B
A
i
-
2
B
⋯
A
B
B
.The subspace identification method relies on determining the state sequence X
f and the extended observability matrix Γ
i directly from the input-output data u
k, y
k and based on them, in a next step, to extract the matrices A, B, C, and D. To this end we define the oblique projection of the row space of Y
f along the row space of U
f on the row space of W
p:(13)
O
i
=
Y
f
/
U
f
W
p
=
Y
f
/
U
f
⊥
,
W
p
/
U
f
⊥
†
W
p
,where † denotes the Moore-Penrose pseudoinverse of a matrix and A
/
B
⊥ is a shorthand for the projection of the row space of matrix A onto the orthogonal complement of the row space of matrix B; that is,(14)
A
/
B
⊥
=
I
-
A
B
T
B
B
T
†
B
.As it is proved in [10] the matrix O
i is equal to the product of the extended observability matrix and the “future” state vector(15)
O
i
=
Γ
i
X
f
.Having in our disposal the matrix O
i we use its singular value decomposition(16)
O
i
=
U
1
U
2
S
1
0
0
0
V
1
T
V
2
T
=
U
1
S
1
V
1
Tand we identify the extended observation matrix Γ
i and the state vector X
f as(17)
Γ
i
=
U
1
S
1
1
/
2
T
,
X
f
=
T
-
1
S
1
1
/
2
V
1
T
=
Γ
i
†
O
i
,where T is an n
×
n arbitrary nonsingular matrix representing a similarity transformation.One method [10] to extract the matrices A, B, C, and D uses in addition the oblique projection(18)
O
i
-
1
=
Y
f
-
/
U
f
-
W
p
†
=
Γ
i
-
1
X
i
+
1
,where Γ
i
-
1 is obtained from Γ
i by deleting the last l rows (l
=
1 in our case). Similar to the previous derivation we have(19)
X
i
+
1
=
Γ
i
-
1
†
O
i
-
1and the matrices A, B, C, and D are obtained by solving (in least squares fashion) the system(20)
X
i
+
1
Y
i
∣
i
=
A
B
C
D
X
i
U
i
∣
i
.After completion of the training phase the state-space model matricesA, B, C, and D for each zone are communicated back to the node responsible for the zone. Upon reception of the matrices the nodes enter a “convergence” phase; that is, starting from an initial state x
0, they update the state (5) for a predetermined period of time, using the measurements of the surrounding zones and the zone’s heater status. Note that the training and the convergence phase need to be executed only once prior to the normal operation of the nodes. After the convergence phase the nodes enter the “decision” phase. During this phase the nodes update the state of the model (as in the convergence phase) and periodically they take decisions on whether to turn the zone’s heater on or not. Suppose that at the decision epoch t the state of the tagged unoccupied node is x
t. The node bases its decision on the estimates S
i
o
n
(
t
) and S
i
o
f
f
(
t
). Recall that S
i
o
n
(
t
) is the time needed to reach the comfort level T
c if the heater is on. That is,(21)
S
i
o
n
t
=
min
⁡
l
:
y
t
+
l
≥
T
c
.The node “freezes” the input-vector u
k to its current value u
t (more recent temperatures of the surrounding zones and the tagged zone’s heater is on) and starting from state x
t repeats (5), (6) until y
t
+
k exceeds the value T
c. Note that(22)
y
t
+
k
=
C
A
k
x
t
+
C
A
k
-
1
B
u
t
+
C
A
k
-
2
B
u
t
+
1
+
⋯
+
C
A
B
u
t
+
k
-
2
+
C
B
u
t
+
k
-
1
+
D
u
t
+
kand for u
t
=
u
t
+
1
=
⋯
=
u
t
+
k we obtain(23)
y
t
+
k
=
C
A
k
x
t
+
C
I
+
A
+
⋯
+
A
k
-
1
B
u
t
+
D
u
t
=
C
A
k
x
t
+
C
I
-
A
-
1
I
-
A
k
B
u
t
+
D
u
t
.If we set y
t
+
k
=
T
c and(24)
w
=
T
c
-
C
I
-
A
-
1
B
u
t
-
D
u
t
,then (23) takes the form(25)
w
=
C
A
k
x
t
-
C
I
-
A
-
1
A
k
B
u
t
=
C
V
Λ
k
V
-
1
x
t
-
C
I
-
A
-
1
V
Λ
k
V
-
1
B
u
t
,where we have considered the diagonalization of the matrix A as A
=
V
Λ
V
-
1, with Λ being the diagonal matrix of the eigenvalues and V being the matrix of the corresponding eigenvectors. Numerical methods can be used to solve (25) for k. In a similar manner we compute S
i
o
f
f
(
t
). In this case, however, we first repeat (5) D times with input-vector u
k reflecting the fact that the heater is off, and then starting from the new state x
t
+
D we estimate the time needed to reach T
c if the heater is on.


### 3.2. On the Distribution ofY
n
(
t
)

The random variableY
n
(
t
) expresses the remaining time from the current time instant t until entrance to zone z
n. That is, we assume that occupants move around as a swarm, visiting zones either in a predetermined order or in a random fashion and finally they end up in zone z
n. Each time the occupants visit a zone z
l they stay in it for a nonzero time period which is itself a random variable with distribution F
l
(
x
). We use the variable t to designate calendar time, the time elapsed since the process started, and τ to designate the time elapsed since entry to a zone. The zone occupation time distribution F
l
(
x
) may be given parametrically either a priori or inferred from examination of historic data by suitable fitting of the distribution’s parameters. We may use the empirical distribution function (e.d.f.) instead which is the nonparametric estimate of F
l
(
x
). That is, if a sample of N occupation periods of zone z
l, {
x
1
,
x
2
,
…
,
x
N
} is available, then(26)
F
^
l
x
=
1
N
∑
m
=
1
N
I
x
m
∈
0
,
x
,where I
(
·
) is the indicator function.Let us consider first the simple case of two adjusting zonesz
l, z
n with zone z
l being occupied and z
n unoccupied as it is shown in Figure 2. We assume that the occupants remain to zone z
l for a random period of time X
l which is distributed according to F
l
(
x
) and then upon departure they enter zone z
n. In this case Y
n
(
t
) expresses the remaining waiting time to zone z
l and follows the residual life distribution which is defined by(27)
R
¯
Y
n
t
y
=
1
-
R
Y
n
t
y
=
P
X
l
>
τ
+
y
∣
X
l
>
τ
=
F
¯
l
τ
+
y
F
¯
l
τ
.Figure 2
Movement of occupants from zonez
l to the adjacent unoccupied zone z
n.Next we consider the case that occupants are currently in a zone (call itz
0) and upon departure they visit in cascade zones z
1
,
z
2
,
…
,
z
n as it is shown in Figure 3. Then(28)
Y
n
t
=
Y
0
+
X
1
+
X
2
+
⋯
+
X
n
-
1
,where Y
0 is the remaining time in zone z
0 and X
i is the occupation period of zone z
i. It is well known that the distribution of Y
n
(
t
) is the convolution of the distributions of the random variables Y
0
,
X
1
,
…
,
X
n
-
1. That is(29)
F
Y
n
t
y
=
P
Y
n
t
≤
y
=
R
0
y
⋆
F
1
y
⋆
F
2
y
⋆
⋯
⋆
F
n
-
1
ywith R
0
(
y
)
=
1
-
F
¯
0
(
τ
+
y
)
/
F
¯
0
(
τ
). A more compact form is obtained by considering a transform of the distributions, that is, the Laplace transform, the moment generating function, or the characteristic function. Using the moment generating function (mgf) which for a distribution F
X
(
x
) of a nonnegative random variable X is defined by(30)
F
X
s
=
∫
0
∞
e
s
x
d
F
X
x
, (29) takes the form(31)
F
Y
n
t
s
=
R
0
s
·
F
1
s
·
F
2
s
·
…
·
F
n
-
1
s
.In some cases (unfortunately few) (26) can be inversed analytically in order to calculate the distribution F
Y
n
(
t
)
(
y
). For example, if the occupation periods X
i, i
=
0
,
…
,
n
-
1, are exponentially distributed with corresponding parameters λ
i, then Y
0 is also exponentially distributed with parameter λ
0 due to the memoryless property of the exponential distribution, and the probability density function of Y
n
(
t
) is (assuming that the parameters λ
i are all distinct)(32)
f
Y
n
t
x
=
∑
i
=
0
n
-
1
λ
0
⋯
λ
n
-
1
∏
j
=
0
,
j
≠
i
n
-
1
λ
j
-
λ
i
exp
⁡
-
x
λ
i
.(The expression for the general case, of nondistinct parameters, is slightly more involved but still exists.) If a closed form for the density f
Y
n
(
t
)
(
x
) does not exist or it is difficult to be obtained, we resort to the saddlepoint approximation [11] to invert F
Y
n
(
t
)
(
s
). To this end, for a random variable X with cumulative distribution function (cdf) F
X
(
x
), we define the cumulant generating function (cgf) K
X
(
s
) as the logarithm of the corresponding mfg F
X
(
s
); that is,(33)
K
X
s
=
log
⁡
F
X
s
.Then, the saddlepoint approximation of the probability density function (pdf) f
X
(
x
) is(34)
f
X
x
≈
1
2
π
K
x
′
′
s
^
x
1
/
2
exp
⁡
K
X
s
^
x
-
s
^
x
x
,where s
^
(
x
) is the solution to the saddlepoint equation(35)
K
X
′
s
=
x
.Estimation (34) relies on numerical computations based on cgf K
X
(
S
), which in turn depends on F
X
(
s
). Note that we may use the empirical mgf F
^
X
(
s
) instead, so that a nonparametric estimation of the pdf f
X
(
x
) is possible.Figure 3
Movement of occupants in cascade from zonez
0 to zone z
n.The more general case is to model the movement of the occupants in the multizone system as a semi-Markov process. The states of the process are the zones of the system and at the transition times they form an embedded Markov chain with transition probabilitiesp
i
j. Given a transition from zone z
i to zone z
j, the occupation time in z
i has distribution function F
i
j
(
x
). The matrix Q
(
x
) with elements Q
i
j
(
x
)
=
p
i
j
F
i
j
(
x
) is the so called semi-Markov kernel. The occupation time distribution in zone z
i independent of the state to which a transition is made is F
i
(
x
)
=
∑
j
Q
i
j
(
x
). Consider now the case that at time t occupants exist in zone z
i and the node of zone z
n has to make a decision on whether to start heating the zone or not. For the decision process the distribution of time until entrance to zone z
n, Y
n
(
t
), is needed. Y
n
(
t
) is the sum of two terms. The first term is the residual occupation time at zone z
i and the second term is the time to reach z
n upon departure from zone z
i. The distribution R
i
(
x
) of the first term is given by (27), whereas for the second term we have to consider the various paths taken upon departure from state z
i. Thus, if zone z
i neighbors N
i zones, namely, z
k, k
=
1
,
…
,
N
i, then the time to reach zone z
n upon departure is distributed according to(36)
G
i
n
x
=
∑
k
=
1
N
i
p
i
k
H
k
n
x
,where H
k
n
(
x
) is the distribution of the first passage from zone z
k to zone z
n. Pyke [12] and Mason [13, 14] provided solutions for the first passage distribution in a semi-Markov process. For example, Pyke’s formula states that H
(
s
), the elementwise transform of the matrix H
(
x
)
=
{
H
i
j
(
x
)
}, is given by(37)
H
s
=
T
s
I
-
T
s
-
1
I
-
T
s
-
1
d
-
1
,where T
(
s
) is the transmittance matrix with the transforms T
i
j
(
s
) of the distributions F
i
j
(
x
). Notation M
d denotes the matrix that is formed by the diagonal elements of M. A simplified version of (37) is obtained if we let Δ
(
s
)
=
det
⁡
(
I
-
T
(
s
)
) and Δ
i
j
=
(
-
1
)
i
+
j
det
⁡
(
I
-
T
(
s
)
)
i
j, the i
jth cofactor of I
-
T
(
s
). Using this notation, for a semi-Markov process of n states, we get(38)
H
1
n
s
=
Δ
n
1
s
Δ
n
n
s
.By relabeling the states 1
,
…
,
n
-
1 of the process we can find the transform of the first passage from any zone to zone z
n.In summary,(39)
F
Y
n
t
s
=
R
i
s
·
G
i
n
s
=
R
i
s
·
∑
k
=
1
N
i
p
i
k
H
k
n
s
,with H
k
n
(
s
) computed using (38). Next, the saddlepoint approximation technique can be used to calculate the density f
Y
n
(
t
)
(
x
).We have to mention that the process of estimating the first passage distributions needs to be performed only prior to the decision process. Thus, given the topology of the multizone system, we posit a parametric family for the occupation time distributions and we estimate the distributions’ parameters from sample data. Then, the transition probabilities are estimated based on the data and the occupation mgfF
i
j
(
s
) is obtained. Having in our disposal the transforms F
i
j
(
s
) we can form the transmittance matrix T
(
s
) and solve for the first passage transform H
i
j
(
s
) using (38). Finally, numerical methods may be used to invert these transforms in order to obtain the necessary probability density functions.


## 3.1. Temperature-Time Predictions Based on SID

It is conceivable that in order to make predictions about a zone’s temperature evolution we need a model able to capture the thermal dynamics of the tagged zone. To this end we resort to the deterministic subspace identification method. Each zone is treated separately and is represented by a discrete time, linear, time invariant, state-space model. That is,(5)
x
k
+
1
=
A
x
k
+
B
u
k
+
w
k
,
(6)
y
k
=
C
x
k
+
D
u
k
+
v
k
,where y
k is the zone temperature at time k and is in general an l-dimensional vector if l temperature sensors scattered in the zone have been deployed. Assuming a uniform zone temperature we use only one sensor node and thus y
k is scalar. v
k represents the measurement noise due to imperfections of the sensor. u
k is a vector of input measurements and its dimensionality varies from zone to zone. This vector contains the measurable output of sources that influence zone’s temperature. Such sources are the power of a heater located in the zone, the air temperature at various positions outside the zone, and so forth. The vector x
k is the state vector of the process at the discrete time k. The states have a conceptual relevance only and they are not assigned physical interpretation. Of course a similarity transform can convert the states to physical meaningful ones. The unmeasurable vector signal w
k represents noise due to the presence of humans in the zone, lamps’ switching on and off, and so forth. For the rest of the paper we neglect the effects of state and output noise and thus we deal with a deterministic identification problem, which is the computation of the matrices A, B, C, and D from the given input-output data. The no noise assumption makes sense since we need temperature-time predictions for the empty zones only which are “free” from occupants and other disturbances.In order to estimate the state-space matricesA, B, C, and D of each zone, we need a training phase during which all the nodes in the structure report their relative measurements to a central computing system. These measurements consist of the temperature of the zone and the current power of a heater, if the sensor node is in charge of a heated zone, or only the temperature if the sensor node monitors an area immaterial to the decision process but relevant to other zones, that is, the exterior of a building. The training phase should be performed only once but under controlled conditions, for example, no extra sources of heating and closed windows. Moreover, the length of this period should be quite large in order to exhibit several variations in the input signals to excite the modes of the system. In the simulation part we considered a training period of 24 h, that is, 86400 samples with sampling period 1 sec. The sampling period of 1 sec is too small for all practical reasons. However, as the dominant factor of performance degradation is the estimation of the residual arrival times we consider an ideal temperature prediction model (high accuracy sensors and temperature noiseless zones) and we access only the uncertainty of the occupation model.After receiving all the measurements, the central system organizes them to input-output data for each zone. For example, if a zone is namedA, neighbors zones are named B, C, and D and the exterior to the building zone is named E, then the temperature of zone A is the output of the system to be identified, whereas the temperatures of zones B, C, D, and E as well as the power profile of the heater of zone A are the input data for the system. A deterministic subspace identification process, described in the following paragraphs, is then run for each zone. The relevant matrices A, B, C, and D of the state-space model of each zone are communicated back to the sensor node responsible for the monitoring of the zone.Following the notation and the derivation in [10] we define the block Hankel matrices(7)
U
p
=
U
0
∣
i
-
1
=
u
0
u
1
⋯
u
j
-
1
u
1
u
2
⋯
u
j
⋮
⋮
⋮
⋮
u
i
-
1
u
i
⋯
u
i
+
j
-
2
,
U
f
=
U
i
∣
2
i
-
1
=
u
i
u
i
+
1
⋯
u
i
+
j
-
1
u
i
+
1
u
i
+
2
⋯
u
i
+
j
⋮
⋮
⋮
⋮
u
2
i
-
1
u
2
i
⋯
u
2
i
+
j
-
2which represent the “past” and the “future” input with respect to the present time instant i. Stacking U
p on top of U
f we have the block Hankel matrix(8)
U
0
∣
2
i
-
1
=
U
0
∣
i
-
1
U
i
∣
2
i
-
1
=
U
p
U
fwhich can also be partitioned as(9)
U
0
∣
2
i
-
1
=
U
0
∣
i
U
i
+
1
∣
2
i
-
1
=
U
p
+
U
f
-by moving the “present” time i one step ahead to time i
+
1. Similarly, we define the output block Hankel matrices Y
0
∣
2
i
-
1, Y
p, Y
f, Y
p
+, Y
f
- and the input-output data matrices(10)
W
p
=
W
0
∣
i
-
1
=
U
0
∣
i
-
1
Y
0
∣
i
-
1
=
U
p
Y
p
,
W
p
+
=
U
p
+
Y
p
+
.The state sequence X
m is defined as(11)
X
m
=
x
m
x
m
+
1
⋯
x
m
+
j
-
1and the “past” and “future” state sequences are X
p
=
X
0 and X
f
=
X
i, respectively. With the state-space model (13), (15) we associate the extended observability matrix Γ
i and the reversed extended controllability matrix Δ
i, where(12)
Γ
i
=
C
C
A
C
A
2
⋮
C
A
i
-
1
,
Δ
i
=
A
i
-
1
B
A
i
-
2
B
⋯
A
B
B
.The subspace identification method relies on determining the state sequence X
f and the extended observability matrix Γ
i directly from the input-output data u
k, y
k and based on them, in a next step, to extract the matrices A, B, C, and D. To this end we define the oblique projection of the row space of Y
f along the row space of U
f on the row space of W
p:(13)
O
i
=
Y
f
/
U
f
W
p
=
Y
f
/
U
f
⊥
,
W
p
/
U
f
⊥
†
W
p
,where † denotes the Moore-Penrose pseudoinverse of a matrix and A
/
B
⊥ is a shorthand for the projection of the row space of matrix A onto the orthogonal complement of the row space of matrix B; that is,(14)
A
/
B
⊥
=
I
-
A
B
T
B
B
T
†
B
.As it is proved in [10] the matrix O
i is equal to the product of the extended observability matrix and the “future” state vector(15)
O
i
=
Γ
i
X
f
.Having in our disposal the matrix O
i we use its singular value decomposition(16)
O
i
=
U
1
U
2
S
1
0
0
0
V
1
T
V
2
T
=
U
1
S
1
V
1
Tand we identify the extended observation matrix Γ
i and the state vector X
f as(17)
Γ
i
=
U
1
S
1
1
/
2
T
,
X
f
=
T
-
1
S
1
1
/
2
V
1
T
=
Γ
i
†
O
i
,where T is an n
×
n arbitrary nonsingular matrix representing a similarity transformation.One method [10] to extract the matrices A, B, C, and D uses in addition the oblique projection(18)
O
i
-
1
=
Y
f
-
/
U
f
-
W
p
†
=
Γ
i
-
1
X
i
+
1
,where Γ
i
-
1 is obtained from Γ
i by deleting the last l rows (l
=
1 in our case). Similar to the previous derivation we have(19)
X
i
+
1
=
Γ
i
-
1
†
O
i
-
1and the matrices A, B, C, and D are obtained by solving (in least squares fashion) the system(20)
X
i
+
1
Y
i
∣
i
=
A
B
C
D
X
i
U
i
∣
i
.After completion of the training phase the state-space model matricesA, B, C, and D for each zone are communicated back to the node responsible for the zone. Upon reception of the matrices the nodes enter a “convergence” phase; that is, starting from an initial state x
0, they update the state (5) for a predetermined period of time, using the measurements of the surrounding zones and the zone’s heater status. Note that the training and the convergence phase need to be executed only once prior to the normal operation of the nodes. After the convergence phase the nodes enter the “decision” phase. During this phase the nodes update the state of the model (as in the convergence phase) and periodically they take decisions on whether to turn the zone’s heater on or not. Suppose that at the decision epoch t the state of the tagged unoccupied node is x
t. The node bases its decision on the estimates S
i
o
n
(
t
) and S
i
o
f
f
(
t
). Recall that S
i
o
n
(
t
) is the time needed to reach the comfort level T
c if the heater is on. That is,(21)
S
i
o
n
t
=
min
⁡
l
:
y
t
+
l
≥
T
c
.The node “freezes” the input-vector u
k to its current value u
t (more recent temperatures of the surrounding zones and the tagged zone’s heater is on) and starting from state x
t repeats (5), (6) until y
t
+
k exceeds the value T
c. Note that(22)
y
t
+
k
=
C
A
k
x
t
+
C
A
k
-
1
B
u
t
+
C
A
k
-
2
B
u
t
+
1
+
⋯
+
C
A
B
u
t
+
k
-
2
+
C
B
u
t
+
k
-
1
+
D
u
t
+
kand for u
t
=
u
t
+
1
=
⋯
=
u
t
+
k we obtain(23)
y
t
+
k
=
C
A
k
x
t
+
C
I
+
A
+
⋯
+
A
k
-
1
B
u
t
+
D
u
t
=
C
A
k
x
t
+
C
I
-
A
-
1
I
-
A
k
B
u
t
+
D
u
t
.If we set y
t
+
k
=
T
c and(24)
w
=
T
c
-
C
I
-
A
-
1
B
u
t
-
D
u
t
,then (23) takes the form(25)
w
=
C
A
k
x
t
-
C
I
-
A
-
1
A
k
B
u
t
=
C
V
Λ
k
V
-
1
x
t
-
C
I
-
A
-
1
V
Λ
k
V
-
1
B
u
t
,where we have considered the diagonalization of the matrix A as A
=
V
Λ
V
-
1, with Λ being the diagonal matrix of the eigenvalues and V being the matrix of the corresponding eigenvectors. Numerical methods can be used to solve (25) for k. In a similar manner we compute S
i
o
f
f
(
t
). In this case, however, we first repeat (5) D times with input-vector u
k reflecting the fact that the heater is off, and then starting from the new state x
t
+
D we estimate the time needed to reach T
c if the heater is on.


## 3.2. On the Distribution ofY
n
(
t
)

The random variableY
n
(
t
) expresses the remaining time from the current time instant t until entrance to zone z
n. That is, we assume that occupants move around as a swarm, visiting zones either in a predetermined order or in a random fashion and finally they end up in zone z
n. Each time the occupants visit a zone z
l they stay in it for a nonzero time period which is itself a random variable with distribution F
l
(
x
). We use the variable t to designate calendar time, the time elapsed since the process started, and τ to designate the time elapsed since entry to a zone. The zone occupation time distribution F
l
(
x
) may be given parametrically either a priori or inferred from examination of historic data by suitable fitting of the distribution’s parameters. We may use the empirical distribution function (e.d.f.) instead which is the nonparametric estimate of F
l
(
x
). That is, if a sample of N occupation periods of zone z
l, {
x
1
,
x
2
,
…
,
x
N
} is available, then(26)
F
^
l
x
=
1
N
∑
m
=
1
N
I
x
m
∈
0
,
x
,where I
(
·
) is the indicator function.Let us consider first the simple case of two adjusting zonesz
l, z
n with zone z
l being occupied and z
n unoccupied as it is shown in Figure 2. We assume that the occupants remain to zone z
l for a random period of time X
l which is distributed according to F
l
(
x
) and then upon departure they enter zone z
n. In this case Y
n
(
t
) expresses the remaining waiting time to zone z
l and follows the residual life distribution which is defined by(27)
R
¯
Y
n
t
y
=
1
-
R
Y
n
t
y
=
P
X
l
>
τ
+
y
∣
X
l
>
τ
=
F
¯
l
τ
+
y
F
¯
l
τ
.Figure 2
Movement of occupants from zonez
l to the adjacent unoccupied zone z
n.Next we consider the case that occupants are currently in a zone (call itz
0) and upon departure they visit in cascade zones z
1
,
z
2
,
…
,
z
n as it is shown in Figure 3. Then(28)
Y
n
t
=
Y
0
+
X
1
+
X
2
+
⋯
+
X
n
-
1
,where Y
0 is the remaining time in zone z
0 and X
i is the occupation period of zone z
i. It is well known that the distribution of Y
n
(
t
) is the convolution of the distributions of the random variables Y
0
,
X
1
,
…
,
X
n
-
1. That is(29)
F
Y
n
t
y
=
P
Y
n
t
≤
y
=
R
0
y
⋆
F
1
y
⋆
F
2
y
⋆
⋯
⋆
F
n
-
1
ywith R
0
(
y
)
=
1
-
F
¯
0
(
τ
+
y
)
/
F
¯
0
(
τ
). A more compact form is obtained by considering a transform of the distributions, that is, the Laplace transform, the moment generating function, or the characteristic function. Using the moment generating function (mgf) which for a distribution F
X
(
x
) of a nonnegative random variable X is defined by(30)
F
X
s
=
∫
0
∞
e
s
x
d
F
X
x
, (29) takes the form(31)
F
Y
n
t
s
=
R
0
s
·
F
1
s
·
F
2
s
·
…
·
F
n
-
1
s
.In some cases (unfortunately few) (26) can be inversed analytically in order to calculate the distribution F
Y
n
(
t
)
(
y
). For example, if the occupation periods X
i, i
=
0
,
…
,
n
-
1, are exponentially distributed with corresponding parameters λ
i, then Y
0 is also exponentially distributed with parameter λ
0 due to the memoryless property of the exponential distribution, and the probability density function of Y
n
(
t
) is (assuming that the parameters λ
i are all distinct)(32)
f
Y
n
t
x
=
∑
i
=
0
n
-
1
λ
0
⋯
λ
n
-
1
∏
j
=
0
,
j
≠
i
n
-
1
λ
j
-
λ
i
exp
⁡
-
x
λ
i
.(The expression for the general case, of nondistinct parameters, is slightly more involved but still exists.) If a closed form for the density f
Y
n
(
t
)
(
x
) does not exist or it is difficult to be obtained, we resort to the saddlepoint approximation [11] to invert F
Y
n
(
t
)
(
s
). To this end, for a random variable X with cumulative distribution function (cdf) F
X
(
x
), we define the cumulant generating function (cgf) K
X
(
s
) as the logarithm of the corresponding mfg F
X
(
s
); that is,(33)
K
X
s
=
log
⁡
F
X
s
.Then, the saddlepoint approximation of the probability density function (pdf) f
X
(
x
) is(34)
f
X
x
≈
1
2
π
K
x
′
′
s
^
x
1
/
2
exp
⁡
K
X
s
^
x
-
s
^
x
x
,where s
^
(
x
) is the solution to the saddlepoint equation(35)
K
X
′
s
=
x
.Estimation (34) relies on numerical computations based on cgf K
X
(
S
), which in turn depends on F
X
(
s
). Note that we may use the empirical mgf F
^
X
(
s
) instead, so that a nonparametric estimation of the pdf f
X
(
x
) is possible.Figure 3
Movement of occupants in cascade from zonez
0 to zone z
n.The more general case is to model the movement of the occupants in the multizone system as a semi-Markov process. The states of the process are the zones of the system and at the transition times they form an embedded Markov chain with transition probabilitiesp
i
j. Given a transition from zone z
i to zone z
j, the occupation time in z
i has distribution function F
i
j
(
x
). The matrix Q
(
x
) with elements Q
i
j
(
x
)
=
p
i
j
F
i
j
(
x
) is the so called semi-Markov kernel. The occupation time distribution in zone z
i independent of the state to which a transition is made is F
i
(
x
)
=
∑
j
Q
i
j
(
x
). Consider now the case that at time t occupants exist in zone z
i and the node of zone z
n has to make a decision on whether to start heating the zone or not. For the decision process the distribution of time until entrance to zone z
n, Y
n
(
t
), is needed. Y
n
(
t
) is the sum of two terms. The first term is the residual occupation time at zone z
i and the second term is the time to reach z
n upon departure from zone z
i. The distribution R
i
(
x
) of the first term is given by (27), whereas for the second term we have to consider the various paths taken upon departure from state z
i. Thus, if zone z
i neighbors N
i zones, namely, z
k, k
=
1
,
…
,
N
i, then the time to reach zone z
n upon departure is distributed according to(36)
G
i
n
x
=
∑
k
=
1
N
i
p
i
k
H
k
n
x
,where H
k
n
(
x
) is the distribution of the first passage from zone z
k to zone z
n. Pyke [12] and Mason [13, 14] provided solutions for the first passage distribution in a semi-Markov process. For example, Pyke’s formula states that H
(
s
), the elementwise transform of the matrix H
(
x
)
=
{
H
i
j
(
x
)
}, is given by(37)
H
s
=
T
s
I
-
T
s
-
1
I
-
T
s
-
1
d
-
1
,where T
(
s
) is the transmittance matrix with the transforms T
i
j
(
s
) of the distributions F
i
j
(
x
). Notation M
d denotes the matrix that is formed by the diagonal elements of M. A simplified version of (37) is obtained if we let Δ
(
s
)
=
det
⁡
(
I
-
T
(
s
)
) and Δ
i
j
=
(
-
1
)
i
+
j
det
⁡
(
I
-
T
(
s
)
)
i
j, the i
jth cofactor of I
-
T
(
s
). Using this notation, for a semi-Markov process of n states, we get(38)
H
1
n
s
=
Δ
n
1
s
Δ
n
n
s
.By relabeling the states 1
,
…
,
n
-
1 of the process we can find the transform of the first passage from any zone to zone z
n.In summary,(39)
F
Y
n
t
s
=
R
i
s
·
G
i
n
s
=
R
i
s
·
∑
k
=
1
N
i
p
i
k
H
k
n
s
,with H
k
n
(
s
) computed using (38). Next, the saddlepoint approximation technique can be used to calculate the density f
Y
n
(
t
)
(
x
).We have to mention that the process of estimating the first passage distributions needs to be performed only prior to the decision process. Thus, given the topology of the multizone system, we posit a parametric family for the occupation time distributions and we estimate the distributions’ parameters from sample data. Then, the transition probabilities are estimated based on the data and the occupation mgfF
i
j
(
s
) is obtained. Having in our disposal the transforms F
i
j
(
s
) we can form the transmittance matrix T
(
s
) and solve for the first passage transform H
i
j
(
s
) using (38). Finally, numerical methods may be used to invert these transforms in order to obtain the necessary probability density functions.


## 4. Simulation Results

We simulated a multizone system consisting of a squared arrangement of rooms (zones) where each room (zone) is equipped with a wireless sensor node as shown in Figure4. We assume that heat transfer takes place between room’s air and the walls as well as between room’s air and the roof. Furthermore, we neglect the effects of ground on the room temperature. North and south walls have the same effect on the room temperature and we assume the same for the east and west walls. According to these assumptions, there is a symmetry in the dynamics of the zones. For example, zones I, III, VII, and IX exhibit the same behavior. Note that, for the multizone system of Figure 4 several state variables, such as wall temperatures, are common to the individual zone systems. In the simulations we treat the multizone system as a single system with 42 states and nine outputs (the zones’ temperatures). The parameter values for this lumped capacity zone model were taken from [15] (see also [16]).Figure 4
Multizone model.In Figure4, T
o denotes the outside temperature which is assumed to be uniform with no loss of generality. For simulation purposes daily outside temperature variations are obtained using Walters’ model [17]. For this model, the maximum temperature of day T
m
a
x, the minimum temperature of day T
m
i
n, and the mean of the 24 hourly temperatures T
m need to be provided in order to estimate T
o. As an example, if T
m
a
x
=
15°C, T
m
i
n
=
2°C, and T
m
=
7°C, the daily variations of temperature in Figure 5 are obtained, which shows that temperature is higher at 2:00–3:00 pm.Figure 5
Empirical estimation ofT
o.


### 4.1. Training Phase

For the deterministic subsystem identification we used 86400 samples with sampling periodT
s
=
1 sec. For the outside temperature we used Walter’s model with parameters given in the first row of Table 1 and a heater gain equal to 600 W. We set the zones’ target temperatures T
t
g to those of Table 2 (first row) with a margin equal to T
m
g
=
1°C. That is, for a specific zone the heater is on until temperature T
t
g
+
T
m
g is reached. After this point the heater is turned off until temperature hits the lower threshold T
t
g
-
T
m
g and the whole process is repeated.Table 1
Training and testing parameters for outer temperature model and heat gain.
Walter’s model
Heat gain
Phase
T
m
a
x
T
m
i
n
T
m
q
(
t
)
Training
15
2
7
600
Testing
16
1
6
700Table 2
Training and testing zones’ target temperatures.
Target temperatures
Zone #
I
II
III
IV
V
VI
VII
VIII
IX
Training
16
18
16
15
14
13
12
16
14
Testing
17
20
17
17
17
14
18
19
21Having collected the data over the period of 24 hours, the measurements of the zones’ temperature, the outer temperature, and the heat gains of each zone are organized into input-output data for the subsystem identification process. For example, for zone I, the outer temperature, the heating gain of zone I, and the temperatures of zones II and IV are the input data to the identification process whereas the temperature of zone I itself is the output data. Note that the number of input signals differs from zone to zone. Zones I, III, VII, and IX use 4 input signals, whereas the rest of the zones use 5 input signals. Based on the input-output data we identify the matricesA, B, C, and D for each zone as described in Section 3. During this process we have to decide on the order of the subsystems. This is achieved by looking at the singular values of the SVD decomposition in (20) and deciding on the number of dominant ones. Figure 6 depicts the singular values for zones I, II, and V. All subsystems exhibit similar behavior regarding the profile of their singular values and therefore we set the order of all subsystems equal to 2.Figure 6
Singular values of zones I, II, and V.Next we run a test on the obtained state-space models. We set the outer temperature parameters as in the second row of Table1 and the target temperatures as in the second row of Table 2. The initial state of the subsystem model was set to -
10
0, whereas the elements of the 42-state vector of the multizone system were set equal to 3 (the initial output temperature). Figure 7 shows the evolution of the actual and predicted temperatures for zones I and IV. For zone I the heater was on (heat gain equal to 700) whereas the heater for zone IV was off. As it is observed, after 5000 samples (appoximately 80 minutes period) the state of the subsystems has converged to one that produces almost the same output as the original system. After this point the WSN nodes can enter in the “decision” phase.Figure 7
Real and predicted temperature of zones I and IV.We should notice at this point that the “original” simulated system (42-state-space model) is also linear and thus the use of linear subspace identification may be questionable for more complex and possibly nonlinear systems. As the simulation results indicate, although the dynamics of each zone are more complex (including the roof temperature,e.g.), a low dimensionality subsystem of order 2 can capture its behavior. For more complex systems we can choose the order of the identified subsystem to be high enough. For nonlinear systems the identification of a time varying system is possible by using a recursive update of the model.Next we simulate the scenario of Figure2. We assume that occupants enter zone I (at time instant 5000) and then after a random period of time they enter zone IV, where they remain for 3600 sec. Zone I is heated since the beginning of the process (heat gain 900 W) until the occupants leave the zone to enter zone IV. Zone IV uses a heater with gain 1200 W, which is turned on (or off) according to the decisions of its node. The target temperatures for zones I and IV are 17°C and 16°C, respectively. For the first set of experiments, the occupation period of zone I is Gamma distributed, that is,(40)
f
X
x
=
β
α
Γ
α
x
α
-
1
e
-
β
x
,with the shape parameter α
=
5 and β such that the expected value of X, E
[
X
]
=
α
/
β, is 3600 and 7200. Decisions were taken every D
=
67 seconds. Figure 8 shows the total comfort cost achieved for various values of the comfort weight C
i. The results are averages over 200 runs of the simulation. The case of C
i
=
0 is the full reactive case where heating of zone IV is postponed until occupants enter the zone. On the other extreme, for high values of C
i. the decision process acts proactively by spending energy in order to reduce the discomfort level.Figure 8
Total comfort cost versusC
i. The occupation time of zone I is Gamma distributed.Figure9 shows the total energy cost in KWh for the two cases of Figure 8. As it is observed, the higher the value of C
i the more the energy consumed since then heating of zone IV starts earlier. The difference between the two curves of Figure 9 is justified by the difference of the means of the occupation period, E
[
X
]
=
3600 and E
[
X
]
=
7200, for the two cases.Figure 9
Total energy cost versusC
i. The occupation time of zone I is Gamma distributed.In Figure10 we plot the temperature profiles for zones I, IV and for two different values of C. As it is observed the value of C
=
0 corresponds to the reactive case; that is the heater of zone IV remains off until occupants are detected in the zone. On the other extreme a large value of C (10000) will force the heater of zone IV to be turned on as soon as possible in an effort to reduce the discomfort penalty.Figure 10
Temperature versus time forC
=
10000 and C
=
0.The scenario of Figure2 was also simulated for Pareto distributed occupation times of zone I. The pdf for Pareto distribution is(41)
f
X
x
=
α
x
m
α
x
α
+
1
,
x
≥
x
m
0
,
x
<
x
mwith expected value(42)
E
X
=
α
x
m
α
-
1
,
α
>
1
∞
,
α
≤
1
.Figure 11 shows the Pareto pdf for four values of x
m (x
m
=
300,900,1800,3600) and α suitably chosen so that the mean value of the occupation period is 7200 in all cases.Figure 11
Pareto pdf for various values ofx
m.Figures12 and 13 show the total comfort cost and the total energy consumed, respectively, for various values of the comfort gain C
i. As it is observed, a 50% improvement on the comfort cost (compared to the reactive case) can be achieved for a value of C equal to 1000. However, there is a performance floor (this is evident for the case of x
m
=
300), which is due to the nature of the Pareto distribution. Most of the samples for the occupation time period are concentrated close to x
m and therefore there is not enough time to preheat the next visiting zone, which in turn implies high discomfort costs. Moreover, the heavy tail of the Pareto distribution causes some extremely high values of the occupation period which result in an increase of the total consumed energy, as it can be observed clearly from Figure 13.Figure 12
Total comfort cost for Pareto distributed occupation period.Figure 13
Total energy cost for Pareto distributed occupation period.Next, we simulate the scenario of Figure3. The occupants move in cascade from zone to zone as it is depicted in Figure 14. The target temperatures for the zones are given in the first row of Table 3, whereas the heat gain of each zone is given in the second row of Table 3. For the occupancy time of each zone we considered an exponential random variable with mean provided in the last row of Table 3. Note that this assumption is the least favorable for adjacent zones due to the memoryless property of the exponential distribution. The total comfort and energy costs for this scenario are plotted in Figure 15 for various values of the comfort weight C. The results are averages over 100 runs. The two curves for each cost correspond to D
=
67 and D
=
167, respectively. As it is observed, considerable comfort gains are achieved for a moderate increase of the consumed energy. In Figure 15 we also show the costs for a “fixed” proactive scenario in which the heating of a zone starts as soon as occupants are detected to its neighbor zone. For example, when occupants enter zone IV, the heater of zone V is turned on, and so on.Table 3
Parameters used in cascade movement scenario.
Zone #
I
IV
V
VIII
IX
VII
III
II
Target  °C
17
16
17
17
16
14
17
18
Power (KW)
0.9
1.2
0.9
0.9
0.9
1.2
0.9
1.2
Mean time × 102
36
27
18
40
54
45
30
20Figure 14
Movement of occupants in cascade (zones I, IV, V, VIII, IX, VI, III, and II).Figure 15
Total energy and comfort costs (cascade movement scenario).Table4 shows the comfort and energy costs obtained when there is no randomness on the occupancy times for various values of the comfort weight C and for two different values of D (D
=
67 and 167). The occupancy times of the zones were set to the mean values provided in the last row of Table 3. As it is observed, extremely low comfort costs can be achieved in this case which implies that the right modelling of the occupancy times is of paramount importance to the process.Table 4
Energy and comfort costs for constant occupancy times.
Comfort weightC
500
1000
1500
2000
D
=
67
Comfort cost
2615
72
67
67
Energy cost
7.75
9.07
9.73
10.20
D
=
167
Comfort cost
2856
72
71
64
Energy cost
7.79
9.11
9.77
10.25Finally, Figure16 shows a sample of the zones’ temperature profile. Note that the temperature of zone II, although zone II is the last zone visited, starts increasing (at a slower rate) much earlier. This is because zone II neighbors zones I and V which are heated in earlier stages of the process.Figure 16
Temperature profiles of zones (cascade movement scenario).


## 4.1. Training Phase

For the deterministic subsystem identification we used 86400 samples with sampling periodT
s
=
1 sec. For the outside temperature we used Walter’s model with parameters given in the first row of Table 1 and a heater gain equal to 600 W. We set the zones’ target temperatures T
t
g to those of Table 2 (first row) with a margin equal to T
m
g
=
1°C. That is, for a specific zone the heater is on until temperature T
t
g
+
T
m
g is reached. After this point the heater is turned off until temperature hits the lower threshold T
t
g
-
T
m
g and the whole process is repeated.Table 1
Training and testing parameters for outer temperature model and heat gain.
Walter’s model
Heat gain
Phase
T
m
a
x
T
m
i
n
T
m
q
(
t
)
Training
15
2
7
600
Testing
16
1
6
700Table 2
Training and testing zones’ target temperatures.
Target temperatures
Zone #
I
II
III
IV
V
VI
VII
VIII
IX
Training
16
18
16
15
14
13
12
16
14
Testing
17
20
17
17
17
14
18
19
21Having collected the data over the period of 24 hours, the measurements of the zones’ temperature, the outer temperature, and the heat gains of each zone are organized into input-output data for the subsystem identification process. For example, for zone I, the outer temperature, the heating gain of zone I, and the temperatures of zones II and IV are the input data to the identification process whereas the temperature of zone I itself is the output data. Note that the number of input signals differs from zone to zone. Zones I, III, VII, and IX use 4 input signals, whereas the rest of the zones use 5 input signals. Based on the input-output data we identify the matricesA, B, C, and D for each zone as described in Section 3. During this process we have to decide on the order of the subsystems. This is achieved by looking at the singular values of the SVD decomposition in (20) and deciding on the number of dominant ones. Figure 6 depicts the singular values for zones I, II, and V. All subsystems exhibit similar behavior regarding the profile of their singular values and therefore we set the order of all subsystems equal to 2.Figure 6
Singular values of zones I, II, and V.Next we run a test on the obtained state-space models. We set the outer temperature parameters as in the second row of Table1 and the target temperatures as in the second row of Table 2. The initial state of the subsystem model was set to -
10
0, whereas the elements of the 42-state vector of the multizone system were set equal to 3 (the initial output temperature). Figure 7 shows the evolution of the actual and predicted temperatures for zones I and IV. For zone I the heater was on (heat gain equal to 700) whereas the heater for zone IV was off. As it is observed, after 5000 samples (appoximately 80 minutes period) the state of the subsystems has converged to one that produces almost the same output as the original system. After this point the WSN nodes can enter in the “decision” phase.Figure 7
Real and predicted temperature of zones I and IV.We should notice at this point that the “original” simulated system (42-state-space model) is also linear and thus the use of linear subspace identification may be questionable for more complex and possibly nonlinear systems. As the simulation results indicate, although the dynamics of each zone are more complex (including the roof temperature,e.g.), a low dimensionality subsystem of order 2 can capture its behavior. For more complex systems we can choose the order of the identified subsystem to be high enough. For nonlinear systems the identification of a time varying system is possible by using a recursive update of the model.Next we simulate the scenario of Figure2. We assume that occupants enter zone I (at time instant 5000) and then after a random period of time they enter zone IV, where they remain for 3600 sec. Zone I is heated since the beginning of the process (heat gain 900 W) until the occupants leave the zone to enter zone IV. Zone IV uses a heater with gain 1200 W, which is turned on (or off) according to the decisions of its node. The target temperatures for zones I and IV are 17°C and 16°C, respectively. For the first set of experiments, the occupation period of zone I is Gamma distributed, that is,(40)
f
X
x
=
β
α
Γ
α
x
α
-
1
e
-
β
x
,with the shape parameter α
=
5 and β such that the expected value of X, E
[
X
]
=
α
/
β, is 3600 and 7200. Decisions were taken every D
=
67 seconds. Figure 8 shows the total comfort cost achieved for various values of the comfort weight C
i. The results are averages over 200 runs of the simulation. The case of C
i
=
0 is the full reactive case where heating of zone IV is postponed until occupants enter the zone. On the other extreme, for high values of C
i. the decision process acts proactively by spending energy in order to reduce the discomfort level.Figure 8
Total comfort cost versusC
i. The occupation time of zone I is Gamma distributed.Figure9 shows the total energy cost in KWh for the two cases of Figure 8. As it is observed, the higher the value of C
i the more the energy consumed since then heating of zone IV starts earlier. The difference between the two curves of Figure 9 is justified by the difference of the means of the occupation period, E
[
X
]
=
3600 and E
[
X
]
=
7200, for the two cases.Figure 9
Total energy cost versusC
i. The occupation time of zone I is Gamma distributed.In Figure10 we plot the temperature profiles for zones I, IV and for two different values of C. As it is observed the value of C
=
0 corresponds to the reactive case; that is the heater of zone IV remains off until occupants are detected in the zone. On the other extreme a large value of C (10000) will force the heater of zone IV to be turned on as soon as possible in an effort to reduce the discomfort penalty.Figure 10
Temperature versus time forC
=
10000 and C
=
0.The scenario of Figure2 was also simulated for Pareto distributed occupation times of zone I. The pdf for Pareto distribution is(41)
f
X
x
=
α
x
m
α
x
α
+
1
,
x
≥
x
m
0
,
x
<
x
mwith expected value(42)
E
X
=
α
x
m
α
-
1
,
α
>
1
∞
,
α
≤
1
.Figure 11 shows the Pareto pdf for four values of x
m (x
m
=
300,900,1800,3600) and α suitably chosen so that the mean value of the occupation period is 7200 in all cases.Figure 11
Pareto pdf for various values ofx
m.Figures12 and 13 show the total comfort cost and the total energy consumed, respectively, for various values of the comfort gain C
i. As it is observed, a 50% improvement on the comfort cost (compared to the reactive case) can be achieved for a value of C equal to 1000. However, there is a performance floor (this is evident for the case of x
m
=
300), which is due to the nature of the Pareto distribution. Most of the samples for the occupation time period are concentrated close to x
m and therefore there is not enough time to preheat the next visiting zone, which in turn implies high discomfort costs. Moreover, the heavy tail of the Pareto distribution causes some extremely high values of the occupation period which result in an increase of the total consumed energy, as it can be observed clearly from Figure 13.Figure 12
Total comfort cost for Pareto distributed occupation period.Figure 13
Total energy cost for Pareto distributed occupation period.Next, we simulate the scenario of Figure3. The occupants move in cascade from zone to zone as it is depicted in Figure 14. The target temperatures for the zones are given in the first row of Table 3, whereas the heat gain of each zone is given in the second row of Table 3. For the occupancy time of each zone we considered an exponential random variable with mean provided in the last row of Table 3. Note that this assumption is the least favorable for adjacent zones due to the memoryless property of the exponential distribution. The total comfort and energy costs for this scenario are plotted in Figure 15 for various values of the comfort weight C. The results are averages over 100 runs. The two curves for each cost correspond to D
=
67 and D
=
167, respectively. As it is observed, considerable comfort gains are achieved for a moderate increase of the consumed energy. In Figure 15 we also show the costs for a “fixed” proactive scenario in which the heating of a zone starts as soon as occupants are detected to its neighbor zone. For example, when occupants enter zone IV, the heater of zone V is turned on, and so on.Table 3
Parameters used in cascade movement scenario.
Zone #
I
IV
V
VIII
IX
VII
III
II
Target  °C
17
16
17
17
16
14
17
18
Power (KW)
0.9
1.2
0.9
0.9
0.9
1.2
0.9
1.2
Mean time × 102
36
27
18
40
54
45
30
20Figure 14
Movement of occupants in cascade (zones I, IV, V, VIII, IX, VI, III, and II).Figure 15
Total energy and comfort costs (cascade movement scenario).Table4 shows the comfort and energy costs obtained when there is no randomness on the occupancy times for various values of the comfort weight C and for two different values of D (D
=
67 and 167). The occupancy times of the zones were set to the mean values provided in the last row of Table 3. As it is observed, extremely low comfort costs can be achieved in this case which implies that the right modelling of the occupancy times is of paramount importance to the process.Table 4
Energy and comfort costs for constant occupancy times.
Comfort weightC
500
1000
1500
2000
D
=
67
Comfort cost
2615
72
67
67
Energy cost
7.75
9.07
9.73
10.20
D
=
167
Comfort cost
2856
72
71
64
Energy cost
7.79
9.11
9.77
10.25Finally, Figure16 shows a sample of the zones’ temperature profile. Note that the temperature of zone II, although zone II is the last zone visited, starts increasing (at a slower rate) much earlier. This is because zone II neighbors zones I and V which are heated in earlier stages of the process.Figure 16
Temperature profiles of zones (cascade movement scenario).


## 5. Conclusions and Future Work

A method that fuses temperature and occupancy predictions in an effort to balance the thermal comfort condition and the consumed energy in a multizone HVAC system was presented. Temperature predictions are based on a subspace identification technique used to model the thermal dynamics of each zone independently. This technique is quite robust and accurate predictions are possible after a convergence period which may last 1-2 hours. With regard to the occupancy predictions the decision process makes use of the distribution of the first passage time from an occupied zone to an unoccupied one. Simulation results for different first passage distributions demonstrated the effectiveness of the proposed technique.Several extensions and modifications of the proposed method are possible. The comfort parameterC
i may depend on the zone, the number of occupants, and the current status of the zones. Moreover, this parameter may be time variable, that is, different values may be used for day and night hours. Regarding the decision process itself, there is no need that it be executed at regular time epochs. If the energy consumed by the wireless sensor nodes is an issue, we may take decisions at irregular time epochs, depending on the occupancy predictions of the zones. The computation of the risks, used by the decision process, may take into consideration additional parameters that affect energy consumption and/or the thermal comfort of the occupants. For example, open windows situation may be easily detected and incorporated suitably in the decision process.



---

*Source: 5074846-2016-03-02.xml*