# A Subspace Identification Method for Detecting Abnormal Behavior in HVAC Systems

Authors: Dimitris Sklavounos; Evangelos Zervas; Odysseas Tsakiridis; John Stonham
Journal: Journal of Energy
 (2015)

---

## Abstract

A method for the detection of abnormal behavior in HVAC systems is presented. The method combines deterministic subspace identification for each zone independently to create a system model that produces the anticipated zone’s temperature and the sequential test CUSUM algorithm to detect drifts of the rate of change of the difference between the real and the anticipated measurements. Simulation results regarding the detection of infiltration heat losses and the detection of exogenous heat gains such as fire demonstrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Energy consumption control in buildings remains nowadays one of the most important issues in the total energy savings. About 40% of the total energy consumption is due to energy requirements for buildings and this percentage tends to be increasing from about 0.5% to 5% per year [1]. Heat ventilation and air conditioning systems (HVAC) have become one of the main and widely equipped environmental technologies for the buildings. HVAC systems of the commercial and residential buildings consume about 57% of the required energy and they are wasting more than 20% of energy due to various faults, insufficient control, or improper positioning. Methods of faults detection (FDD) as well as abnormal situations detection are potentially able to save about 10–40% of the HVAC energy consumption [2]. Thus, sophisticated building automation systems, which utilize such methods of detection and improve the control mechanisms of HVAC systems, can significantly contribute to the total energy savings. To this end, the fault detection diagnosis as well as the detection of abnormal situations in the HVAC systems is a very challenging research area.As the wireless sensor networks’ (WSNs) technology has been improved and widely applied, it has been also utilized in the HVAC systems contributing to more sophisticated control. A promising and effective way of implementing a HVAC system in a building area with respect to improved energy control is the multizone model approach. In such a model, each area of a heating/cooling space such as room, corridor, or others is treated as a single zone. With the use of a wireless node installed in each zone, a wireless network of scattered nodes (WSN) could be formed where each node would employ a variety of sensors for temperature, light, occupancy, and so forth. An appropriate utilization of a WSN would significantly contribute to the control of the energy consumption of the building. WSNs are easily installed and deployed allowing cost effective retrofits. The HVAC application requires a sensor network to process data cooperatively and combines information from multiple sources. In traditional centralized systems, measurements collected by sensor nodes are relayed to a central unit for further processing. In the decentralized or distributed systems, all wireless sensor nodes are autonomous and they perform their assigned task locally utilizing information from neighboring nodes.A hybrid scheme of operation is used in the proposed algorithm. In a first phase, called “training phase,” the nodes installed in each room (zone) of the multizone system send their readings to a central computing unit. These readings are the measurements of the zone temperature as well as the current power of the heater that exists in each zone. Upon reception, the central unit arranges the received data to an input-output form in order to run a subspace identification (SID) process for each zone. The inputs, for each zone SID process, are the temperature measurements of its adjacent zones as well as the power profile of the heater of the zone. The output of the SID process is the zone’s temperature itself. The central unit identifies a linear state space system for each zone and its parameters are communicated back to the wireless sensor node that monitors the zone. After this point, the system enters in the second phase (detection phase) and the operation is turned to a decentralized mode. Each node collects temperature measurements by its surrounding sensor nodes and power measurements of the zone’s heater. Based on these measurements the identified state space subsystem predicts the temperature of the zone. A suitable detection algorithm is then applied to detect possible deviations of the predicted values from real measurements. Deviations may be due to high infiltration heat gains or other exogenous factors such as fire. The detection phase is spit into cycles and it utilizes the CUSUM algorithm to detect possible deviations from the normal operation.The rest of the paper is organized as follows: Section2 presents related work on the detection of faults and abnormalities as well as the use of SID to HVAC systems. Section 3 describes the system model and the subspace identification process as well as the CUSUM algorithm that is used to detect abnormal system operation. 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

Several methods have been proposed for the detection of abnormal energy consumption or fault detection diagnosis (FDD) in HVAC systems and they are divided in two main categories: the statistical methods and the computational approaches. The statistical methods are mainly based on fault detection algorithms that compare data under normal operation conditions with the data under current conditions in order to detect any abnormal behavior. The authors in [3] proposed a method that is based on the principal component analysis (PCA) detection of sensor faults in air handling units (AHU). The Q-statistic method is used for the sensor fault detection, and in addition with the use of the Q contribution plot, the faulty sensors can be isolated. In [4], a statistical based method of detecting abnormal energy consumption in buildings is proposed, including the detection of HVAC abnormalities causing outliers. In a first stage, features like the average daily energy consumption or peak demands for a day are determined using data like the total energy consumption of the building. Then these features are sorted according to similar energy consumption and thus groups of days with similar energy consumption are formed. Thereafter, an outlier identifier is applied to determine features of the same day type with significant difference from the normal ones and if detected, a modified z-score is used to determine the amount and direction of the variation. In [5], the authors estimated the appropriate power consumption by approximating the minimum cooling demands of a building (National Taiwan University) and comparing them with the real cooling supply. The results showed a high discrepancy and the authors proposed two types of statistical methods to reduce the energy consumption: polynomial regression and feature based regression are the methods used to model the behavior in the building and a Hampel identifier is applied to test the consistency of the data.In the category of the computational approaches earlier works have introduced computer simulations as embedded mechanisms within the control methods of energy consumption of HVAC systems. In [6] the authors, based on the idea of encapsulating simulation programs within building energy management systems (BEMs), proposed a prototype simulation-assisted controller with an embedded simulation program in order to provide real time control decisions. In the same category also falls the work in [7], where a WSN was utilized for the abnormal situation detection. Each node of the WSN, covering a single zone, act as a controller (PI), tuning the heating supply to a predetermined temperature value. A lumped capacity model was used to predict the hypothetical normal operation and the CUSUM sequential algorithm detected possible divergences of the energy consumption from the anticipated one. In [8] a method using a multilevel fault detection diagnosis (FDD) algorithm is proposed, with an energy description of all units in a HVAC and a spatial temporal partition strategy, as the two main elements of the method. The energy performance signals of the HVAC units are becoming the inputs to the FDD algorithm, and possible hardware faults within the HVAC system are captured. The work in [9] proposes a method that detects actuator faults in a HVAC system. It is a software based fault detection diagnosis mechanism with a two-tiered detection approach. At the first tier, the method utilizes a quantitative model-based approach, which relies on a simple thermodynamic model and it does not require full knowledge of the system model. At the second tier, a qualitative model-based approach is utilized which, based on the air temperature, provides a quick decision whether an actuator is working properly or not. In [10] a Model Predictive Controller is presented which uses both weather forecasts as well as the thermal model of a building in order to maintain indoor temperature independent of the outdoor conditions. An accurate model of the building was indispensable and thus they created a simplified model of the crittall type ceiling of radiant heating and applied a subspace identification algorithm giving the appropriate inputs. In [11] a strategy is presented consisting of two schemes, the FDD of HVAC systems and the sensors fault detection. In the first scheme the indication of each system performance is taken by one or more performance indices (PIs) which they are validated from the actual measurements by the use of regression models. The detection and diagnosis of faulty sensors are achieved with the use of a PCA method. In [12] the authors have combined the benefits of the model predictive control (MPC) technique as well as building simulation software such as TRNSYS, Energy plus, and so forth, in order to create a physical model, as close as possible to a real building. To achieve this, they applied a subspace identification method which is appropriate to identify a multiple input multiple output (MIMO) system. In [13] the authors have focused on the preventive maintenance of the HVAC systems and they proposed a fault detection method that combines the model FDD method and a support vector machine classifier (SVM). The authors worked on the detection of components sensitive to faults. Using computer simulations they investigated three major faults such as recirculation damper stuck, the block of the cooling coil, and the decreasing of the supply fan.


