MATLAB Examples

Continuous monitoring of wireless network of temperature sensors using MATLAB® and XBee® (Part 4)

Contents

Introduction

In part one, part two, and part three of this series of blog post I went through the process of reading voltages from a network of XBee® modules, building and testing my wireless network of temperature sensors, calibrating the sensors, and gathering temperature data from the network of sensors in my apartment. In this post I will show you some of the basic analyses I performed on the data. In particular, I have two goals:

  1. determine how long it takes for my apartment to warm up after the heat has been turned on, and
  2. find out what rooms in my apartment take the longest to heat up, so I can make the necessary adjustments to my radiators so that all the rooms warm up at the same rate.

Overview of Data

As I mentioned in a previous post, I collected the temperature every two minutes over the course of 9 days. I placed 14 sensors in my apartment: 9 located inside, 2 located outside, and 3 located in radiators. The data is stored in the file twoweekstemplog.txt.

[tempF,ts,location,lineSpecs] = XBeeReadLog('twoweekstemplog.txt',60);
tempF = calibrateTemperatures(tempF);
plotTemps(ts,tempF,location,lineSpecs)
legend('off')
xlabel('Date')
title('All Data Overview')

Figure 1: All temperature data from a 9 day period.

That graph is a bit too cluttered to be very meaningful. Let me remove the radiator data and the legend and see if that helps.

notradiator = [1 2 3 5 6 7 8 9 10 12 13];
plotTemps(ts,tempF(:,notradiator),location(notradiator),lineSpecs(notradiator,:))
legend('off')
xlabel('Date')
title('All Inside and Outside Temperature Data')

Figure 2: Just inside and outside temperature data, with radiator data removed.

Now I can see some places where one of the outdoor temperature sensors (blue line) gave erroneous data, so let's remove those data points. This data was collected in March in Massachusetts, so I can safely assume the outdoor temperature never reached 80 F. I replaced any values above 80 F with NaN (not-a-number) so they are ignored in further analysis.

outside = [3 10];
outsideTemps = tempF(:,outside);
toohot = outsideTemps>80;
outsideTemps(toohot) = NaN;
tempF(:,outside) = outsideTemps;

plotTemps(ts,tempF(:,notradiator),location(notradiator),lineSpecs(notradiator,:))
legend('off')
xlabel('Date')
title('Cleaned-up Inside and Outside Temperature Data')

Figure 3: Inside and outside temperature data with erroneous data removed.

I'll also remove all but one inside sensor per room, and give the remaining sensors shorter names, to keep the graph from getting too cluttered.

show = [1 5 9 12 10 3];
location(show) = ...
    {'Bedroom','Kitchen','Living Room','Office','Front Porch','Side Yard'}';

plotTemps(ts,tempF(:,show),location(show),lineSpecs(show,:))
ylim([0 90])
legend('Location','SouthEast')
xlabel('Date')
title('Summary of Temperature Data')

Figure 4: Summary of temperature data with only one inside temperature sensor per room with outside temperatures.

That looks much better. This data was collected over the course of 9 days, and the first thing that stands out to me is the periodic outdoor temperature, which peaks every day at around noon. I also notice a sharp spike in the side yard (green) temperature on most days. My front porch (blue) is located on the north side of my apartment, and does not get much sun. My side yard is on the east side of my apartment, and that spike probably corresponds to when the sun hits the sensor from between my apartment and the building next door.

When do my radiators start to heat up?

The radiator temperature can be used to measure how long it takes for my boiler and radiators to warm up after the heat has been turned on. Let's take a look at 1 day of data from the living room radiator:

% Grab the Living Room Radiator Temperature (index 11) from the |tempF| matrix.
radiatorTemp = tempF(:,11);

% Fill in any missing values:
validts = ts(~isnan(radiatorTemp));
validtemp = radiatorTemp(~isnan(radiatorTemp));
nants = ts(isnan(radiatorTemp));
radiatorTemp(isnan(radiatorTemp)) = interp1(validts,validtemp,nants);

% Plot the data
oneday = [ts(1) ts(1)+1];
figure
plot(ts,radiatorTemp,'k.-')
xlim(oneday)
xlabel('Time')
ylabel('Radiator Temperature (\circF)')
title('Living Room Radiator Temperature')
datetick('keeplimits')
snapnow

Figure 5: One day of temperature data from the living room radiator.

As expected, I see a sharp rise in the radiator temperature, followed by a short leveling off (when the radiator temperature maxes out the temperature sensor), and finally a gradual cooling of the radiator. Let me superimpose the rate of change in temperature onto the plot.

