Estimating air changes per hour with a blower door

When I was trying to figure out how big a gas boiler we needed for our house a few months ago, I tried to account for both the insulation in our walls as well as the air leaks that let warm air escape as cold air sneaks in. I had read that an old, drafty house turns over its volume in air once per hour. That seemed high to me, but I was looking for a conservative estimate, so that's what I used in my calculations. Since then, I've been hoping to find a way to make a better estimate.

Solution: Colin's blower door

A friend of mine from Stanford, Colin Fay, runs a company with his dad, David Fay, called Energy Metrics. Last weekend, Colin and his nearly homonymic associate Cullen were kind enough to bring Colin's blower door over to our house to run a test to see how drafty our house is.

The basic idea of a blower door is this: you fill your front door with a curtain and a massive fan that forces air out of the house. While it's doing that, a small sensor measures the air pressure difference between the inside and outside of the house. There are a few different tests you can run, but the standard test that the fan controller runs is to automatically adjust the fan speed until the pressure inside is 50 Pa lower than outside. For comparison: 50 Pa is roughly equivalent to the pressure from a windspeed of 20 mph, but blowing at your house uniformly from all directions. Atmospheric pressure is around 100,000 Pa.

When the fan reaches a steady state, air is whistling in through all the gaps around your windows, doors and foundation, and you can tell where the problems are. For us, the largest draft was coming under the basement door. The next worst were the gaps between the sashes in our larger, older double-hung windows. In real life, I suspect that the gap under the basement door is not so bad-- the thermal gradient keeps the colder, denser air sunk down in the basement. I didn't realize it at the time, but most of the draft was probably coming down through our vestigial chimney.

Colin's blower door, a Retrotec 2000 with a DM-2 Mark II controller, pulled air through our house at 3900-4000 ft³ per minute to generate a pressure difference of 50 Pa. Our house has a volume of around 18000 ft³, so with the fan blowing, we were replacing all the air in our house every 4.5 minutes, or 13.3 times per hour.

Assembling the blower door frame

The blower door from the inside

Once you know how drafty your house is with a fan pulling the heavens through it, you need to scale that to match the typical conditions for your house. As a rough rule of thumb: just divide by 20. With the fan, we had 13.3 air changes per hour, so that's about 0.7 air changes per hour without the fan.

But if you want to ascend to the peak of Mount Energygeek, and you're willing to do it unashamedly, you can use the empirical corrections of Max Sherman of the Energy Performance of Buildings Group at Lawrence Berkeley National Lab, who completed his thesis on modeling building air infiltration in 1980, when oil rolled down like waters and righteousness like acid rain. You look up correction factors for climate (~18 for Boston), building height (0.7 for our house), wind shielding (1) and leakiness (1), multiply them together, and you've got a better correction factor than the rough guess of 20. For our house, we end up with 13.3/(18 x 0.7 x 1 x 1) = 1.06 air changes per hour.

With that knowledge, you can calculate the power required to offset the drafts cooling or heating your house. Our house, nominally a 1900 ft² Victorian, has an internal volume of 18000 ft³, or 510 m³, so when it's 0 C outside, we're heating about 1.06 x 510 m³ of air per hour by around 20 C. The heat capacity of air is around 1200 J/m³C. That means we need to pour in 1200 J/m³C x 540 m³ x 20 C every 3600 seconds. By my calculation, that's about 3.6 kW.

The conductive heat loss model I developed for our house a few months ago when we were installing the boiler predicts that the conductive heat loss at the same temperature will be around 18 kW, so we lose about 1/6th of our heat from air infiltration.

Colin suggested we could reduce our draftiness by around 2x before we'd have to worry about the effects of too little fresh air (farts, basically). He suggested picking up a tube of transparent silicone caulk in the fall to fill the gap between the window sashes, as that's where our worst leaks are. In the spring, when it's time to open the windows again, the silicone peels off.

After seeing fellow energy geek Holly's sexy basement windows last weekend, I think I might look into replacing those as well.

Sizing a gas boiler