## 3. System Modeling

### 3.1. Detection of Abnormal System Behavior

As it has been already mentioned, we consider a multizone system with a WSN deployed consisting of temperature sensor nodes. Each wireless node measures the temperature of the zone it covers and conveys this information to its neighboring nodes. Additionally, the nodes have the functionality to detect abnormal operation, for example, slower temperature rising than the anticipated one due to open windows during winter, or high temperature values due to the onset of a fire, and signal it to an operation center. The detection mechanism relies on knowledge of the surrounding zones’ temperatures and the dynamics of the covered zone. The dynamics of each zone are learned during a training period using the subspace identification procedure presented in the next subsection.In general, the assumption is made that each zone is represented by a discrete model of the form(1)
x
k
+
1
=
F
x
k
,
u
k
,
ψ
y
k
=
G
x
k
,
u
k
,
ζ
,where x
k is the state vector at discrete time k, u
k is a vector of external inputs, y
k is the predicted zone temperature, and ψ, ζ are specific parameters to the systems F
· and G
·, respectively. Knowing the dynamics of the system and the external inputs u
k, the node is able to predict the temperature of the zone y
k. The inputs u
k in this case are the temperatures of the surrounding zones, the outdoor temperature, and the power of heaters located in the zone. Comparing the predicted values to the actual ones, as measured by the temperature sensors, possible changes in the dynamics of the system can be detected, which signal an abnormal operation. Optimal detection theory deals with the problem of detection of changes in the distribution of measurement data. Control charts [14] and the CUSUM algorithm [15] are some of the simplest and most applied solutions. Alternatively, the subsystem model can be periodically reidentified looking for possible changes in the parameter space ψ
,
ζ. However, such a detection process is computationally demanding and overloads the wireless nodes. Moreover, more data are needed to be processed to draw safe conclusions and thus large delays are unavoidable. Shewhart charts cannot detect small shifts; that is, the probability of detecting small shifts is rather small. On the contrary CUSUM charts can detect easily small systematic shifts but their response to large shifts is relatively slow. For all these reasons, in this work, the CUSUM technique is used as the basic change detection algorithm since even small deviations from strict operating requirements should be detected. The detection process will be based on the rate of variation of the prediction error between the real and the predicted temperature. To this end the sequence of prediction errors E
k
=
Y
k
r
-
Y
k is defined where Y
k denotes the predicted temperature process and Y
k
r is the “real” temperature process of the zone. The derivative of this discrete process is the random variable Z
k, where(2)
Z
k
=
E
k
-
E
k
-
1
.The assumption is made that Z
k have density f
(
z
k
;
q
n
,
σ
) for k
=
1
,
…
,
τ
-
1 and density f
(
z
k
;
q
a
,
σ
) for k
≥
1 where the parameter q
n is known (or it is estimated) and q
a and σ are generally unknown. The time index τ signals an event that changes the distribution of the measurements. In terms of the proposed algorithm, q
n is the mean of the rate of the prediction error in normal operation and σ
2 is the variance of the rate due to uncertainties of the measurement devices. σ is a nuisance parameter and it is generally unknown. The parameter q
a denotes the mean of the rate of the prediction error in case of an abnormal situation, and it is considered unknown. The parameter τ is the time index when a change of densities occurs and sequential tests deal with the detection of this change. Although, q
a is in general a function of time k, the detection algorithm will be applied using the random variables Q
k
=
Z
k
-
q
n. Thus, under normal operation (null hypothesis) Q
k has a constant mean value equal to 0, whereas under the alternative hypothesis Q
k has a nonzero unknown mean value q
=
q
a
-
q
n. Note that the measurement noise is only due to the sensors and not due to exogenous factors such as lamps and the presence of people that may influence the detection process. These sources of noise affect both the estimate of q
n, the mean of the rate of the prediction error on normal operation, and the drift of the test statistic (see below) so that one cannot draw safe conclusions about an abnormal situation. Nevertheless, in the simulation section, we include an experiment with such noise processes modeled, primarily to demonstrate the effectiveness of the SID method and to provide directions for future research.As it has already been stated, one of the most promising algorithms to sequentially detect the change is the CUSUM test. Gombay and Serban [16] adapted Page’s CUSUM test for change detection in the presence of nuisance parameters. Gombay proposed statistics based on the efficient score (Rao’s statistics), on the maximum likelihood estimator (Wald’s statistics), or on the log likelihood ratio. The efficient score vector is defined as(3)
V
k
q
,
σ
=
∑
i
=
1
k
∇
ν
log
⁡
f
q
i
;
q
,
σ
,
v
=
q
,
σ
.As it can be proved, if the density f
(
·
) belongs to the exponential family, that is, Gaussian, then if some regularity conditions hold under the null hypothesis (normal operation), there exists a Wiener process W
(
t
) that approximates(4)
W
k
=
Γ
-
1
/
2
q
,
σ
V
k
q
,
σ
^
k
,where σ
^
k is the maximum likelihood estimation of σ and Γ
q
,
σ is the Fisher information matrix.The test statisticW
k in (4) can be used to check if a change in densities has occurred at some time instant τ
≤
k. Under the alternative hypothesis (abnormal behavior) this statistic drifts for k
≥
τ with the size of the drift proportional to the rate at which the test statistic moves in the direction of the alternative density. Moreover, in order to make decisions after n observations have been obtained, the following result is used (Darling and Erdös [17])(5)
exp
⁡
⁡
-
e
-
s
=
lim
n
→
∞
⁡
P
{
α
(
log
⁡
⁡
n
)
max
⁡
1
<
k
≤
n
⁡
k
-
1
/
2
W
k
≤
s
+
b
(
log
⁡
⁡
n
)
}
,where(6)
a
x
=
(
2
log
⁡
⁡
x
)
1
/
2
,
b
x
=
2
log
⁡
⁡
x
+
0.5
log
⁡
⁡
log
⁡
⁡
x
-
0.5
log
⁡
⁡
π
.To make use of this result a false alarm rate f is set, that is, f
=
0.001, where 1
-
f
=
exp
⁡
⁡
(
-
e
-
s
) and the threshold is computed(7)
T
f
=
2
log
⁡
⁡
log
⁡
⁡
n
-
1
/
2
·
-
log
⁡
⁡
-
log
⁡
⁡
1
-
f
+
2
log
⁡
⁡
log
⁡
⁡
n
+
0.5
log
⁡
⁡
log
⁡
⁡
log
⁡
⁡
n
-
0.5
log
⁡
⁡
π
.Then, the alternative hypothesis (abnormal operation) is supported by the data at the first k, when(8)
k
-
1
/
2
W
k
≥
T
f
.If no such k exists for k
≤
n the hypothesis of normal operation is not rejected. For n
=
300 (that is 5 min with sampling period 1 sec) and the two indicative values of f
=
0.01 and f
=
0.001
T
f
=
4.17 and T
f
=
5.42 are obtained, respectively.In what follows the assumption is made that all measurementsQ
i, i
≥
1, are independent random variables. In this case the test statistic in (4) takes the form(9)
k
-
1
/
2
W
k
=
k
-
1
/
2
∑
i
=
1
k
Q
i
∑
i
=
1
k
Q
i
2
/
k
1
/
2
.Under the alternative, the drift of k
-
1
/
2
W
k after change is (10)
k
-
1
/
2
k
-
τ
-
1
q
σ
2
+
k
-
t
-
1
/
k
q
2
1
/
2
.Figure 1 shows the drift for τ
=
1, ν
=
300, σ
=
5, q
=
1 (blue), q
=
2 (green), and q
=
3 (red). As it is observed the greater the excess difference of the means of the rates (
q
=
q
a
-
q
n
) the largest the slope of the drift.Figure 1
Drift of test statistic forσ
=
5 and q
=
1,2, and 3.The only parameter needed to run the detection algorithm isq
n, that is, the mean of the rate of the prediction error during normal operation. This parameter is estimated for each node using the model described by (1).


