3
$\begingroup$

I have a huge dataset containing temperature data inside a building. I want to extract the time that the building starts and stops controlling the temperature (approximately around the vertical black lines in the pictures). In this case the temperature control starts and stops at the same time. However, I want to detect the the same for a building where the start and stopping time is variable, per day and per area inside the building. How would you go about finding a solutions for such a problem?

I can quite reliably detect the first line by detecting when the absolute velocity / acceleration reaches above a certain threshold, but this method seems less reliable for detecting the stopping point.

enter image description here

Smaller subset to better show the behavior of the data: enter image description here

$\endgroup$
  • 3
    $\begingroup$ Good question! I don't have time for an answer just yet (I'll try later today), but you might want to look at this answer and the references in it. The CUSUM algorithm is the generic answer for any change detection problem. The issue will be how to quantify it for your specific signal case. $\endgroup$ – Peter K. Sep 9 at 14:43
  • 1
    $\begingroup$ Thanks a lot, I will take a look there. Looking forward to your answer! $\endgroup$ – Regnav Sep 9 at 14:47
2
$\begingroup$

One option: It's pretty clear that in the uncontrolled range the temperature follows a simple exponential heating/cooling curve. Something like $$T_{n+1} = T_n + (T_{Ambient} - T_n) \cdot e^{-\Delta t / t_{Building}} $$

where $\Delta t$ is the time difference between samples and $ t_{Building}$ the thermal time constant of your building (which may be slightly different at different sensor locations).

You could simply test the difference of each new sample to the model prediction and set an threshold of the allowable error. The only two model parameters are ambient temperature and thermal time constant of the building.

The thermal time constant is a function of the thermal resistance of the building and it really doesn't change unless someone opens/closes windows or there is some actual structural work happening. Looking at your data that should be easy to estimate. Ambient you can try to get directly either with a sensor or looking at a reliable weather report. You can also just use the data from the first few sensor to estimate ambient and then test whether the other sensors agree or not.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you, I implemented this and it is working great! I am still struggling a little bit to find a good measure / threshold when the measurements are too far off from the estimated e-curve, but I will figure that out along the way. Thanks a lot! $\endgroup$ – Regnav Sep 10 at 15:45
0
$\begingroup$

So this is an attempt to do it, but it's not really based on how your "controlled" section is.

The idea is to model the beginning and ending sections as a first order system using a Kalman Filter (as Hilmar's good answer suggests). The controlled section will be different from this.

You can look at the innovations (error term) in the Kalman Filter to see whether the innovations is noise, or something more structured.

The first code is just generating the signal.

from numpy import log10, asarray, polyfit, ceil, arange, exp, sin, pi, log, random, sum, diff
import matplotlib.pyplot as plt



T = 1000
Ton = 300
Toff = 650


#
# First period: temperature rising or falling as a first order system.
#
# IC @ 1 = min FC @ Ton = max
# f(t) = K1 + K2 exp(-t/tau)
# f(1) = K1 + K2 exp(-1/tau) = min        (1)
# f(Ton) = K1 + K2 exp(-Ton/tau) = max    (2)
#
# (1) - (2) --> K2 ( exp(-1/tau) - exp(-Ton/tau) ) = min - max --> K2 = (min - max) / (exp(-1/tau) - exp(-Ton/tau) )
mx = 100
mn = 10
tau = 150
time_period_1 = list(arange(1,Ton))
K2 = (mn - mx) / (exp(-1/tau) - exp(-Ton/tau))
print(K2)
K1 = mn - K2*exp(-1/tau)
K1_2 = mx - K2*exp(-Ton/tau)
print(str(K1) + " " + str(K1_2))

temperature = [K1 + K2*exp(-x/tau) + random.normal(0,0.001) for x in time_period_1]

plt.figure(1)
plt.plot(time_period_1, temperature)

#
# Second period: being controlled.
#
time_period_2 = list(arange(Ton, Toff))

variation = 50
mean_value = temperature[Ton-2]
tau2 = 120

temperature2 = [variation*sin(2.0*pi*(x/100))*exp(-(x-Ton)/tau2) + mean_value for x in time_period_2]
plt.plot(time_period_2, temperature2)

#
# Third period: back to first order.
#
# IC @ Toff = last value of previous period FC @ T = mx3
# f(t) = K1 + K2 exp(-t/tau)
# f(Toff) = K1 + K2 exp(-Toff/tau) = min        (1)
# f(T) = K1 + K2 exp(-T/tau) = max              (2)
#
# (1) - (2) --> K2 ( exp(-Toff/tau) - exp(-T/tau) ) = last value - mx3
#           --> K2 = (last value - mx) / (exp(-Toff/tau) - exp(-T/tau) )
mx3 = 110
mn3 = temperature2[Toff-Ton-2]
tau2 = 50
time_period_3 = list(arange(Toff, T))

K23 = (mn3 - mx3) / (exp(-Toff/tau2) - exp(-T/tau2))
print(K23)
K13 = mn3 - K23*exp(-Toff/tau2)
K13_2 = mx3 - K23*exp(-T/tau2)
print(str(K13) + " " + str(K13_2))

temperature3 = [K13 + K23*exp(-x/tau2) for x in time_period_3]
plt.plot(time_period_3, temperature3)

all_temps = list(temperature) + list(temperature2) + list(temperature3)

plt.figure(2)
plt.plot(arange(1,T), all_temps)

And then setting up the Kalman filter:

import matplotlib.pyplot as plt
import numpy as np
from filterpy.kalman import KalmanFilter
from filterpy.common import Q_discrete_white_noise, Saver

dt = 0.1
r_std = 0.1
q_std = 0.1

cv = KalmanFilter(dim_x=2, dim_z=1)
cv.x = np.array([[all_temps[0]], [10.]]) # position, velocity
cv.F = np.array([[1, dt],[0, 1]])
cv.R = np.array([[r_std**2]])
cv.H = np.array([[1., 0.]])
cv.P = np.diag([.1**2, .03**2])
cv.Q = Q_discrete_white_noise(2, dt, q_std**2)

saver = Saver(cv)
for z in range(len(all_temps)):
    cv.predict()
    cv.update([all_temps[z] + random.randn()*q_std ])
    saver.save() # save the filter's state

saver.to_array()
plt.figure(figsize=(10,10))
plt.plot(saver.x[:, 0], 'b.')
plt.plot(saver.x[:, 1], 'go')
plt.plot(all_temps,'k.')

# plot all of the priors
plt.plot(saver.x_prior[:, 0], 'r+')

# plot mahalanobis distance
plt.figure()
plt.figure(figsize=(10,10))
plt.plot(saver.P[:,0,0])
plt.plot(saver.P[:,0,1])
plt.plot(saver.P[:,1,0])
plt.plot(saver.P[:,1,1])

plt.figure()
plt.figure(figsize=(10,10))
plt.plot(abs(saver.y[:,0,0]))

N = 50
smoothed_innovations = np.convolve(abs(saver.y[:,0,0]), np.ones((N,))/N, mode='valid')
plt.plot(smoothed_innovations)
threshold = np.mean(smoothed_innovations[100:200])
standard_deviation = np.std(smoothed_innovations[100:200])

plt.plot(8*(smoothed_innovations > threshold + 3*standard_deviation))

plt.savefig('Q70221.png')

The result is shown below.

Attempt at thresholding

The blue line is the absolute value of the innovations. The orange line is a smoothed version of it. The green line indicates when the orange line is above or below the selected threshold.

Not really CUSUM, but I'll work on making it closer.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.