When you own a house, you are forced into bizarre contortions to reduce the cost of maintenance. For example, our roof needs to be replaced, but part of that is repairing the chimney. We don't really even need the chimney, because we don't have any fireplaces left-- they were all removed decades ago, well before we owned the house. But, until late this fall, when we installed a new boiler and water tank, we still needed a way to vent our ancient and inefficient gas boiler and water heater, or the basement would fill up with carbon dioxide and soot. We realized that if we replace those first, we can buy efficient systems that need only a plastic pipe through a basement wall to vent, rather than a chimney. Our heating cost would drop, and we could just toss the chimney in a dumpster, which is cheaper than repairing it, and need only be done once. Also, when we do replace the roof, we don't have to worry about finding someone with the skills to make a reliable flashing around the chimney.

As a result, we decided to replace our gas boiler and water heater. (The water heater actually had some life left in it, but we just wanted to do this once and be done with it.) An added advantage is that we could tear out the matrix of steel pipes that undergirds our first floor and replace them with PEX tubing, so I no longer have to wear a helmet when I walk around the basement.

Knowing that we're going to replace the boiler and water heater, someone had to figure out how big the replacements should be. The water heater is fairly easy-- we had a 40 gallon tank, which was plenty for the two of us, given that our shower uses around 6 L (1.5 gallons) per minute. Given that we were replacing the boiler and the water heater at the same time, it made sense to install an indirect storage tank, rather than a separate water heater. With an indirect setup, the super-efficient boiler heats water to warm the radiators in the house, but it also has a second loop that runs through a heat exchanger in a water storage tank. The heat exchanger is necessary so you don't end up drinking water filled with radiator rust. The boiler manufacturer we considered most seriously, Viessmann, doesn't make a water storage tank smaller than 42 gallons, so that's what we chose. We also considered getting a boiler with an integral on-demand system (the WB2 6-24C), but the output looked a little shy of what we wanted. (She-who-takes-hot-showers overruled he-who-takes-long-showers on this one.)

Sizing the gas boiler is more difficult. Our old boiler was rated at 48 kW (164,000 BTU/h or 164 MBH). Our old boiler didn't modulate-- it either burned at full power or turned off. Using our gas meter and a stopwatch, I measured the full-power consumption at 3.0 ft³ per minute, which is 53.5 kW for natural gas, slightly above the nominal rating of the boiler. On a cold evening recently, when the temperature was holding constant for several hours around 33° F outside and 66° F inside, and there was no sun and little wind, I measured the duty cycle of our old boiler. Over a 133-minute period, it was burning 26 minutes, which is about 20%. This means that average power going into the boiler with a temperature difference of 33° F between inside and outside was 11 kW. Maybe 8 or 9 kW of that actually makes it into the house, while the rest goes up the chimney.

The choices

I was trying to choose between the Viessmann WB2 6-24 and the larger 8-32. (The parts of the names of the models after the hyphens are close to, but not equal to, their peak output in kW-- perhaps they were the design targets.) Both models are extremely efficient condensing boilers; depending on how they're loaded, they operate in the 90-98% efficient range. As they're loaded closer to peak power, their efficiency drops. It also drops with the returning water temperature-- if the water comes back from the radiators cold, it's draws in heat from the burner more efficiently than if it comes back warm. The radiators in our house were installed before insulation and storm windows were standard, so they're quite effective at distributing heat.

This suggested to me that our new boiler would probably be operating near the high end of its efficiency range; the real danger was that in the fall and spring, when just a little heat was needed, it would inefficiently cycle on and off, because its lowest level of modulation would be too high. The larger of the two models I was choosing between cut off at 11 kW, which, you'll recall, was what we were using at 33° F with a much less efficient boiler, so that made the smaller one seem like a better fit. (For those of you unfamiliar with the weather in Boston, the average daily high temperature in the coldest month is 36° F in January.) My estimation was that on a really cold night, we might have a temperatures outside in the single digits, and we'd keep the house in the upper 50's or low 60's, so call that a delta of 50 degrees, about 1.5 times the 8-9 kW we needed at 33 F-- call it 14 kW to be safe.

The last factor I thought about was wind. Our house is not particularly exposed, but I suspect a cold wind could double the heat loss, which would definitely max out the smaller boiler, which peaked at about 25 kW, allowing for boiler inefficiency. If that were the end of the story, I would have chosen the larger boiler, but there were three other factors. The first was that the new boiler would zone the radiators, so I could turn off less critical rooms if we were short on heat. The second was that I'd rather invest money in insulation than a bigger furnace. Finally, we were thinking about installing a woodstove (for strictly recreational purposes-- Sundays, crossword puzzles, and so forth).

