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.

Signal Processing Asked by Regnav on January 1, 2021

2 AnswersOne 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[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.

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

1 Asked on October 24, 2021 by mohammadsadeq-borjiyan

digital filters filter design infinite impulse response lowpass filter matlab

0 Asked on October 24, 2021 by rohitm

1 Asked on October 24, 2021

0 Asked on October 24, 2021 by jonah-f

2 Asked on October 24, 2021

1 Asked on March 10, 2021 by bl-lov

discrete signals filter design filters finite impulse response lowpass filter

2 Asked on March 7, 2021 by digi1

1 Asked on February 24, 2021 by ankit-chudasama

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

0 Asked on February 23, 2021 by blue_electronx

1 Asked on February 21, 2021 by uxkqez7

1 Asked on February 19, 2021 by anand-kulkarni

0 Asked on February 11, 2021 by tobalt

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

1 Asked on February 10, 2021 by jisbon

convolution differential equation laplace transform linear systems math

1 Asked on February 7, 2021

0 Asked on February 6, 2021 by jay-patel

1 Asked on February 5, 2021 by malik12

Get help from others!

Recent Answers

- Philipp on How do i draw a ray in unity
- DMGregory on MouseLook Script “Pops” back to the last value when the script is enabled after being disabled or destroyed
- eric_kernfeld on How to test consistency of responses?
- Justin Markwell on Unity app crashes when using unmodified custom Android manifest (didn’t find class “UnityPlayerActivity”)
- kjetil b halvorsen on How to test consistency of responses?

Recent Questions

- MouseLook Script “Pops” back to the last value when the script is enabled after being disabled or destroyed
- Unity app crashes when using unmodified custom Android manifest (didn’t find class “UnityPlayerActivity”)
- How do i draw a ray in unity
- How to test consistency of responses?
- How can I understand these variograms?

© 2022 AnswerBun.com. All rights reserved.