### 3.2. Deterministic Subspace Identification

Each zone is treated separately and is represented by a discrete time, linear, time invariant, and state space model. That is,F
· and G
· are linear functions and for a specific zone we can write(11)
x
k
+
1
=
A
x
k
+
B
u
k
+
w
k
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 only one sensor node is needed and thus y
k becomes 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 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, switching lamps on and off, and so forth. For the rest of the paper the effects of state and output noise are neglected and thus a deterministic identification problem is obtained, which is the computation of the matrices A, B, C, and D from the given input-output data.Following the notation and the derivation in [18] the block Hankel matrices are defined:(12)
U
p
=
U
0
∣
i
-
1
=
u
0
u
1
⋯
u
j
-
1
u
1
u
2
⋯
u
j
⋮
⋮
⋮
u
i
-
1
⋯
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
+
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 the block Hankel matrix is obtained(13)
U
0
∣
2
i
-
1
=
U
0
∣
i
-
1
U
i
∣
2
i
-
1
=
U
p
U
fwhich can also partitioned as(14)
U
0
∣
2
i
-
1
=
U
0
∣
i
U
i
+
1
∣
2
i
-
1
=
U
p
+
U
fby moving the “present” time one step ahead to time i
+
1. Similarly, the output block Hankel matrices Y
0
∣
2
i
-
1, Y
p, Y
f, Y
p
+, Y
f
- and the input-output data matrices are defined(15)
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(16)
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 (11), the extended observability matrix Γ
i and the reversed extended controllability matrix Δ
i are associated where(17)
Γ
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 is based on determining the state sequence X
f and the extended observability matrix Γ
i directly from the input-output data u
k, y
k and then based on them, in a next step, extracting the matrices A, B, C, and D. To this end the oblique projection of the row space of Y
f along the row space of U
f on the row space of W
p is defined:(18)
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 the matrix B; that is,(19)
A
B
⊥
=
I
-
A
B
T
B
B
T
†
B
.As it is proved in [19] the matrix O
i is equal to the product of the extended observability matrix and the “future” state vector, (20)
O
i
=
Γ
i
X
f
.Having the matrix O
i its singular value decomposition (SVD) is computed(21)
O
i
=
U
1
U
2
S
1
0
0
0
V
1
T
V
2
T
=
U
1
S
1
V
1
Tand the extended observation matrix Γ
i and the state vector X
f are obtained(22)
Γ
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 [18] to extract the matrices A, B, C, and D uses in addition the oblique projection(23)
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, the following equation is obtained: (24)
X
i
+
1
=
Γ
i
-
1
†
O
i
-
1and the matrices A, B, C, and D are obtained by solving (in a least square fashion) the system(25)
X
i
+
1
Y
i
∣
i
=
A
B
C
D
X
i
U
i
∣
i
.


### 3.3. Proposed Algorithm for Abnormal Operation Detection

As it has been already stated deterministic subspace system identification and the CUSUM algorithm are the key ingredients in the proposed algorithm, which undergoes two phases: the training and the detection phase.During the first phase (training) all the nodes in the structure report their 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 detection 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 that several variations in the input signals to excite the modes of the system are present. In the simulation part a training period of 24 h, that is, 86400 samples with sampling period 1 sec are considered. After receiving all the measurements, the central system arranges them to input-output data for each zone. For example, in the case of a zone namedA, with neighbors zones named B, C, and D and the exterior to the building zone named E, then the temperature values of zone A are considered the output of the system to be identified, whereas the temperatures of the 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. Then, a deterministic subspace identification process is run for each zone as described in the previous subsection. The relevant matrices A, B, C, and D of the state space model of each zone are communicated back to the sensor node that is responsible for monitoring the zone. Figure 2 depicts the flow diagram of the training phase. After reception of these matrices, the sensor node enters the detection phase.Figure 2
Training phase steps of the algorithm.The detection phase is split into cycles of possibly unequal lengths. In the beginning of each cycle the valueq
n, that is, the mean of the rate of the prediction error is estimated using a few samples. For the simulation part in this paper short periods of 1 min (60 samples) are used to estimate the value q
n. Then, for the rest of the cycle, the CUSUM algorithm is run to detect possible deviations from the normal operation. The length of the detection interval is quite large, in the order of 1 hour, and it should depend on the time of the day. During periods (of the day or night) with small (large) expected changes of the exterior temperature, the detection interval can be increased (decreased).Note that a positive drift of the test statistic (q
>
0) designates an abnormal increase of temperature that may be due to the onset of a fire, whereas a negative drift of the test statistic (q
<
0) designates leakage of heat that may be due to open windows during winter, or malfunctioning of the heater. The detection phase flow diagram is depicted in Figure 3.Figure 3
Detection phase steps of the algorithm.


## 3.1. Detection of Abnormal System Behavior

