How to detect start and finish of temperature control in temperature time series

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. Smaller subset to better show the behavior of the data: Signal Processing Asked by Regnav on January 1, 2021

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.

Correct answer by Hilmar on January 1, 2021

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], [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. 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.

Answered by Peter K. on January 1, 2021

Related Questions

Frequency conversion in IIR lowpass-to-lowpass transformation

0  Asked on October 24, 2021 by rohitm

How to estimate the modulation transfer function of images?

1  Asked on October 24, 2021

Any Ideas for Extracting the lines on a Football Pitch?

0  Asked on October 24, 2021 by jonah-f

Output of marginally stable systems

2  Asked on October 24, 2021

Converting a lowpass filter to a highpass filter. FIR filter-type 1

1  Asked on March 10, 2021 by bl-lov

What are the practical constraints on designing Sensing matrix in compressed Sensing?

2  Asked on March 7, 2021 by digi1

High pass filter not attenuating signal at start of input

1  Asked on February 24, 2021 by ankit-chudasama

Why do I have a high peak at the beginning of my FFT (not DC)?

1  Asked on February 24, 2021 by gabriel-galeote-checa

Response time of residual current device and sampling times

0  Asked on February 23, 2021 by blue_electronx

What is the difference between a range image and a depth map?

1  Asked on February 21, 2021 by uxkqez7

Weird ringing at signal start and end for default Matlab lowpass filter

1  Asked on February 19, 2021 by anand-kulkarni

1  Asked on February 12, 2021 by jumbot

Most efficient way to find single dominant frequency (without amplitude) in analog signal

0  Asked on February 11, 2021 by tobalt

Blind Estimation of Signal Parameter and Noise Variance

1  Asked on February 10, 2021 by dsp-guy-sam

Finding the system output by convolution

1  Asked on February 10, 2021 by jisbon

Estimating a Filter’s Envelope

0  Asked on February 8, 2021

Audio Activity Detection and Classification

1  Asked on February 7, 2021

How can i identify number of return of the ping in the sonar?

0  Asked on February 6, 2021 by jay-patel

Difference between Analog IQ sampling and creating IQ digitally

1  Asked on February 5, 2021 by malik12