The results, in which we learn whether Brandon is an idiot

Against the advice of the heating contractor, we got the smaller boiler. They did a beautiful job installing it in November. Last week, we had a very cold (sub-10° F), very windy night, and the system kept up with demand. We just got our first gas bill covering a month with the new boiler, and our gas usage dropped to 150 therms from around 250 therms for the same time last year, and last year was slightly warmer, so I'm pretty happy. The payback time on the added expense above the cost of an inefficient system will be less than 10 years, which is less than the warranty on the system (except the control electronics). The water storage tank, which is stainless steel, is warranted for the lifetime of the original owner in the original house.

I'd estimate that about half of the savings are from reductions in heated area (spare bedroom, attic) and half are from higher efficiency.

Shortcomings

But nothing is perfect. Because we replaced all the large, cast iron pipes in the basement with PEX tubing, the basement hovers in the low 50s in the winter, rather than in the low 60s. This was expected, but I decided I'd rather wear a winter hat than a helmet. Also, because we're not losing heat throughout the system, the boiler runs a little cooler than the old beast, which means that the radiant floor heat in the kitchen has trouble keeping up on the coldest days. I think I can probably fix that by tweaking the boiler target temperature and the circulation pump flow rates.

One other wrinkle is the control scheme used by the new system. Viessmann boilers come with an external temperature sensor, which is mounted outside on the north side of the house. Viessmann's intent (I think) is that the installer calibrates the system to drive the temperature of the hot water supplied to the radiators as a linear function of the outside temperature and the return temperature of the water. The idea is that when the temperature outside changes drastically, the boiler can react earlier than if it had to wait for the effect to be detectable in the radiator water or thermostats. Strangely, our boiler is set up to use the two programmable thermostats we used with the old boiler. The brain of the boiler targets 75° F, but the circulation pumps for the radiators and the floor heat only turn on when the thermostats call for heat. I'm not sure that this is actually worse, but it sure is strange to install an external sensor if you're not planning on using the signal.

Pictures

Here are the old system and the new system. If you click on the pictures of the new system, you'll get to Flickr, where I've added a few notes about what the different parts in the picture do.

Old system with headbanging pipes

New system

Supporting plumbing for new system

Least squares fit of a surface to a 3D cloud of points in Python (with ridiculous application)

The floor in the room above the kitchen in our house has a floor that slopes almost 1 inch per foot for half of the room. Experimentally, we have found that this is steep enough to make a desk chair roll– kind of irritating, particularly to my special ladyfriend who happens to occupy such a chair in that zone.

To compensate for the slope, I decided to fit a stack of masonite sheets to the curve of the floor. Unfortunately, the floor slopes nonlinearly in two directions, like the rounded corner of a swimming pool. After making a series of measurements of the floor, I decided to fit a polynomial in two variables to the cloud of points using a least squares estimate.

(In retrospect, the floor was close enough to singly-curved that I could have gotten away with a linear fit.) The blue thing in the picture is a level. I shimmed the level until it was worthy of its name, and then measured the distance to the floor with calipers.

Estimating the flatness of the floor

I’ll recount the basics from my earlier post about least squares fitting in Python. Skip ahead to the next section if you read that already.

As before, the first step is to arrange the equations in canonical form: Ax=y where:

* A is an M x N full rank matrix with M > N (skinny)
* x is a vector of length N
* y is a vector of length M

As matrices, this looks like

In polynomial fitting, A is called the Vandermonde matrix and takes the form:

The 3D case

In the 2D case, we’re trying to find polynomial in x such that f(x) approximates y. In the 3D case at hand, we have two independent variables, so we’re looking for a polynomial in x and y such that f(x, y) approximates z. Rather than the 2D case:

we want the final output to look like this:

(Spike Curtis astutely notes in the comments that I am omitting the cross terms, such as xy, and refers us to a Matlab script. Spike is right, but the floor’s already fixed.)

We can use two Vandermonde matrices next to each other.