As it has been already mentioned, we consider a multizone system with a WSN deployed consisting of temperature sensor nodes. Each wireless node measures the temperature of the zone it covers and conveys this information to its neighboring nodes. Additionally, the nodes have the functionality to detect abnormal operation, for example, slower temperature rising than the anticipated one due to open windows during winter, or high temperature values due to the onset of a fire, and signal it to an operation center. The detection mechanism relies on knowledge of the surrounding zones’ temperatures and the dynamics of the covered zone. The dynamics of each zone are learned during a training period using the subspace identification procedure presented in the next subsection.In general, the assumption is made that each zone is represented by a discrete model of the form(1)
x
k
+
1
=
F
x
k
,
u
k
,
ψ
y
k
=
G
x
k
,
u
k
,
ζ
,where x
k is the state vector at discrete time k, u
k is a vector of external inputs, y
k is the predicted zone temperature, and ψ, ζ are specific parameters to the systems F
· and G
·, respectively. Knowing the dynamics of the system and the external inputs u
k, the node is able to predict the temperature of the zone y
k. The inputs u
k in this case are the temperatures of the surrounding zones, the outdoor temperature, and the power of heaters located in the zone. Comparing the predicted values to the actual ones, as measured by the temperature sensors, possible changes in the dynamics of the system can be detected, which signal an abnormal operation. Optimal detection theory deals with the problem of detection of changes in the distribution of measurement data. Control charts [14] and the CUSUM algorithm [15] are some of the simplest and most applied solutions. Alternatively, the subsystem model can be periodically reidentified looking for possible changes in the parameter space ψ
,
ζ. However, such a detection process is computationally demanding and overloads the wireless nodes. Moreover, more data are needed to be processed to draw safe conclusions and thus large delays are unavoidable. Shewhart charts cannot detect small shifts; that is, the probability of detecting small shifts is rather small. On the contrary CUSUM charts can detect easily small systematic shifts but their response to large shifts is relatively slow. For all these reasons, in this work, the CUSUM technique is used as the basic change detection algorithm since even small deviations from strict operating requirements should be detected. The detection process will be based on the rate of variation of the prediction error between the real and the predicted temperature. To this end the sequence of prediction errors E
k
=
Y
k
r
-
Y
k is defined where Y
k denotes the predicted temperature process and Y
k
r is the “real” temperature process of the zone. The derivative of this discrete process is the random variable Z
k, where(2)
Z
k
=
E
k
-
E
k
-
1
.The assumption is made that Z
k have density f
(
z
k
;
q
n
,
σ
) for k
=
1
,
…
,
τ
-
1 and density f
(
z
k
;
q
a
,
σ
) for k
≥
1 where the parameter q
n is known (or it is estimated) and q
a and σ are generally unknown. The time index τ signals an event that changes the distribution of the measurements. In terms of the proposed algorithm, q
n is the mean of the rate of the prediction error in normal operation and σ
2 is the variance of the rate due to uncertainties of the measurement devices. σ is a nuisance parameter and it is generally unknown. The parameter q
a denotes the mean of the rate of the prediction error in case of an abnormal situation, and it is considered unknown. The parameter τ is the time index when a change of densities occurs and sequential tests deal with the detection of this change. Although, q
a is in general a function of time k, the detection algorithm will be applied using the random variables Q
k
=
Z
k
-
q
n. Thus, under normal operation (null hypothesis) Q
k has a constant mean value equal to 0, whereas under the alternative hypothesis Q
k has a nonzero unknown mean value q
=
q
a
-
q
n. Note that the measurement noise is only due to the sensors and not due to exogenous factors such as lamps and the presence of people that may influence the detection process. These sources of noise affect both the estimate of q
n, the mean of the rate of the prediction error on normal operation, and the drift of the test statistic (see below) so that one cannot draw safe conclusions about an abnormal situation. Nevertheless, in the simulation section, we include an experiment with such noise processes modeled, primarily to demonstrate the effectiveness of the SID method and to provide directions for future research.As it has already been stated, one of the most promising algorithms to sequentially detect the change is the CUSUM test. Gombay and Serban [16] adapted Page’s CUSUM test for change detection in the presence of nuisance parameters. Gombay proposed statistics based on the efficient score (Rao’s statistics), on the maximum likelihood estimator (Wald’s statistics), or on the log likelihood ratio. The efficient score vector is defined as(3)
V
k
q
,
σ
=
∑
i
=
1
k
∇
ν
log
⁡
f
q
i
;
q
,
σ
,
v
=
q
,
σ
.As it can be proved, if the density f
(
·
) belongs to the exponential family, that is, Gaussian, then if some regularity conditions hold under the null hypothesis (normal operation), there exists a Wiener process W
(
t
) that approximates(4)
W
k
=
Γ
-
1
/
2
q
,
σ
V
k
q
,
σ
^
k
,where σ
^
k is the maximum likelihood estimation of σ and Γ
q
,
σ is the Fisher information matrix.The test statisticW
k in (4) can be used to check if a change in densities has occurred at some time instant τ
≤
k. Under the alternative hypothesis (abnormal behavior) this statistic drifts for k
≥
τ with the size of the drift proportional to the rate at which the test statistic moves in the direction of the alternative density. Moreover, in order to make decisions after n observations have been obtained, the following result is used (Darling and Erdös [17])(5)
exp
⁡
⁡
-
e
-
s
=
lim
n
→
∞
⁡
P
{
α
(
log
⁡
⁡
n
)
max
⁡
1
<
k
≤
n
⁡
k
-
1
/
2
W
k
≤
s
+
b
(
log
⁡
⁡
n
)
}
,where(6)
a
x
=
(
2
log
⁡
⁡
x
)
1
/
2
,
b
x
=
2
log
⁡
⁡
x
+
0.5
log
⁡
⁡
log
⁡
⁡
x
-
0.5
log
⁡
⁡
π
.To make use of this result a false alarm rate f is set, that is, f
=
0.001, where 1
-
f
=
exp
⁡
⁡
(
-
e
-
s
) and the threshold is computed(7)
T
f
=
2
log
⁡
⁡
log
⁡
⁡
n
-
1
/
2
·
-
log
⁡
⁡
-
log
⁡
⁡
1
-
f
+
2
log
⁡
⁡
log
⁡
⁡
n
+
0.5
log
⁡
⁡
log
⁡
⁡
log
⁡
⁡
n
-
0.5
log
⁡
⁡
π
.Then, the alternative hypothesis (abnormal operation) is supported by the data at the first k, when(8)
k
-
1
/
2
W
k
≥
T
f
.If no such k exists for k
≤
n the hypothesis of normal operation is not rejected. For n
=
300 (that is 5 min with sampling period 1 sec) and the two indicative values of f
=
0.01 and f
=
0.001
T
f
=
4.17 and T
f
=
5.42 are obtained, respectively.In what follows the assumption is made that all measurementsQ
i, i
≥
1, are independent random variables. In this case the test statistic in (4) takes the form(9)
k
-
1
/
2
W
k
=
k
-
1
/
2
∑
i
=
1
k
Q
i
∑
i
=
1
k
Q
i
2
/
k
1
/
2
.Under the alternative, the drift of k
-
1
/
2
W
k after change is (10)
k
-
1
/
2
k
-
τ
-
1
q
σ
2
+
k
-
t
-
1
/
k
q
2
1
/
2
.Figure 1 shows the drift for τ
=
1, ν
=
300, σ
=
5, q
=
1 (blue), q
=
2 (green), and q
=
3 (red). As it is observed the greater the excess difference of the means of the rates (
q
=
q
a
-
q
n
) the largest the slope of the drift.Figure 1
Drift of test statistic forσ
=
5 and q
=
1,2, and 3.The only parameter needed to run the detection algorithm isq
n, that is, the mean of the rate of the prediction error during normal operation. This parameter is estimated for each node using the model described by (1).


## 3.2. Deterministic Subspace Identification

