Anyone around Boston have some free/broken hedge trimmers?

June 28th, 2009

Greetings, fellow internet denizen! I am concocting a device that I call The Electric Lawn Mower That Doesn’t Suck.

I have an electric lawn mower that is large, heavy, almost as loud as a gas mower, and only works with a cord attached. I also have a reel mower that is pretty sweet, but it can’t cut grass taller than approximately the radius of the blade cylinder, so it leaves stragglers.

I think if I could mount a hedge trimmer on something like a reel mower and power it off a cordless drill battery, I’d be golden. Hedge trimmers are pretty loud, but I think I could rework the drive mechanism so that it was quiet (though it would cost a bit more).

To this end, if you have a crappy hedge trimmer that you would be willing to surrender to my prototyping efforts, I would be happy to take it off your hands.

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

June 15th, 2009

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 \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

 \left[<br />
  \begin{array}{c c c}<br />
  \\<br />
  \\<br />
  \quad & A & \quad \\<br />
  \\<br />
  \\<br />
  \end{array}<br />
\right]<br />
\left[<br />
  \begin{array}{c}<br />
  \\<br />
  x \\<br />
  \\<br />
  \end{array}<br />
\right] = \left[<br />
  \begin{array}{c}<br />
  \\<br />
  \\<br />
  y \\<br />
  \\<br />
  \\<br />
  \end{array}<br />
\right]

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

 \left[<br />
  \begin{array}{ c c c c c}<br />
     t_1^N & \cdots & t_1^2 & t_1 & 1 \\<br />
     t_2^N & \cdots & t_2^2 & t_2 & 1 \\<br />
     \vdots \\<br />
     t_M^N & \cdots & t_M^2 & t_M & 1 \\<br />
  \end{array}<br />
\right]

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:

y = x_Nt^N + . . . + x_2 t^2 + x_1t + x_0

we want the final output to look like this:

z = a_Nx^N + . . . + a_2 x^2 + a_1x + a_0 + b_Ny^N + . . . + b_2 y^2 + b_1y + b_0