Here’s the Python code that creates the two Vandermonde matrices and joins them into one matrix. x, y, and z are lists of corresponding coordinates, so, for example, x[5], y[5] and z[5] are the coordinates of one point that the surface should approximate. The order of the points is not important.

    import numpy as np
     
    z = [0.0, 0.695, 1.345, 1.865, 2.225, 2.590, 0.0, 0.719, 1.405, 1.978, 2.398, 2.730, 0.0, 0.789, 1.474, 2.064, 2.472, 2.775, 0.0, 0.763, 1.453, 1.968, 2.356, 2.649]
     
    x = [0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 0.0, 12.0, 24.0, 36.0, 48.0, 60.0, 0.0, 12.0, 24.0, 36.0, 48.0, 60.0]
     
    y = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 48.0, 48.0, 48.0, 48.0, 48.0, 48.0, 72.0, 72.0, 72.0, 72.0, 72.0, 72.0]
     
    degree = 3
    thickness = 0.167
     
    # Set up the canonical least squares form
    Ax = np.vander(x, degree)
    Ay = np.vander(y, degree)
    A = np.hstack((Ax, Ay))
     
    # Solve for a least squares estimate
    (coeffs, residuals, rank, sing_vals) = np.linalg.lstsq(A, z)
     
    # Extract coefficients and create polynomials in x and y
    xcoeffs = coeffs[0:degree]
    ycoeffs = coeffs[degree:2 * degree]
     
    fx = np.poly1d(xcoeffs)
    fy = np.poly1d(ycoeffs)

Once I knew the coefficients of the polynomial approximation, I could calculate the contours of the masonite. From a measurement of the stack of masonite, I knew the average thickness is 0.167 inches. (Thanks to Dr. Alex T. Tung, Ph.D., for helping me get the masonite home from Home Depot.) To find out where the first layer should end, I picked a series of stations spaced every 12 inches along the length of the floor in the y-direction. Along those stations, I solved f(x, y) = 0.167, f(x, y) = 2 * 0.167, f(x, y) = 3 * 0.167 and so forth.

In practice, solving f(x, y) = c, where c is a constant, means finding the roots of the equation f(x, y) - c = 0. (The mathematicians call this solving the homogeneous equation.) In Python, the numpy.roots method solves the homogeneous case. For each contour/section crossing, I generated a polynomial of the form f(x, y) - c and solved it with numpy.roots.

    ystations = range(0, 84, 12)
    sections = [[np.poly1d(xcoeffs - [0,0,zoffset - fy(ypos)]) for zoffset in np.arange(thickness, max(z), thickness).tolist()] for ypos in ystations]
    pts = [[min(func.roots) for func in list_of_fs] for list_of_fs in sections]

For fabrication, I printed out a list of the locations where the masonite contours crossed the stations.

    for (pt_list, ystation) in zip(pts, ystations):
        print('\nBoundaries at station y = {0} inches:'.format(ystation))
        print('\t'.join(['{0:.3}'.format(pt) for pt in pt_list]))

Armed with my list of measurements, I headed to the garage and set up some sawhorses with a sheet of plywood to keep the masonite from bowing and flopping around. It took a few hours of marking points and cutting gentle curves with a jigsaw, but the results were delightful.

The masonite platform (non-impressive view)

Level and gleeful

Side view of the masonite platform

Another side view of the masonite platform

As you can see above, I haven’t fastened the layers together yet. They seem to be sticking together reasonably well.

In the end, I think the slope went from about a 2.75″ drop across the 48″ span, to within 1/8″ of flat across the same distance. Sharon was delighted, and I found the whole project both more time-consuming and more satisfying than I expected.

Least squares polynomial fitting in Python

A few weeks ago at work, in the course of reverse engineering a dissolved oxygen sensor calibration routine, Jon needed to fit a curve to measured data so he could calculate calibration values from sensor readings, or something like that. I don’t know how he did the fitting, but it got me thinking about least squares and mathematical modeling.

The tools I’ve used for modeling (Matlab, Excel and the like) and the programming languages I’ve used for embedded systems (C and Python, mostly) are disjoint sets. Who wants to use Matlab for network programming, when you could use Python instead? Who wants to use C for calculating eigenvalues, when you could use Matlab instead? Generally, nobody.

Unfortunately, sometimes you need to do both in the same system. I’m not sure that I’m making a wise investment, but I think that Python’s numerical libraries have matured enough that Python can straddle the divide.

So, here’s what I’ve figured out about least squares in Python. My apologies to most of you regular readers, as you will find this boring as hell. This one’s for me and the internet.

Least squares estimation