Each zone is treated separately and is represented by a discrete time, linear, time invariant, and state space model. That is,F
· and G
· are linear functions and for a specific zone we can write(11)
x
k
+
1
=
A
x
k
+
B
u
k
+
w
k
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 only one sensor node is needed and thus y
k becomes 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 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, switching lamps on and off, and so forth. For the rest of the paper the effects of state and output noise are neglected and thus a deterministic identification problem is obtained, which is the computation of the matrices A, B, C, and D from the given input-output data.Following the notation and the derivation in [18] the block Hankel matrices are defined:(12)
U
p
=
U
0
∣
i
-
1
=
u
0
u
1
⋯
u
j
-
1
u
1
u
2
⋯
u
j
⋮
⋮
⋮
u
i
-
1
⋯
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
+
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 the block Hankel matrix is obtained(13)
U
0
∣
2
i
-
1
=
U
0
∣
i
-
1
U
i
∣
2
i
-
1
=
U
p
U
fwhich can also partitioned as(14)
U
0
∣
2
i
-
1
=
U
0
∣
i
U
i
+
1
∣
2
i
-
1
=
U
p
+
U
fby moving the “present” time one step ahead to time i
+
1. Similarly, the output block Hankel matrices Y
0
∣
2
i
-
1, Y
p, Y
f, Y
p
+, Y
f
- and the input-output data matrices are defined(15)
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(16)
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 (11), the extended observability matrix Γ
i and the reversed extended controllability matrix Δ
i are associated where(17)
Γ
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 is based on determining the state sequence X
f and the extended observability matrix Γ
i directly from the input-output data u
k, y
k and then based on them, in a next step, extracting the matrices A, B, C, and D. To this end the oblique projection of the row space of Y
f along the row space of U
f on the row space of W
p is defined:(18)
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 the matrix B; that is,(19)
A
B
⊥
=
I
-
A
B
T
B
B
T
†
B
.As it is proved in [19] the matrix O
i is equal to the product of the extended observability matrix and the “future” state vector, (20)
O
i
=
Γ
i
X
f
.Having the matrix O
i its singular value decomposition (SVD) is computed(21)
O
i
=
U
1
U
2
S
1
0
0
0
V
1
T
V
2
T
=
U
1
S
1
V
1
Tand the extended observation matrix Γ
i and the state vector X
f are obtained(22)
Γ
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 [18] to extract the matrices A, B, C, and D uses in addition the oblique projection(23)
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, the following equation is obtained: (24)
X
i
+
1
=
Γ
i
-
1
†
O
i
-
1and the matrices A, B, C, and D are obtained by solving (in a least square fashion) the system(25)
X
i
+
1
Y
i
∣
i
=
A
B
C
D
X
i
U
i
∣
i
.


## 3.3. Proposed Algorithm for Abnormal Operation Detection

As it has been already stated deterministic subspace system identification and the CUSUM algorithm are the key ingredients in the proposed algorithm, which undergoes two phases: the training and the detection phase.During the first phase (training) all the nodes in the structure report their 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 detection 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 that several variations in the input signals to excite the modes of the system are present. In the simulation part a training period of 24 h, that is, 86400 samples with sampling period 1 sec are considered. After receiving all the measurements, the central system arranges them to input-output data for each zone. For example, in the case of a zone namedA, with neighbors zones named B, C, and D and the exterior to the building zone named E, then the temperature values of zone A are considered the output of the system to be identified, whereas the temperatures of the 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. Then, a deterministic subspace identification process is run for each zone as described in the previous subsection. The relevant matrices A, B, C, and D of the state space model of each zone are communicated back to the sensor node that is responsible for monitoring the zone. Figure 2 depicts the flow diagram of the training phase. After reception of these matrices, the sensor node enters the detection phase.Figure 2
Training phase steps of the algorithm.The detection phase is split into cycles of possibly unequal lengths. In the beginning of each cycle the valueq
n, that is, the mean of the rate of the prediction error is estimated using a few samples. For the simulation part in this paper short periods of 1 min (60 samples) are used to estimate the value q
n. Then, for the rest of the cycle, the CUSUM algorithm is run to detect possible deviations from the normal operation. The length of the detection interval is quite large, in the order of 1 hour, and it should depend on the time of the day. During periods (of the day or night) with small (large) expected changes of the exterior temperature, the detection interval can be increased (decreased).Note that a positive drift of the test statistic (q
>
0) designates an abnormal increase of temperature that may be due to the onset of a fire, whereas a negative drift of the test statistic (q
<
0) designates leakage of heat that may be due to open windows during winter, or malfunctioning of the heater. The detection phase flow diagram is depicted in Figure 3.Figure 3
Detection phase steps of the algorithm.


## 4. Simulation Results

In this paper a multizone system consisting of a squared arrangement of rooms (zones) is simulated where each room (zone) is equipped with a wireless sensor node as shown in Figure4. Each node is capable to measure the temperature of the zone and convey this information to the WSN. Several protocols exist for information dissemination but their presentation is out of the scope of this paper. A suitable scheme is to let the nodes to transmit their measurements in regular periods of time with probability p and therefore the information about temperature of zones is spread in an “epidemic” fashion. In general the probability p varies from node to node. In power constrained environments this probability may depend on the available for nodes battery energy. A node with limited power supply transmits with small probability in order to elongate its life. In the proposed application it is more meaningful to adjust the transmission probability depending on the difference of successive measurements. Thus, a node that senses a noticeable difference of its measurements increases p in order to communicate faster this change to its neighbors.Figure 4
Multizone model.


### 4.1. Outside Temperature Modeling

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 [19]. This model provides the “average unit curve” given by(26)
X
a
h
=
0.463
sin
⁡
B
+
232
°
38
′
+
0.121
sin
⁡
2
B
+
55
°
21
′
+
0.031
sin
⁡
3
B
+
73
°
19
′
,where X
a
h is the unit ordinate at time h (local apparent time, LAT), B is angular measure of the time of day, h (LAT), 0
<
B
<
360
°, B
=
15
h. Using this formula along with (27)
X
a
h
=
T
o
-
T
m
T
max
⁡
-
T
min
⁡the screen temperature T
o at time h is estimated. T
max
⁡ is the maximum temperature, T
min
⁡ is the minimum temperature of the day, and T
m is the mean of the 24 hourly temperatures approximated by the daily mean temperature plus a correction factor of -
0.054
T
max
⁡
-
T
min
⁡. For example, if T
max
⁡ = 15°C, T
min
⁡ = 2°C, and T
m = 7°C, the daily variations of temperature of Figure 5 are obtained, which shows that temperature is highest at 2:00–3:00 pm.Figure 5
Empirical estimation ofT
o.


### 4.2. Single Zone Modeling

