Designing embedded systems with web frameworks

I have a prediction.

We're about to see a shift in embedded systems development. Around 2008 or 2009, embedded microprocessors like the one in your cellphone, reached a threshold where they can perform as decent webservers without special tuning. Even with a slow database query or some inefficient templates, they've got the speed to handle real web serving. This means that suddenly the most efficient development pattern for embedded systems is not the proprietary hardware and software tools that have dominated the industry for the last 30 years, but the open source web development tools that have arisen in the last decade.

The transition will be painfully slow, and of course there are some niches where specialized hardware and real-time control will preclude the use of generic web tools-- motion control springs to mind. But I think that the combination of cheap hardware and modern web frameworks will crush the industrial controller market, like digital cameras did to film cameras.

First, some background on tiny computers

There are millions of tiny computers used for monitoring and control systems around the world. Let's break them into two categories: small (microcontrollers like the Arduino or a PIC development board, which runs $10-500) and large (embedded controllers like National Instruments hardware, which cost $500-5000).

To use the small ones, you write code from scratch, perhaps with some pre-written libraries to talk to certain peripherals and a bootloader to run your code on power-up. The vast majority cannot be connected to the internet without substantial effort, and when connected to the internet, they aren't powerful enough to work like most servers on the internet. For example, a webserver built on a small microcontroller would be overwhelmed by the background noise on the internet, i.e. traffic from bots and viruses. This kind of system is perfect if you want to log temperature in your basement, or turn on a light whenever the garage door is open. They're no good for running an inventory management system in a warehouse with 20 workers.

Below: an Arduino

The large ones come with an operating system, like Windows CE, Linux, or VxWorks. Most of the devices are reworked versions of hardware from the pre-internet days that have had Ethernet ports added to them. They can handle real internet traffic, but they usually use proprietary software to do it. They're shaped like a long, skinny shoebox.

Below: a larger controller

The change is that something equivalent to the iPhone, which uses a 600 MHz ARM processor, can store years of data and serve it up to anyone with a web browser in seconds with hardware cost of a few hundred dollars. Industrial controllers lose on cost of both hardware and programming time, and performance. The hardware has to be expensive to support the R&D costs, which are borne entirely by each manufacturer. There's just no way that even a large industrial control company, which might have a dozen dedicated programmers at best, can compete with the thousands of developers working on open source web software worldwide.

How web software development has changed

In the 90's, software development for the internet meant either writing server software or designing static web pages. Starting around 2000 (give or take a few years), websites started incorporating dynamic data, which was stored in databases. Around 2005, a new kind of web software gained popularity-- the web framework-- with the release of Ruby on Rails.

With old-style web development, a web programmer would write code that inserted data into a database, more code that updated the database with new data, and more code that retrieved and sorted the data for presentation in a web page. With a web framework, the programmer writes out a template for how the data should be presented on a webpage, and the framework figures out what to request from the database. Web frameworks can't scale to the level of a big website like Amazon, but for low traffic systems, they work fine and reduce the programming time needed dramatically.

Most of the time, what people want to do with microcontrollers is log some data from sensors and maybe trigger some actuators in response. After they log the data, people want to analyze the data, make graphs with it, and then do it again, maybe with a different sensor. This matches well with the typical database-backed website. The only substantial additions are code to interact with hardware-- read sensors and trigger actuators. I think this is the least developed part of the systems I expect we'll see in the next few years.

Proof of concept

It's definitely possible that I'm just some blowhard on the internet. I mean, I'm at least some blowhard on the internet, but I might still be right in predicting this change. To that end, I've done some testing to see whether I'm going in the right direction.

Using an off-the-shelf microcontroller kit that cost $339 plus shipping, I installed a Python web framework called Django and wrote code to make it act like a thermostat to replace the one in our house. It took about 2.5 weekends to write the code, which is much faster than I could write a similar application for something like an Arduino with an Ethernet carrier board, and this was my first try. I had played around with Django a few times previously, but this was my largest effort by far.

Below: the proof of concept

The hardware interface is just a cron job that runs once per minute. It queries a temperature sensor on the PCB using a short C program called by a Python script, which also logs the data to a SQLite database. Retrieval of a webpage that queries the database for the last 300 datapoints (5 hours worth) and generates a chart of the data using Javascript takes around 1.5 seconds. That's with a processor running at 250 MHz (relatively slow) and the database stored on an SD card. Most of the time is spent converting Python datetimes to Javascript timestamps.

I could have written a similar application even faster using a GUI tool like LabView, but that requires specialized hardware that cost 3-10 times more-- either a dedicated PC with a USB device for sensors, or a industrial controller with a sensor module. With Django, I got an administrative interface, allowing different users and groups different levels of access, for free. If I want a central repository of users with LabView, I'm looking at another $800 for the "Datalogging and Supervisory Control Module." If I want to talk to a SQL database, I'll need the "Database Connectivity Toolkit" for another $1000. This is on top of the $1500 I would have paid for LabView already, plus the hardware.

Embedded control systems running web services are still an immature technology at best, but I think they'll grow up quickly in the next few years.

Why are you writing this?

I've been thinking about this change for a year or two now. I'm working on developing something like the Arduino, but a little heavier duty, so it can run a web framework, but with the hardware drivers pre-integrated. Send me an email at brandon at pingswept.org or post a comment if you want to know when it exists for real.

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.

Calculating solar panel shading in Python

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.

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.

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

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

Here's the finite difference code.

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

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

Here's the full code.

# 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

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!