Generally, least squares estimation is useful when you have a set of unknown parameters in a mathematical model. You want to estimate parameters that minimize the errors between the model and your data. The canonical form is:

Ax=y where:

* A is an M \times N full rank matrix with M > N (skinny)
* x is a vector of length N
* y is a vector of length M

As matrices, this looks like

You want to minimize the L2-norm, ||Ax-y||. To be explicit, this means that if you define the residuals r = Ax - y, you want to choose the that minimizes sqrt{r_1^2 + ... r_M^2}. This is where the “least squares” name comes from– you’re going to choose the that results in the squares of r being least.

Once you know A and y, actually finding the best x (assuming you have a computer at your disposal) is quick. The derivation sets the gradient of ||r||^2 with respect to x to zero and then solves for x. Mathematicians call the solution the Moore-Penrose pseudo-inverse. I call it “what the least squares function returned.” In any case, here it is:

Fortunately, one need not know the correct name and derivation of a thing for it to be useful.

So now, how do you go about using this technique to fit a polynomial to some data?

Polynomial fitting

In the case of polynomial fitting, let’s say you have a series of data points, (t_i, y_i) that you got from an experiment. The t_i are the times at which you made your observations. Each y_i is the reading from your sensor that you noted at the time with the same i.

Now you want to find a polynomial y = f(t) that approximates your data as closely as possible. The final output will look something like

where you have cleverly chosen the x_i for the best fit to your data. You’ll be able to plug in whatever value of t you want and get an estimate of y in return.

In the conventional form, we are trying to find an approximate solution to Ax=y

In this case, A is called the Vandermonde matrix; it takes the form:

(This presumes that you want to use powers of t as your basis functions– not a good idea for all modeling tasks, but this blog post is about polynomial fitting. Another popular choice is sinusoidal functions; perhaps Google can tell you more on that topic.)

(Also, the Vandermonde matrix is sometimes defined with the powers ascending from left to right, rather than descending. Descending works better with Python because of the coefficient order that numpy.poly1d() expects, but mathematically either way works.)

(Back to polynomial fitting.)

Each observation in your experiment corresponds to a row of A; A_i * x = y_i. Presumably, since every experiment has some noise in it, there is no x that fits every observation exactly. You know all the t_i, so you know A.

As a reminder, our model looks like this:

You can check that A makes sense by imagining the matrix multiplication of each row of A with x to equal the corresponding entry in y. For the second row of A, we have y_2 = x_Nt_2 + . . . + x_2t_2^2 + x_1t_2^1 + x_0

which is just the polynomial we’re looking to fit to the data.

Now, how can we actually do this in Python?

Python code

Let’s suppose that we have some data that looks like a noisy parabola, and we want to fit a polynomial of degree 5 to it. (I chose 5 randomly; it’s a stupid choice. More on selection of polynomial degree in another post.)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

degree = 5

# generate a noisy parabola
   t = np.linspace(0,100,200)
   parabola = t2
   noise = np.random.normal(0,300,200)
   y = parabola + noise
   

So now we have some fake data– an array of times and an array of noisy sensor readings. Next, we form the Vandermonde matrix and find an approximate solution for x. Note that numpy.linalg.lstsq() returns a tuple; we’re really only interested in the first element, which is the array of coefficients.

# form the Vandermonde matrix
A = np.vander(t, degree)
 
# find the x that minimizes the norm of Ax-y
(coeffs, residuals, rank, sing_vals) = np.linalg.lstsq(A, y)
 
# create a polynomial using coefficients
f = np.poly1d(coeffs)

Three lines of code later, we have a solution!

I could have also used numpy.polyfit() in place of numpy.linalg.lstsq(), but then I wouldn’t be able to write a blog post about regularized least squares later.

Now we’ll plot the data and the fitted curve to see if the function actually works!

# for plot, estimate y for each observation time
y_est = f(t)
 
# create plot
plt.plot(t, y, '.', label = 'original data', markersize=5)
plt.plot(t, y_est, 'o-', label = 'estimate', markersize=1)
plt.xlabel('time')
plt.ylabel('sensor readings')
plt.title('least squares fit of degree 5')
plt.savefig('sample.png')

This post took two-thirds of eternity to write. Feel free to ask questions or point out how confused I am below.

Next: regularized least squares! Or maybe selection of polynomial degree! Or maybe forgetting the blog and melting brain with TV!