Isolating a single zone of the multizone system its dynamics can be expressed using a lumped capacity model.The assumption is made that the zone temperature distribution is uniform and therefore a single sensor is capable of providing the zone temperature values (Figure6). Another assumption is that heat transfer takes place between room’s air and the walls as well as between room’s air and the roof. North and south walls have the same effect on the room temperature and that are to be assumed the same for the east and west walls. Furthermore, the effects of ground on the room temperature are neglected. According to these assumptions, the zone model is characterized by six state variables, the room’s temperature T
r, the east-room wall temperature T
er, the west-room wall temperature T
wr, the north-room wall temperature T
nr, the south-room wall temperature T
sr
⁡, and the roof temperature T
R. Adopting the notation used in [20] the energy balance equations of the model can be written in a continuous time form as(28)
C
r
d
T
r
d
t
=
U
ew
A
ew
T
er
-
T
r
+
U
ew
A
ew
T
wr
-
T
r
+
U
ns
A
ns
T
nr
-
T
r
+
U
ns
A
ns
T
sr
⁡
-
T
r
+
U
R
A
R
T
R
-
T
r
+
q
(
t
)
C
ew
d
T
er
d
t
=
U
ew
A
ew
T
r
-
T
er
+
U
ew
A
ew
T
e
-
T
er
C
ew
d
T
wr
d
t
=
U
ew
A
ew
T
r
-
T
wr
+
U
ew
A
ew
T
w
-
T
wr
C
ns
d
T
nr
d
t
=
U
ns
A
ns
T
r
-
T
nr
+
U
ns
A
ns
T
n
-
T
nr
C
ns
d
T
sr
⁡
d
t
=
U
ns
A
ns
T
r
-
T
sr
⁡
+
U
ns
A
ns
T
s
-
T
sr
⁡
C
R
d
T
R
d
t
=
U
R
A
R
T
r
-
T
R
+
U
R
A
R
T
o
-
T
R
.Equation (28) expresses the fact that the difference between the heat transferred to the zone and the heat removed from the zone is equal to the rate change of heat in the zone. Similarly, the rest of the equations state that the heat transferred through the walls due to difference of temperatures between the outside and inside air is equal to rate change of heat through walls. In this simple model the effects of moisture, doors, windows, and pressure losses across the zone are neglected.Figure 6
Single zone model consisting of a heat source and a wireless node.The parameters used in the previous equations are as follows:A
ew: area of the east-west wall,
A
ns: area of the north-south wall,
A
R: area of the roof,
C
r: thermal capacity of the room,
C
ew: thermal capacity of the east-west wall,
C
ns: thermal capacity of the north-south wall,
C
R: thermal capacity of the roof,
U
ew: heat transfer coefficient of east-west walls,
U
ns: heat transfer coefficient of north-south walls,
U
R: heat transfer coefficient of roof,
q
(
t
): heat gain due to heating sources,
T
r: temperature of the room,
T
er: temperature of the east-room wall,
T
wr: temperature of the west-room wall,
T
nr: temperature of the north-room wall,
T
sr
⁡: temperature of the south-room wall,
T
R: temperature of the roof.The zone model parameter values are taken from [20] and are summarized as follows.A
ew: 9 m2,
A
ns: 12 m2,
A
R: 9 m2,
C
r: 47.1 KJ/°C,
C
ew: 70 KJ/°C,
C
ns: 70 KJ/°C,
C
R: 70 KJ/°C,
U
ew: 2 W/m2,
U
ns: 2 W/m2,
U
R: 1 W/m2.Equations (28) can be written compactly using the state space model(29)
d
d
t
x
t
=
A
z
x
t
+
B
z
u
(
t
)
y
t
=
C
z
x
t
,where(30)
x
t
=
T
r
T
er
T
wr
T
ns
T
sr
⁡
T
R
T
u
t
=
T
e
T
w
T
n
T
s
T
o
q
t
T
y
t
=
T
r
,where A
z, B
z are 6
×
6 matrices and C
z is an 1
×
6 matrix implicitly defined by (28).The discrete version of the model, with sampling periodT
s, is given by the equations(31)
x
k
+
1
=
A
z
d
x
k
+
B
z
d
u
k
y
k
=
C
z
x
k
,where(32)
A
z
d
=
e
A
z
T
s
B
z
d
=
A
z
-
1
A
z
d
-
I
B
z
.Note that for the multizone system of Figure2 several state variables, such as wall temperatures, are common to the individual systems. In the simulations the multizone system is treated as a single system with 42 states and nine outputs (the zones’ temperatures).


### 4.3. Training Phase

For the deterministic subsystem identification 86400 samples are used with sampling periodT
s = 1 sec. For the outside temperature Walter’s model is used with parameters given in the first row of Table 1 and a heater gain equal to 600 W. The target temperatures of the zones T
tg are set to those of Table 2 (first row) with a margin equal to T
mg = 1°C. That is, for a specific zone, the heater is on until temperature T
tg
+
T
mg is reached. After this point the heater is turned off until temperature hits the lower threshold T
tg
-
T
mg 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
max
⁡
T
min
⁡
T
m
q
(
t
)
Training
13
2
7
600
Testing
16
1
6
700Table 2
Training and testing zones’ target temperature.
Target temperatures
Zone
I
II
III
IV
V
VI
VII
VIII
IX
Training
16
18
13
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 arranged 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 the 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 the matricesA, B, C, and D for each zone as described in Section 3 are identified. During this process we must be decided on the order of the subsystems. This is achieved by looking at the singular values of the SVD decomposition in (21) and deciding on the number of dominant ones. Figure 7 depicts the singular values for zones I 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 7
Singular values of zones I and V.Next a test is run on the obtained state space models. The outer temperature parameters are set as in the second row of Table1 the heat gain equal to 700 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 5 (the initial output temperature). Figure 8 shows the evolution of the actual and predicted temperature for zone I. As it is observed, after 7000 samples (approximately 2-hour period), the state of the subsystem has converged to one that produces almost the same output as the original system. After this point the WSN nodes can enter in the detection phase.Figure 8
Real and predicted temperature of zone I.


### 4.4. Detection Phase