tempChange = diff([NaN; radiatorTemp]);

hold on
plot(ts,tempChange,'b.-')
legend({'Temperature', 'Temperature Change'},'Location','Best')

Figure 6: One day of data from the living room radiator with temperature change.

It looks like I can detect those peaks by looking for large jumps in the temperature. After some trial and error, I settled on three criteria to identify when the heat comes on:

  1. Change in temperature greater than four times the previous change in temperature.
  2. Change in temperature of more than 1 degree F.
  3. Keep the first in a sequence of matching points (remove doubles)
fourtimes = [tempChange(2:end)>abs(4*tempChange(1:end-1)); false];
greaterthanone = [tempChange(2:end)>1; false];
heaton = fourtimes & greaterthanone;
doubles = [false; heaton(2:end) & heaton(1:end-1)];
heaton(doubles) = false;

Let's see how well I detected those peaks by superimposing red dots over the times I detected.

figure
plot(ts,radiatorTemp,'k.-')
hold on
plot(ts(heaton),radiatorTemp(heaton),'r.','MarkerSize',20)
xlim(oneday);
datetick('keeplimits')
xlabel('Time')
ylabel('Radiator Temperature (\circF)')
title('Heat On Event Detection')
legend({'Temperature', 'Heat On Event'},'Location','Best')

Figure 7: Radiator temperature with heating events marked with red dots.

Looks pretty good, which means now I have a list of all the times that the heat came on in my apartment.

heatontimes = ts(heaton);

How long does it take for my heat to turn on?

I currently have a programmable 5/2 thermostat, which means I can set one program for weekdays (Monday through Friday) and one program for both Saturday and Sunday. I know my thermostat is set to go down to 62 at night, and back up to 68 at 6:15am Monday through Friday and 10:00am on Saturday and Sunday. I used that knowledge to determine how long after my thermostat activates that my radiators warm up.

I started by creating a vector of all the days in the test period. I removed Monday because I manually turned on the thermostat early that day.

mornings = floor(min(ts)):floor(max(ts));
mornings(2) = []; % Remove Monday

Then I added either 6:15am or 10:00am to each day depending on whether it was a weekday or a weekend.

isweekend = weekday(mornings) == 1 | weekday(mornings) == 7;
mornings(isweekend) = mornings(isweekend)+10/24; % 10:00 AM
mornings(~isweekend) = mornings(~isweekend)+6.25/24; % 6:15 AM

Next I looked for the first time the heat came on after the programmed time each morning.

heatontimes_mat = repmat(heatontimes,1,length(mornings));
mornings_mat = repmat(mornings,length(heatontimes),1);
timelag = heatontimes_mat - mornings_mat;
timelag(timelag<=0) = NaN;
[delay, heatind] = min(timelag);
delay = round(delay*24*60);

Let's take a look at those times to make sure we found the right ones. In this plot, I'll circle in blue the first time the heat comes on each morning, and plot blue vertical lines indicating when the thermostat turns on each morning.

heatontemp = radiatorTemp(heaton);
onemorning = mornings(3)+[-1/24 5/24];

