Fitting a model to data

Least squares in plain Python, and a robot that drives by the model it fitted to itself.

U12.2Learning, and the capstoneUniversity35 min

Do this lesson in the simulator

A single ratio is a model with one parameter. The next step is a model with two, fitted to many measurements at once, which is where least squares earns its place.

The problem

Collect pairs: a command you sent, and the speed the robot went. Find the line

v = a * c + b

that fits them best. Best in the least squares sense means the a and b that make the sum of squared residuals as small as possible:

loss(a, b) = sum over samples of (a * c + b - v)^2

Squared, so that errors above and below both count, and so that the loss has one derivative everywhere and one minimum. Setting the two partial derivatives to zero gives two linear equations, and their solution is the formula every statistics package implements:

n = len(data)
sx = sum(c for c, v in data)
sv = sum(v for c, v in data)
sxx = sum(c * c for c, v in data)
sxv = sum(c * v for c, v in data)
slope = (n * sxv - sx * sv) / (n * sxx - sx * sx)
intercept = (sv - slope * sx) / n

Six sums and a division. No iteration, no learning rate, no local minima. When a problem can be posed as least squares on parameters that enter linearly, solve it this way and do not reach for anything cleverer.

Collecting data on a robot that moves

Data collection is the part that is actually difficult, and it is a robotics problem rather than a statistics one.

  • Let it settle. The drive has a lag of about a quarter of a second. Read the speed too soon and you have measured the ramp, not the speed. Two or three time constants is enough after a step from rest, and about five after a step that reverses the command, which is what a forwards and backwards sweep does every time. Get this wrong and every reading is low by the same fraction, which is worse than noise: it is a slope that is quietly 10 percent out.
  • Stay on the mat. A sweep from command 20 to 100 at a second each would cross the whole mat. Run every command forwards and then backwards, and the robot ends where it began.
  • Measure what you can measure. The robot does not know its true speed. It knows flow(), which is the true speed times a scale factor of its own, plus noise. So the model you fit is a model of the flow sensor's opinion, not of the world, and it inherits the flow sensor's few percent of scale error. That is honest and it is usually enough, and it is worth saying out loud in your report.
  • Average the two directions. Taking (forward - backward) / 2 cancels any fixed offset in the sensor and halves the noise, for free.
from bugbot import *
connect()

DT = 0.1

def turn_cmd():
    err = (imu()[0] + 180) % 360 - 180
    if abs(err) < 2.0:
        return 0
    cmd = max(18.0, min(30.0, abs(1.5 * err)))   # below 15 the drive does nothing at all
    return -cmd if err > 0 else cmd

def hold(cmd, ticks, collect=False):
    vs = []
    for i in range(ticks):
        drive(cmd, 0, turn_cmd())
        if collect:
            vs.append(flow()[1])
        wait(DT)
    return vs

def speed_at(cmd):
    hold(cmd, 12)                   # the command reverses between samples: five time constants, not two
    return sum(hold(cmd, 3, True)) / 3.0

data = []
for c in (20, 40, 60, 80, 100):
    out, back = speed_at(c), speed_at(-c)
    data.append((c, (out - back) / 2.0))
    print("command", c, "->", round(data[-1][1], 2), "cm/s")
stop()

Run this in the simulator

Two details in there are worth more than they look.

Hold the heading. The sweep takes half a minute, and over half a minute a robot quietly turns. Every reading after that is of a robot pointing somewhere else, and the final drive goes off at an angle. A correction of a degree or two costs nothing and removes a whole class of confusion.

Mind the dead band in the correction. A proportional heading correction of 8 percent command does absolutely nothing, because the drive ignores anything below 15. The controller looks right, the plot looks right, and the robot never turns. Give the correction a minimum magnitude above the dead band, which is U5.6 wearing different clothes.

The dead band, and what a model is for

Below about command 15 the vibration drive does not move at all. A straight line fitted across that region will have an intercept that is pure nonsense, and it will mispredict everything near the bottom of the range.

Two reasonable responses, and both are used in practice:

  1. Restrict the domain. Fit only where the model applies, and record the range it is valid over. A model without a stated domain is a trap for whoever uses it next.
  2. Change the model. v = a * max(0, c - d) has a dead band in it. It is no longer linear in d, so least squares in closed form does not apply, and you are into a search. Which is a fair trade when the dead band matters.

Using the fit

The model earns its keep when the robot acts on it. To cover a distance D at command c:

seconds = D / (a * c + b)

and the lag politely cancels: the robot loses about v * tau of travel while the motors come up, and gains about the same coasting after the stop. Drive for D / v and you arrive, without a single sensor reading during the run.

That is worth pausing on. This is open loop control done well: no feedback, no depth sensor, no camera, and it works because the model is good. A better model buys you the same accuracy with less sensing, which on a real robot is less cost, less latency and fewer things to go wrong.

Task: fit the drive

Sweep the command, fit speed against command by least squares, print slope: and intercept:, then use the fit to drive 80 cm open loop and stop there.

from bugbot import *
connect()

DT = 0.1
COMMANDS = [20, 30, 40, 50, 60, 70, 80, 90, 100]
TARGET = 80.0

Challenges

  1. Fit the lateral axis as well. Is the slope the same as forward, and should it be?
  2. Fit with the command 20 point included and excluded. How much does the intercept move?
  3. Work out, from the residuals, how far off your 80 cm drive should be expected to land, and compare with where it landed.