To demonstrate the detection capabilities of the algorithm two scenarios are used. In the first scenario, there is heat leakage, possibly due to an open window. To take into account the infiltration gain the term is added(33)
q
i
=
m
V
(
T
o
-
T
r
)
3to the right hand of (28). V is the room volume and m is an index that indicates the air changes per hour. A conservative value of m (that in the simulations is used) is 1
/
4. Assumption is made that the change in the dynamics of the system takes place at zone I, at time instant n
=
8000. Figure 9 shows the real and the predicted temperatures in this case. The target temperature is 17°C.Figure 9
Real and predicted temperature of zone I, with a heat leakage starting atn
=
8000.As it is observed the predicted values deviate from the real ones after time index 8000. In fact the predicted values do not stay in the zone[
16
°
,
18
°
], that is, the target temperature ±
1
°. This is due to the fact that in this scenario the heater is turned on more frequently and this forces the prediction model to produce higher temperature values. Crossing the upper limit of the operation zone is an indicator of a heat loss and this information could be fused with the results of the CUSUM algorithm for more reliable detection of the abnormal condition. Figure 10 shows the drift of the test statistic of the CUSUM algorithm in case of normal and abnormal operation.Figure 10
Drift of the test statistic for scenario I, under no change and a heat leakage starting atn
=
8000.The CUSUM algorithm started at time 7000 and 60 samples were used to estimateq
n that is the mean of the rate of the prediction error. In the second scenario there is an extra source of heat with power 100 W. This source is activated at time instant 8000 and it remains on independently of the zone’s heater. Figure 11 depicts the real and the predicted temperatures in this case.Figure 11
Real and predicted temperature of zone I, with an extra heat source powered on atn
=
8000.As it is observed, the predicted values in this case fall below the zone[
16
°
,
18
°
]. This is because the heater of the zone is turned on less often since there is an additional heat power in the zone. However, the prediction model is unaware of this fact and “sees” only the power on-off process of the zone’s heater. Similar to the previous scenario, crossing the lower limit of the operation zone is an indicator of an extra heat source in the zone. The CUSUM test results for this scenario are depicted in Figure 12. In both scenarios, the drift test statistic crosses the threshold of T
f
=
5.42 within 300 samples and therefore with a false alarm rate equal to f
=
0.001 it can be determined that a change has occurred or not within 5 min.Figure 12
Drift of the test statistic for scenario II, under no change and an extra heat source powered on atn
=
8000.It should be noticed 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 for example) a low dimensionality subsystem of order 2 can capture its behaviour.For more complex systems the order of the identified subsystem can be chosen high enough. For nonlinear systems the identification of a time varying system is possible using a recursive update of the model.Next, the SID method is tested when an exogenous heat noise is present. This noise is modeled as(34)
n
t
=
∑
i
a
i
w
t
-
T
i
T
i
,where a
i is uniformly distributed in the range [
50,100
] Watt and T
i is exponentially distributed with mean 180 sec. w
(
t
) denotes a rectangular window of unit length and height equal to one. A sample of this process is given in Figure 13.Figure 13
Exogenous heat noise profile.The real and predicted temperature values for this case are shown in Figure14. As it is observed the predicted values are always lower than the real ones and this is due to the fact that the SID method ignores the extra source of heat and it bases its predictions solely on the zone’s heater on/off periods. This is more clear after time instant 8000 where an extra heat source of 100 Watt is always present. Several techniques can be used in this case to detect the discrepancies from the target operating restrictions. Deviations of the mean predicted values from the anticipated ones or the number of crossings of a low threshold are some possible solutions.Figure 14
Real and predicted temperature of zone I, with an exogenous heat noise depicted in Figure13 and an extra heat source powered on at n
=
8000.


## 4.1. Outside Temperature Modeling

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 [19]. This model provides the “average unit curve” given by(26)
X
a
h
=
0.463
sin
⁡
B
+
232
°
38
′
+
0.121
sin
⁡
2
B
+
55
°
21
′
+
0.031
sin
⁡
3
B
+
73
°
19
′
,where X
a
h is the unit ordinate at time h (local apparent time, LAT), B is angular measure of the time of day, h (LAT), 0
<
B
<
360
°, B
=
15
h. Using this formula along with (27)
X
a
h
=
T
o
-
T
m
T
max
⁡
-
T
min
⁡the screen temperature T
o at time h is estimated. T
max
⁡ is the maximum temperature, T
min
⁡ is the minimum temperature of the day, and T
m is the mean of the 24 hourly temperatures approximated by the daily mean temperature plus a correction factor of -
0.054
T
max
⁡
-
T
min
⁡. For example, if T
max
⁡ = 15°C, T
min
⁡ = 2°C, and T
m = 7°C, the daily variations of temperature of Figure 5 are obtained, which shows that temperature is highest at 2:00–3:00 pm.Figure 5
Empirical estimation ofT
o.


## 4.2. Single Zone Modeling

Isolating a single zone of the multizone system its dynamics can be expressed using a lumped capacity model.The assumption is made that the zone temperature distribution is uniform and therefore a single sensor is capable of providing the zone temperature values (Figure6). Another assumption is that heat transfer takes place between room’s air and the walls as well as between room’s air and the roof. North and south walls have the same effect on the room temperature and that are to be assumed the same for the east and west walls. Furthermore, the effects of ground on the room temperature are neglected. According to these assumptions, the zone model is characterized by six state variables, the room’s temperature T
r, the east-room wall temperature T
er, the west-room wall temperature T
wr, the north-room wall temperature T
nr, the south-room wall temperature T
sr
⁡, and the roof temperature T
R. Adopting the notation used in [20] the energy balance equations of the model can be written in a continuous time form as(28)
C
r
d
T
r
d
t
=
U
ew
A
ew
T
er
-
T
r
+
U
ew
A
ew
T
wr
-
T
r
+
U
ns
A
ns
T
nr
-
T
r
+
U
ns
A
ns
T
sr
⁡
-
T
r
+
U
R
A
R
T
R
-
T
r
+
q
(
t
)
C
ew
d
T
er
d
t
=
U
ew
A
ew
T
r
-
T
er
+
U
ew
A
ew
T
e
-
T
er
C
ew
d
T
wr
d
t
=
U
ew
A
ew
T
r
-
T
wr
+
U
ew
A
ew
T
w
-
T
wr
C
ns
d
T
nr
d
t
=
U
ns
A
ns
T
r
-
T
nr
+
U
ns
A
ns
T
n
-
T
nr
C
ns
d
T
sr
⁡
d
t
=
U
ns
A
ns
T
r
-
T
sr
⁡
+
U
ns
A
ns
T
s
-
T
sr
⁡
C
R
d
T
R
d
t
=
U
R
A
R
T
r
-
T
R
+
U
R
A
R
T
o
-
T
R
.Equation (28) expresses the fact that the difference between the heat transferred to the zone and the heat removed from the zone is equal to the rate change of heat in the zone. Similarly, the rest of the equations state that the heat transferred through the walls due to difference of temperatures between the outside and inside air is equal to rate change of heat through walls. In this simple model the effects of moisture, doors, windows, and pressure losses across the zone are neglected.Figure 6
Single zone model consisting of a heat source and a wireless node.The parameters used in the previous equations are as follows:A
ew: area of the east-west wall,
A
ns: area of the north-south wall,
A
R: area of the roof,
C
r: thermal capacity of the room,
C
ew: thermal capacity of the east-west wall,
C
ns: thermal capacity of the north-south wall,
C
R: thermal capacity of the roof,
U
ew: heat transfer coefficient of east-west walls,
U
ns: heat transfer coefficient of north-south walls,
U
R: heat transfer coefficient of roof,
q
(
t
): heat gain due to heating sources,
T
r: temperature of the room,
T
er: temperature of the east-room wall,
T
wr: temperature of the west-room wall,
T
nr: temperature of the north-room wall,
T
sr
⁡: temperature of the south-room wall,
T
R: temperature of the roof.The zone model parameter values are taken from [20] and are summarized as follows.A
ew: 9 m2,
A
ns: 12 m2,
A
R: 9 m2,
C
r: 47.1 KJ/°C,
C
ew: 70 KJ/°C,
C
ns: 70 KJ/°C,
C
R: 70 KJ/°C,
U
ew: 2 W/m2,
U
ns: 2 W/m2,
U
R: 1 W/m2.Equations (28) can be written compactly using the state space model(29)
d
d
t
x
t
=
A
z
x
t
+
B
z
u
(
t
)
y
t
=
C
z
x
t
,where(30)
x
t
=
T
r
T
er
T
wr
T
ns
T
sr
⁡
T
R
T
u
t
=
T
e
T
w
T
n
T
s
T
o
q
t
T
y
t
=
T
r
,where A
z, B
z are 6
×
6 matrices and C
z is an 1
×
6 matrix implicitly defined by (28).The discrete version of the model, with sampling periodT
s, is given by the equations(31)
x
k
+
1
=
A
z
d
x
k
+
B
z
d
u
k
y
k
=
C
z
x
k
,where(32)
A
z
d
=
e
A
z
T
s
B
z
d
=
A
z
-
1
A
z
d
-
I
B
z
.Note that for the multizone system of Figure2 several state variables, such as wall temperatures, are common to the individual systems. In the simulations the multizone system is treated as a single system with 42 states and nine outputs (the zones’ temperatures).