figure
plot(ts,radiatorTemp,'k.-')
hold on
plot(heatontimes,heatontemp,'r.','MarkerSize',20)
plot(heatontimes(heatind),heatontemp(heatind),'bo','MarkerSize',10)
plot([mornings;mornings],repmat(ylim',1,length(mornings)),'b-');
xlim(onemorning);
datetick('keeplimits')
xlabel('Time')
ylabel('Radiator Temperature (\circF)')
title('Detection of Scheduled Heat On Events')
legend({'Temperature', 'Heat On Event', 'Scheduled Heat On Event',...
    'Scheduled Event'},'Location','Best')

Figure 8: Six hours of radiator data, with a blue line indicating when the thermostat turned on in the morning, and blue circle indicating the corresponding heat on event of the radiator.

Let's look at a histogram of those delays:

figure
hist(delay,min(delay):max(delay))
xlabel('Minutes')
ylabel('Frequency')
title('How long before the radiator starts to warm up?')

Figure 9: Histogram showing delay between thermostat activation and the radiators starting to warm up.

It looks like the delay between the thermostat coming on in the morning and the radiators starting to warming up can range from 7 minutes to as high as 24 minutes, but on average this delay is around 12-13 minutes.

heatondelay = 12;

How long does it take for the radiators to warm up?

Once the radiators start to warm up, it takes a few minutes for them to reach full temperature. Let's look at how long this takes. I'll look for times when the radiator temperature first maxes out the temperature sensor after having been below the maximum for at least 10 minutes (5 samples).

maxtemp = max(radiatorTemp);
radiatorhot = radiatorTemp(6:end)==maxtemp & ...
    radiatorTemp(1:end-5)<maxtemp &...
    radiatorTemp(2:end-4)<maxtemp &...
    radiatorTemp(3:end-3)<maxtemp &...
    radiatorTemp(4:end-2)<maxtemp &...
    radiatorTemp(5:end-1)<maxtemp;
radiatorhot = [false(5,1); radiatorhot];
radiatorhottimes = ts(radiatorhot);

Let's see how well that worked:

figure
plot(ts,radiatorTemp,'k.-')
hold on
plot(ts(heaton),radiatorTemp(heaton),'r.','MarkerSize',20)
plot(ts(radiatorhot),radiatorTemp(radiatorhot),'b.','MarkerSize',20)
xlim(onemorning);
datetick('keeplimits')
xlabel('Time')
ylabel('Radiator Temperature (\circF)')
title('Radiator Hot Event Detection')
legend({'Temperature', 'Heat On Event', 'Radiator Hot Event'},...
    'Location','Best')

Figure 10: Six hours of radiator data, with red dots indicating the heat coming on and blue dots indicating the radiator is hot.

Now I'll match the radiatorhottimes to the heatontimes using the same technique I used above.

radiatorhottimes_mat = repmat(radiatorhottimes',length(heatontimes),1);
heatontimes_mat = repmat(heatontimes,1,length(radiatorhottimes));
timelag = radiatorhottimes_mat - heatontimes_mat;
timelag(timelag<=0) = NaN;
[delay, foundmatch] = min(timelag);
delay = round(delay*24*60);

Let's look at a histogram of those delays:

figure
hist(delay,min(delay):2:max(delay))
xlabel('Minutes');
ylabel('Frequency')
title('How long does the radiator take to warm up?')

Figure 11: Histogram showing time required for the radiators to warm up.

It looks like the radiators take between 4 and 8 minutes from when they start to warm up until they are at full temperature.

radiatorheatdelay = 6;

Later on in my analysis, I will only want to use times that the heat came on and the radiators reached full temperature, so I will only keep those values in the heaton vector.

heatonind = find(heaton);
heatonind = heatonind(foundmatch);

heaton = false(size(heaton));
heaton(heatonind) = true;

When do the radiators cool off?

I'll use the same technique as above to detect when the radiators start to cool off, which must mean the heat has gone off.

heatoff = radiatorTemp(1:end-5)==maxtemp & ...
    radiatorTemp(2:end-4)<maxtemp &...
    radiatorTemp(3:end-3)<maxtemp &...
    radiatorTemp(4:end-2)<maxtemp &...
    radiatorTemp(5:end-1)<maxtemp &...
    radiatorTemp(6:end)<maxtemp;
heatoff = [heatoff; false(5,1)];
heatofftimes = ts(heatoff);
heatoffind = find(heatoff);

Let's take a look at the heat on, radiator hot, and heat off times all together in one graph.

figure
plot(ts,radiatorTemp,'k.-')
hold on
plot(ts(heaton),radiatorTemp(heaton),'r.','MarkerSize',20)
plot(ts(radiatorhot),radiatorTemp(radiatorhot),'b.','MarkerSize',20)
plot(ts(heatoff),radiatorTemp(heatoff),'g.','MarkerSize',20)
xlim(onemorning);
datetick('keeplimits')
xlabel('Time')
ylabel('Radiator Temperature (\circF)')
title('Radiator Cooling Event Detection')
legend({'Temperature', 'Heat On Event', 'Radiator Hot Event',...
    'Radiator Cooling Event'},'Location','Best')

Figure 12: Six hours of radiator data, with red dots indicating the heat coming on, blue dots indicating the radiator is hot, and green dots indicating the radiator starting to cool.

How long does it take my living room to warm up?

Now that I know it takes about 12-13 minutes from the time my thermostat activates in the morning until the radiators start to warm up, and another 4-8 minutes for the radiators themselves to warm up, let's look at how long it actually takes my apartment to warm up after the radiators heat up. I'll focus on the living room inside temperature for now.

Grab the Living Room Inside Temperature (index 9) from the tempF matrix.

livingroomtemp = tempF(:,9);

Let's look at one day's worth of living room inside temperatures to see how they compare the the living room radiator temperatures.

figure
[ax,p1,p2] = plotyy(ts,radiatorTemp,ts,livingroomtemp);
set(ax,'XLim',oneday)
set(ax(1),'YColor','k','YLim',[60 180],'YTick',60:20:180)
set(ax(2),'YColor','b','YLim',[60  72],'YTick',60:2:72,'XTick',[])
set(p1,'Color','k')
set(p2,'Color','b')
datetick(ax(1),'keeplimits')
xlabel(ax(1),'Time')
ylabel(ax(1),'Radiator Temperature (\circF)')
ylabel(ax(2),'Room Temperature (\circF)')
title('Living Room and Radiator Temperatures')
legend([p1,p2],{'Radiator Temperature','Room Temperature'},'Location','SouthEast')

Figure 13: Living room and radiator temperatures plotted overlapping but with different y-axes scaling.

As you would expect, once the radiator temperatures rise, so do the room temperatures, but I would like to find out how quickly the room temperatures rise. To do that, I will first break up the temperature data into segments delimited by the times the heat came on. I'll keep only the rising portion of each segment (up to the first occurrence of the maximum temperature in the segment, or until the heat is off). I'll also keep track of the minimum and maximum temperature in the segment, and at what time those temperatures occur. I am initializing a matrix to store the segments (segmentTemps) so that segments that are shorter than the maximum length are padded with NaN values instead of zeros.

numSegments = sum(heaton);
segmentSizes = heatoffind-heatonind+1;
segmentTemps = NaN(numSegments, max(segmentSizes));
for c = 1:numSegments
    segmentTemps(c,1:segmentSizes(c)) = livingroomtemp(heatonind(c):heatoffind(c));
    [segmentMaxTemp(c,1),segmentTimeToMax(c,1)] = max(segmentTemps(c,:));
    [segmentMinTemp(c,1),segmentTimeToMin(c,1)] = min(segmentTemps(c,1:segmentTimeToMax(c)));
    segmentTemps(c,segmentTimeToMax(c)+1:end) = NaN;
end

% Clip off columns that are all NaN at the end of the matrix.
maxSegmentSize = max(segmentTimeToMax);
segmentTemps = segmentTemps(:,1:maxSegmentSize);

% Shift from 1-based indexes to 0-based minutes.
segmentTimeToMin = (segmentTimeToMin - 1)*2;
segmentTimeToMax = (segmentTimeToMax - 1)*2;
segmentTimes = (0:(maxSegmentSize-1))*2; % Time since heat came on in minutes

Let's take a look at one of these segments:

figure
plot(segmentTimes, segmentTemps(1,:),'b.-')
hold on
plot(segmentTimeToMin(1),segmentMinTemp(1),'r.','MarkerSize',20)
legend({'Temperature','Minimum Temperature'},'Location','NorthWest')
xlabel('Minutes since radiator started to warm')
ylabel('Temperature (\circF)')
title('Room Temperature During One Heating Event')

Figure 14: Room temperature during one heating event, time zero is when the radiator started to warm up.

In the plot you can see that the temperature continues to decrease for a few minutes, then seems to rise almost linearly. Based on this information, let's shift the plot so that the base of the linear increase in temperature is at the origin.

figure
plot(segmentTimes-segmentTimeToMin(1), segmentTemps(1,:)-segmentMinTemp(1),'b.-')
legend({'Temperature Increase'},'Location','NorthWest')
xlabel('Minutes since minimum temperature')
ylabel('Temperature Increase (\circF)')
title('Shifted Room Temperature During One Heating Event')

Figure 15: Change in room temperature during one heating event, time zero is when the room temperature reached its minimum value.

That was just one segment, let's look at all of them, with the data shifted so that the lowest temperature in each segment occurs at the origin.

First I shift all the times based on the segmentTimeToMin

segmentTimes_mat = repmat(segmentTimes,length(segmentTimeToMin),1);
segmentTimeToMin_mat = repmat(segmentTimeToMin,1,length(segmentTimes));
segmentTimesShifted = segmentTimes_mat - segmentTimeToMin_mat;

Then shift all the temperatures based on the segmentMinTemp

segmentMinTemp_mat = repmat(segmentMinTemp,1,size(segmentTemps,2));
segmentTempsShifted = segmentTemps - segmentMinTemp_mat;

And find the total temperature change and time for each segment.

segmentRiseTime = segmentTimeToMax-segmentTimeToMin;
segmentTempRise = segmentMaxTemp-segmentMinTemp;

Now I can plot the shifted data

figure
h1 = plot(segmentTimesShifted',segmentTempsShifted','b-');
legend(h1(1),{'Temperature Increase'},'Location','NorthWest')
xlim([-10 max(segmentTimesShifted(:))])
xlabel('Minutes since minimum temperature')
ylabel('Temperature Increase (\circF)')
title('Shifted Room Temperature During All Heating Events')

Figure 16: Change in room temperature during all heating events, time zero is when the room temperature reached its minimum value.

Although it isn't perfect, it looks close to a linear relationship. Since I am interested in the time it takes to reach the desired temperature (what could be considered the "specific heat capacity" of the room), let me replot the data with time on the y-axis and temperature on the x-axis (swapping the axes from the previous figure). I'll also plot the data as individual points instead of lines, because that is how the data is going to be fed into polyfit later.

% Remove temperatures occuring before the minimum temperature.
segmentTempsShifted(segmentTimesShifted<0) = NaN;

figure
h1 = plot(segmentTempsShifted',segmentTimesShifted','k.');
xlabel('Temperature Increase (\circF)')
ylabel('Minutes since minimum temperature')
title('Time to Heat Living Room')
snapnow

Figure 17: The time it takes to heat the living room (axes flipped from Figure 16).

Now let me fit a line to the data so I can get an equation for the time it takes to heat the living room.

First I collect all the time and temperature data into a single column vector and remove NaN values.

allTimes = segmentTimesShifted(:);
allTemps = segmentTempsShifted(:);
allTimes(isnan(allTemps)) = [];
allTemps(isnan(allTemps)) = [];

Then I can fit a line to the data.

linfit = polyfit(allTemps,allTimes,1);

Let's see how well we fit the data.

hold on
h2 = plot(xlim,polyval(linfit,xlim),'r-');
linfitstr = sprintf('Linear Fit (y = %.1f*x + %.1f)',linfit(1),linfit(2));
legend([ h1(1), h2(1) ],{'Data',linfitstr},'Location','NorthWest')

Figure 18: The time it takes to heat the living room along with a linear fit to the data.

Not a bad fit. Looking closer at the coefficients from the linear fit, it looks like it takes about 3 minutes after the radiators start to heat up for the room to start to warm up. After that, it takes about 5 minutes for each degree of temperature increase.

What room takes the longest to warm up?

I can apply the techniques above to each room to find out how long each room takes to warm up. I took the code above and put it into a separate function called temperatureAnalysis, and applied that to each inside temperature sensor.

inside = [1 5 9 12];

figure
xl = [0 14];
for s = 1:size(inside,2)
    linfits(s,1:2) = temperatureAnalysis(tempF(:,inside(s)), heaton, heatoff);
    y = polyval(linfits(s,1:2),xl) + heatondelay;
    plot(xl, y, lineSpecs{inside(s),1}, 'Color',lineSpecs{inside(s),2},...
        'DisplayName',location{inside(s)})
    hold on
end
legend('Location','NorthWest')
xlabel('Desired temperature increase (\circF)')
ylabel('Estimated minutes to heat')
title('Estimated Time to Heat Each Room')

Figure 19: The estimated time it takes to heat each room in my apartment.

In part one of this series of blog posts, I mentioned that "I often find that even though the thermostat shows a nice and toasty 73, I'm sitting in my bedroom or the office freezing." Given that my thermostat is located in the living room, the data backs up my anecdotal observation: the office and bedroom take almost twice as long to warm up as the living room, while the kitchen warms up a little faster than the living room.

Conclusions

I started this series of blog posts with the intention of using a wireless network of temperature sensors (built using XBee® modules) to get a better understanding of how long my heating system takes to warm up my apartment, and to use that knowledge to help tune my heating system to provide more uniform heating.

The analysis above is just the start of the analysis that can be done on this data, but it is enough to help me to tune my heating system. I'll start by opening the valves on the bedroom and office radiators more, and maybe even close the valves on the living room and kitchen radiators a little. Then I'll collect data for another week and see whether that helps.

I now also know that my heating system takes around 15-20 minutes from when the heat comes on until the temperature starts to rise, and at best it takes about 3-5 minutes for each degree of temperature increase. My thermostat is set to allow the temperature to drop to 62 degrees at night, and warm up to 68 in the morning. Looking up 6 degrees on that last plot, I should adjust my thermostat to come on at least 30-40 minutes before I get up in the morning so it has enough time to warm up. No wonder my heating bill is so high.