(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.

 \left[<br />
  \begin{array}{ c c c c c c c c c c}<br />
     x_1^N & \cdots & x_1^2 & x_1 & 1 & y_1^N & \cdots & y_1^2 & y_1 & 1 \\<br />
     x_2^N & \cdots & x_2^2 & x_2 & 1 & y_2^N & \cdots & y_2^2 & y_2 & 1 \\<br />
     \vdots \\<br />
     x_M^N & \cdots & x_M^2 & x_M & 1 & y_M^N & \cdots & y_M^2 & y_M & 1 \\<br />
  \end{array}<br />
\right]<br />
\left[<br />
  \begin{array}{c}<br />
  a_N \\<br />
  \vdots \\<br />
  a_2 \\<br />
  a_1 \\<br />
  a_0 \\<br />
  b_N \\<br />
  \vdots \\<br />
  b_2 \\<br />
  b_1 \\<br />
  b_0 \\<br />
  \end{array}<br />
\right] = \left[<br />
  \begin{array}{c}<br />
  z_1 \\<br />
  z_2 \\<br />
  z_3 \\<br />
  \vdots \\<br />
  z_M \\<br />
  \end{array}<br />
\right]

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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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.

26
27
28
29
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.

30
31
32
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.

Opening salvo in the Great Skunk War of 2009

June 1st, 2009

Our house in Somerville has two porches that serve as barracks for at least two skunks.

Today, Sharon and I fired the first salvo, using a piece of chicken wire (AKA “hardware cloth”) to prevent Merle (the commander-in-chief of the skunk army) from tunneling under the lattice work under the front porch.

Trench warfare

Trench warfare


Historic photos from the trench warfare occurring at the northern front have reached Flickr.

Update: As of the next morning, the fortifications had not been breached! However, there was a whiff of skunk in the front hall for a few minutes. It may be that we have misjudged which side of enemy lines we were on, or more specifically, which side the enemy was on.

Homemade edging tool

May 3rd, 2009

Yesterday, I welded up an edging tool for making the grass along our sidewalk conform to my strict straightness requirements. I had been using it for 15 or 20 minutes before I realized that using an edging tool on your grass borders on the neurotic, and making your own edging tool is even weirder. But, I had already made it, and it came out pretty well, so here are some pictures.

The handle and the blue steel yoke come from the failed edger that I had, which clogged with grass and dirt almost instantaneously.

Grim power embodied in the form of an edger

Pictures of the construction process are on Flickr.

The garden, after a savage edging:
The garden, after a savage edging

Evaluating renewable energy companies

March 24th, 2009

As part of my job at GreenMountain, I work with a lot of renewable energy startups. As a result, my friends are constantly telling me about new startups that they’re excited about. This often leads to a discussion of how I evaluate renewable energy companies. The companies that my friends think are on the way to the top are rarely the ones that I think will win, so I think I should try to explain myself.

First off, my goal is to predict whether the company at hand will be able to produce huge amounts of energy from renewable sources at prices that rival fossil fuels. In most of the world, this is not currently possible– fossil fuels are still cheap. There are regional exceptions, like residential solar power in ridiculously sunny places like the American southwest, but generally, no renewable energy company has yet succeeded.

Thus, the interesting question is: what do you need to do to win?

You must scale to huge.

A few years ago, I worked at an environmental foundation in Maine. For the last 8 years, the foundation has produced their own biodiesel from excess frialator oil from two nearby seafood restaurants that boomed in the summer. While I worked there, I was able to buy biodiesel for my Passat at whatever price the local gas station was charging for regular diesel.

It was brilliant, but the only reason it worked was that we got the frialator oil for free. I personally burned around 10% (50 gallons) of the annual biodiesel yield from the two seafood restaurants. I was driving fewer miles than the national average, and in a 42 mpg car. It’s still a great idea, but biodiesel from frialator oil will never scale to huge.

Breakthroughs in the lab are not enough.

Every few weeks, I hear news of a breakthrough in solar cell technology in a laboratory setting. Usually, it’s the Grätzel cell again. If not, it’s often a new manufacturing technique that is purported to be cheap in high volume, but which is currently being produced in very low volume. Even the largest cell producers don’t have a good idea of how cheap their cells will be when they scale beyond GW levels. (Maybe they’re pretty good at 10 GW, but that’s pushing it.) Claiming to be able to make those predictions accurately when you’re below 1 kW is not credible.

You will not challenge Vestas soon.

In the more established renewable industries, like big wind, there are just a few players. Vestas, for example, sells about 25% of the large wind turbines in the world. They have around 20,000 employees and have been producing turbines for 30 years. They do over €1 BN in sales per quarter. Designing a MW-scale turbine takes years, never mind testing.

Most startups could not fit one large turbine blade in their offices. If you’re a big wind startup with a plan to compete on that scale in the next few years, you are not credible. You might be able to do it, but you need a ten year plan.

Higher efficiency at higher cost is not an automatic win.

In small wind, a common approach is the shrouded wind turbine. The general idea is to add a shroud around the wind turbine that channels wind into the blades. So far, nobody has been successful with this approach– the economics dictate that the marginal cost of blade extension beats that of shrouding. This is not to say that shrouded turbines will always lose, just that so far, nobody has been able to get an increase in efficiency that outweighs the increase in cost.

Concentrating photovoltaic systems face the same problem– can you add a concentrating mirror cheaply enough that it’s worth the cost? The jury is still out on that one.

Lower cost with lower efficiency is not an automatic win.

For years, solar cell companies have been trying to make “thin film” solar cells, i.e. cells that are manufactured by depositing a coating a few microns thick on a substrate, usually metal. Because thin film cells are so thin, they use very little material, so they’re much cheaper per unit of area than conventional cells. Unfortunately, they’re also less efficient. In the long run, thin film cells may still win, but so far it seems that, despite massive investments in R&D, the efficiency penalty outweighs the cost savings.

Every engineering challenge comes with a counterpart in marketing.

To win in renewable energy, you can’t be just a brilliant engineer. You also need a team that can sell your products. It’s like preventing car theft– it doesn’t do any good to lock one door of your car, even if you lock it really carefully, and the lock is really strong. You have to lock all the doors, every time you leave the car. In the same vein, your heat recovery system may be brilliant, but if you marketing efforts look amateurish, you’re sunk.

The ocean, sun, wind, and rain are harsh.

Most energy systems have to survive outside for their ~20 year lifetime. If they can’t, you have lost. Ocean power systems, for example, have to survive scouring from ice and abuse from drifting trees. Wind turbines have to survive hurricanes; solar panels must endure hail. If you’re a slippery businessman, you might be able to sell a product good enough to give you 5 years to change your identity and escape to a country missing the requisite extradition treaty, but that won’t solve our energy problem.

Please note in the comments respectable axes of evaluation that I have omitted.

Calculating solar panel shading in Python

March 16th, 2009

Eventually, I’d like to install solar panels on our house, but I want to know that it will be worth it before I commit the money. Back in 2007, I wrote a library for calculating the path of the sun given the time and your location on the earth. Since then, I’ve been thinking about the next step, which is calculating how much neighboring houses and trees would obstruct the sun at different times of the day. In a nutshell, I wanted to replicate the calculations performed by devices like the Solmetric Suneye without paying $1500 for the device. To do the job right, I still need a quality circular fisheye lens ($680 new), but I’ve got a first approximation working with a $13 door peephole from McMaster.

(I do realize that I could just get a solar contractor to do all the site evaluation for me for free. That misses the point. If all I wanted was some sort of fiscal efficiency, I’d get a job as a financial executive or maybe something different from that, like a bank robber.)

Back to the task at hand.

Below is the image I started with, which I took back in the fall, before the leaves fell off the trees. This was using a Canon SD450 camera and the aforementioned peephole. I didn’t take particular care to level the camera or center the peephole.

Original peephole image

Using the Python Imaging Library, I cropped the photo to be square and switched to black and white mode. Then I mapped concentric rings in the image below to the subsequent straightened image.

Peephole image after mapping to rectilinear coordinate system

Here’s the straightening code. The guts are in the nested for loops near the end. (If there are any Pythonistas out there who know how to iterate over concentric rings using a list comprehension or map(), please let me know. The code below is functional and reasonable clear, but a little slow.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# Code under GPLv3; see pysolar.org for complete version.
def despherifyImage(im):
    (width, height) = im.size
    half_width = im.size[0]/2
    half_height = im.size[1]/2
 
    inpix = im.load()
    out = Image.new("L", (width, half_height))
    outpix = out.load()
 
    full_circle = 1000.0 * 2 * pi
 
    for r in range(half_width):
        for theta in range(int(full_circle)):
            (inx, iny) = (round(r * cos(theta/1000.0)) + half_width, round(r * sin(theta/1000.0)) + half_width)
            (outx, outy) = (width - width * (theta/full_circle) - 1, r)
            outpix[outx, outy] = inpix[inx, iny]
    return out

The straightening works pretty well, but there is a little distortion that I think is caused by the peephole not being a perfectly spherical lens. Also, because the peephole is not centered on the camera lens, my concentric ring transformation isn’t centered either.

From here, I needed to figure out where the sky ends and the buildings or trees start. Unfortunately, the buildings and trees are both lighter *and* darker than the sky in different places, so I can’t just look for one type of transition. To detect the edge, I scan down each column of pixels and calculate the difference in darkness between consecutive pixels, ignoring whether the change was from light to dark or the reverse. The image below shows those differences, amplified by 10x to increase the contrast. The mathematicians call this calculation a finite forward difference.

The derivative of each column of pixels

Here’s the finite difference code.

1
2
3
4
5
6
7
8
9
10
# Code under GPLv3; see pysolar.org for complete version.
def differentiateImageColumns(im):
    (width, height) = im.size
    pix = im.load()
 
    for x in range(width):
        for y in range(height - 1):
            pix[x, y] = min(10 * abs(pix[x, y] - pix[x, y + 1]), 255)
 
    return im

The last step is to scan down each column looking for the first large value. The first change that crosses a threshold is recorded. The final output is an array of values that measure the angle of the highest obstruction as a function of direction. As a sanity check, I drop a red dot on each value. It’s hard to make out in the thumbnail below, but if you click on the image below, you’ll get a larger version where you can see the red dots work pretty well.

Red dots highlighting the first large change in each column

I think the next step will be to calculate the total energy delivered per year using Pysolar.

Here’s the full code.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
# Code under GPLv3; see pysolar.org for complete version.
from PIL import Image
from math import *
import numpy as np
 
def squareImage(im):
    (width, height) = im.size
    box = ((width - height)/2, 0, (width + height)/2, height)
    return im.crop(box)
 
def despherifyImage(im):
    (width, height) = im.size
    half_width = im.size[0]/2
    half_height = im.size[1]/2
 
    inpix = im.load()
    out = Image.new("L", (width, half_height))
    outpix = out.load()
 
    full_circle = 1000.0 * 2 * pi
 
    for r in range(half_width):
        for theta in range(int(full_circle)):
            (inx, iny) = (round(r * cos(theta/1000.0)) + half_width, round(r * sin(theta/1000.0)) + half_width)
            (outx, outy) = (width - width * (theta/full_circle) - 1, r)
            outpix[outx, outy] = inpix[inx, iny]
    return out
 
def differentiateImageColumns(im):
    (width, height) = im.size
    pix = im.load()
 
    for x in range(width):
        for y in range(height - 1):
            pix[x, y] = min(10 * abs(pix[x, y] - pix[x, y + 1]), 255)
 
    return im
 
def redlineImage(im):
    (width, height) = im.size
    pix = im.load()
 
    threshold = 300
 
    for x in range(width):
        for y in range(height - 1):
            (R, G, B) = pix[x, y]
            if R + G + B > threshold:
                pix[x, y] = (255, 0, 0)
                break
    return im
 
im = Image.open('spherical.jpg').convert("L")
im = squareImage(im)
 
lin = despherifyImage(im)
d = differentiateImageColumns(lin).convert("RGB")
r = redlineImage(d)
 
r.show()

Least squares polynomial fitting in Python

January 24th, 2009

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

 \left[<br />
  \begin{array}{c c c}<br />
  \\<br />
  \\<br />
  \quad & A & \quad \\<br />
  \\<br />
  \\<br />
  \end{array}<br />
\right]<br />
\left[<br />
  \begin{array}{c}<br />
  \\<br />
  x \\<br />
  \\<br />
  \end{array}<br />
\right] = \left[<br />
  \begin{array}{c}<br />
  \\<br />
  \\<br />
  y \\<br />
  \\<br />
  \\<br />
  \end{array}<br />
\right]

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 x that minimizes \sqrt{r_1^2 + ... r_M^2}. This is where the “least squares” name comes from– you’re going to choose the x 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:

x_{ls} = (A^TA)^{-1}A^Ty

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

y = x_Nt^N + . . . + x_2 t^2 + x_1t + x_0

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:

 \left[<br />
  \begin{array}{ c c c c c}<br />
     t_1^N & \cdots & t_1^2 & t_1 & 1 \\<br />
     t_2^N & \cdots & t_2^2 & t_2 & 1 \\<br />
     \vdots \\<br />
     t_M^N & \cdots & t_M^2 & t_M & 1 \\<br />
  \end{array}<br />
\right]

(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:

 \left[<br />
  \begin{array}{c c c}<br />
  \\<br />
  \\<br />
  \quad & A & \quad \\<br />
  \\<br />
  \\<br />
  \end{array}<br />
\right]<br />
\left[<br />
  \begin{array}{c}<br />
  \\<br />
  x \\<br />
  \\<br />
  \end{array}<br />
\right] = \left[<br />
  \begin{array}{c}<br />
  \\<br />
  \\<br />
  y \\<br />
  \\<br />
  \\<br />
  \end{array}<br />
\right]

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.)

1
2
3
4
5
6
7
8
9
10
11
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 = t**2
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.

12
13
14
15
16
17
18
19
# 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!

20
21
22
23
24
25
26
27
28
29
# 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')
Least squares fit of a noisy parabola

Click the graph for a larger version.

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!

Slides from Business and Society Conference

January 17th, 2009

My colleague David Hague and I spoke at the Business and Society Conference at Dartmouth this week. In talking with folks after the panel ended, there was some interest in the slides I presented: 800 kB .ppt file.

Cable pulling lubricant exists

December 30th, 2008

A few weeks ago, I was attempting to run Cat 6 UTP cable for gigabit Ethernet from our house to the garage in a pre-existing conduit. After a very frustrating attempt involving a fishtape, various strings, a flashlight, and lots of scrambling around in the crawlspace under the kitchen, I was convinced that it was not possible to pull the cable through all the corners in the crowded conduit without damaging the cable. My next thought was, “I need . . . cable . . . lubricant. I bet that exists.” A quick trip to McMaster yielded 7431K41, a quart bottle of Ideal industries’ Yellow 77 Plus cable pulling lubricant.

Yellow 77 Plus cable pulling lubricant

The next weekend, I re-pulled a new cable slathered in Yellow 77 Plus. The force required was reduced by an order of magnitude. I was gleeful.

Anyway, I just wanted to help distribute global knowledge a little more evenly: cable pulling lubricant exists! Also: gigabit Ethernet to the garage!

Becoming a renewable energy engineer

December 29th, 2008

Due to an unlikely confluence of world energy trends, my interest in engineering, and the internet, lots of people ask me how they can get my job as a renewable energy engineer. Most of them don’t want my job literally, but they do want to work on the engineering side of renewable energy. I generally don’t mind answering that kind of question, but I’m busy enough that I thought it would be worthwhile to summarize my answers to the basic questions. I realize that this will probably make more people email me rather than fewer, but at least we’ll be able to skip the first few questions; the efficiency of the world should be slightly higher.

What do I do first?

Try as hard as you can to perform engineering immediately. If you can get an internship or an entry-level job, do it now. Your employer will wish that you had waited until you had more skills, but that’s their problem, not yours. Do not wait until you have completed engineering school– you might hate engineering, and school is expensive in time and money.

People hate engineering?

Yes, they do. Engineering is the impossible job, never completed, too complex to understand in full, always held up by details, perverted by marketing, unending in unexpected wrinkles. Years later, you realize your approach was wrong. Plus, it’s boring.

You should assume that engineering is not the right job for you. You have to be indoors almost all of the time. You must sit at a desk and use a computer most of the day. Most of the results of your labor will be thrown away or, if you’re lucky, recycled for the next design. You will be asked to resend spreadsheets to people who will not understand them. If you aren’t truly obsessed with solving hard problems, it’s not for you. If you cannot artificially sustain your interest in dry topics, it’s not for you.

Yeah, yeah. What if I can’t get a job right away?

There are two second-tier classes of tasks: figuring out whether engineering will satisfy you and learning basic engineering skills. If you pay attention, you can figure out whether engineering satisfies you as you develop your skills.

If you want to be a mechanical engineer, you have to learn a CAD program like Solidworks. It is extremely likely that your first job will involve a lot of Solidworks. You might end up using Proengineer or Catia rather than Solidworks, but most renewable energy companies are startups, so Solidworks, the cheapest of the lot, is most common. If you’re well-suited for mechanical engineering, you’re probably already building stuff in your basement or garage or living room; model your next project in CAD before you build it.

If you’re leaning more toward electrical engineering, you have to learn to lay out circuit boards. The most common software packages for PCB layout among startups are, I think, Eagle and Altium Designer (formerly known as Protel). Eagle is popular because the free version is legal, but Altium Designer is more powerful and pleasant to use. Lay out a simple PCB and send it to AP Circuits for fabrication. For $50-100, you’ll have your first set of PCBs.

What else do I need to learn?

Even if you want to be a mechanical engineer, you have to learn the basics of electronics. Read chapters 1-4 and 6 of Electronic Circuits and Applications by Senturia and Wedlock, even though it’s from 1975 (thanks to Max Davis for the recommendation). Other people will recommend The Art of Electronics by Horowitz and Hill; I think Senturia and Wedlock is much clearer. If you can take a class where you get to use a soldering iron, that’s a great step, before or after reading about the topic.

You also have to learn some kind of computer programming. The languages you might learn include: C, C++, Java, C#, Python, PHP, Javascript, Perl, and Ruby. I’d start with Python or Ruby, especially if you have no programming experience. If you have any interest in embedded electronics, like the brains in battlebots, you need to learn C. If you hate or fear computers, learning to use LabView is probably your best bet.

You should also learn how the internet works. Reading the chapter 1 of Richard Stevens’ TCP/IP Illustrated, Volume 1, is a good start. If that goes well, read chapters 2-4, 9, 11, 14, and 17, plus the Wikipedia page on HTTP. I’d also recommend reading C. J. Date’s Database in Depth. If you make it halfway through, you’re ahead of most database programmers.

That sounds like a lot of work.

No, it doesn’t. That sounds like what you want to do with your spare time anyway. If you don’t want to learn this stuff, go do something else. There is a lot of work in the world that needs to be done; engineering is only the optimal solution if you’re obsessed with engineering or bad at optimization. We need good dentists, scientists, nurses, welders, luthiers, lawyers, and cheesemakers as well. Especially the cheesemakers.

But what about engineering school?

You must go to engineering school. It will take 1-5 years, depending on your background and what school you go to. There is a very slight chance that you can gain enough proficiency in engineering to make it without engineering school, but it’s much more likely that you won’t have the skill or discipline to teach yourself engineering faster and more cheaply than engineering school would. You can probably learn basic mechanics of materials or some math on your own, but learning advanced topics on your own is painful and slow.

For renewable energy, make sure you pay attention during thermodynamics and heat transfer. You need to take a statistics course as well.

No, you should not get a PhD.

Wait, you’ve barely mentioned renewable energy.

Renewable energy engineering is a broad field. The major skill you need is the ability to understand the physical principles that drive complex systems. For example, when someone says to you on the street, “I have a remarkable new device that will allow this skyscraper to capture all the energy it needs from the wind or sun,” you listen politely. You do not need to know the details of the geometry of their vertical axis wind turbine or their solar concentrating lens because you know that the power density of wind and solar energy flows is 100-1000 times less than the power density of a skyscraper. When they have stopped talking, you ask them what power density they expect to reach. They stare at you blankly, and you can continue on your way, certain that you have not missed an opportunity.

Even when you are not turning down equity opportunities on the street, you should expect your engineering to run in a pattern that starts with the application of broad principles. Once you have an approach that can work in theory, you develop the details.

Beyond a broad engineering education, you need to learn the basics of solar, wind, and biofuels, plus or minus whatever you’re interested in.

I’d start with Vaclav Smil’s books:

Then I’d read in the subject areas.

Solar: Applied Photovoltaics by Wenham, Green, et al.
Wind: Wind Turbines by Erich Hau is outstanding.
Biofuels: Ahindra Nag’s book is at least decent.

So I just read books for years?

No. As I said at the outset, most important is that you start engineering right away. Build something today, and measure whether its performance is in accordance with theory. The renewable energy field is rife with crackpots; we need well-educated engineers to solve our worsening energy problems. Get to work, now.