## 4.3. Training Phase

For the deterministic subsystem identification 86400 samples are used with sampling periodT
s = 1 sec. For the outside temperature Walter’s model is used with parameters given in the first row of Table 1 and a heater gain equal to 600 W. The target temperatures of the zones T
tg are set to those of Table 2 (first row) with a margin equal to T
mg = 1°C. That is, for a specific zone, the heater is on until temperature T
tg
+
T
mg is reached. After this point the heater is turned off until temperature hits the lower threshold T
tg
-
T
mg 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
max
⁡
T
min
⁡
T
m
q
(
t
)
Training
13
2
7
600
Testing
16
1
6
700Table 2
Training and testing zones’ target temperature.
Target temperatures
Zone
I
II
III
IV
V
VI
VII
VIII
IX
Training
16
18
13
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 arranged 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 the 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 the matricesA, B, C, and D for each zone as described in Section 3 are identified. During this process we must be decided on the order of the subsystems. This is achieved by looking at the singular values of the SVD decomposition in (21) and deciding on the number of dominant ones. Figure 7 depicts the singular values for zones I 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 7
Singular values of zones I and V.Next a test is run on the obtained state space models. The outer temperature parameters are set as in the second row of Table1 the heat gain equal to 700 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 5 (the initial output temperature). Figure 8 shows the evolution of the actual and predicted temperature for zone I. As it is observed, after 7000 samples (approximately 2-hour period), the state of the subsystem has converged to one that produces almost the same output as the original system. After this point the WSN nodes can enter in the detection phase.Figure 8
Real and predicted temperature of zone I.


## 4.4. Detection Phase

To demonstrate the detection capabilities of the algorithm two scenarios are used. In the first scenario, there is heat leakage, possibly due to an open window. To take into account the infiltration gain the term is added(33)
q
i
=
m
V
(
T
o
-
T
r
)
3to the right hand of (28). V is the room volume and m is an index that indicates the air changes per hour. A conservative value of m (that in the simulations is used) is 1
/
4. Assumption is made that the change in the dynamics of the system takes place at zone I, at time instant n
=
8000. Figure 9 shows the real and the predicted temperatures in this case. The target temperature is 17°C.Figure 9
Real and predicted temperature of zone I, with a heat leakage starting atn
=
8000.As it is observed the predicted values deviate from the real ones after time index 8000. In fact the predicted values do not stay in the zone[
16
°
,
18
°
], that is, the target temperature ±
1
°. This is due to the fact that in this scenario the heater is turned on more frequently and this forces the prediction model to produce higher temperature values. Crossing the upper limit of the operation zone is an indicator of a heat loss and this information could be fused with the results of the CUSUM algorithm for more reliable detection of the abnormal condition. Figure 10 shows the drift of the test statistic of the CUSUM algorithm in case of normal and abnormal operation.Figure 10
Drift of the test statistic for scenario I, under no change and a heat leakage starting atn
=
8000.The CUSUM algorithm started at time 7000 and 60 samples were used to estimateq
n that is the mean of the rate of the prediction error. In the second scenario there is an extra source of heat with power 100 W. This source is activated at time instant 8000 and it remains on independently of the zone’s heater. Figure 11 depicts the real and the predicted temperatures in this case.Figure 11
Real and predicted temperature of zone I, with an extra heat source powered on atn
=
8000.As it is observed, the predicted values in this case fall below the zone[
16
°
,
18
°
]. This is because the heater of the zone is turned on less often since there is an additional heat power in the zone. However, the prediction model is unaware of this fact and “sees” only the power on-off process of the zone’s heater. Similar to the previous scenario, crossing the lower limit of the operation zone is an indicator of an extra heat source in the zone. The CUSUM test results for this scenario are depicted in Figure 12. In both scenarios, the drift test statistic crosses the threshold of T
f
=
5.42 within 300 samples and therefore with a false alarm rate equal to f
=
0.001 it can be determined that a change has occurred or not within 5 min.Figure 12
Drift of the test statistic for scenario II, under no change and an extra heat source powered on atn
=
8000.It should be noticed 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 for example) a low dimensionality subsystem of order 2 can capture its behaviour.For more complex systems the order of the identified subsystem can be chosen high enough. For nonlinear systems the identification of a time varying system is possible using a recursive update of the model.Next, the SID method is tested when an exogenous heat noise is present. This noise is modeled as(34)
n
t
=
∑
i
a
i
w
t
-
T
i
T
i
,where a
i is uniformly distributed in the range [
50,100
] Watt and T
i is exponentially distributed with mean 180 sec. w
(
t
) denotes a rectangular window of unit length and height equal to one. A sample of this process is given in Figure 13.Figure 13
Exogenous heat noise profile.The real and predicted temperature values for this case are shown in Figure14. As it is observed the predicted values are always lower than the real ones and this is due to the fact that the SID method ignores the extra source of heat and it bases its predictions solely on the zone’s heater on/off periods. This is more clear after time instant 8000 where an extra heat source of 100 Watt is always present. Several techniques can be used in this case to detect the discrepancies from the target operating restrictions. Deviations of the mean predicted values from the anticipated ones or the number of crossings of a low threshold are some possible solutions.Figure 14
Real and predicted temperature of zone I, with an exogenous heat noise depicted in Figure13 and an extra heat source powered on at n
=
8000.


## 5. Conclusions and Future Work

A method for the detection of abnormal behaviour in HVAC systems based on the multizone principle is proposed. The method uses deterministic subspace identification to obtain state space system models that will provide the reference temperatures for each zone. As it is shown using simulations, simple state space models of order 2 are able to capture the dynamics of each zone and to produce predicted temperature profiles that closely follow real measurements. The CUSUM algorithm is used next to detect possible divergences of the rate of change of the difference between the real and the reference values. The method was tested using two different scenarios, that is, heat leakage due to infiltration losses and heat gains due to an exogenous source, that is, the onset of a fire. In both cases, the abnormal operation is detected within 5 min with a small false alarm rate equal to 0.001.Two directions are foreseen for future work: in the first direction, investigating the use of sophisticated classifiers such as SVM or Neural Nets incorporating more features of the reference signals in the decision process for the early detection of abnormal behaviour and in the second direction, focusing on utilizing heterogeneous sensors (temperature, occupancy) in the WSN aiming at optimizing the energy consumption at each zone while in parallel maintaining the comfort conditions of the occupants. To this end, combining the powerful subspace identification process should be combined with results from Optimal Stopping Theory.



---

*Source: 693749-2015-03-11.xml*