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

Print it: worksheet (PDF) · answers (PDF) · as a web page

Guide: Sim to real: why a policy that works in simulation fails on the robot

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. Three time constants after a step from rest leaves you 5 percent short; a step that reverses the command has twice as far to go and wants about five, and a forwards and backwards sweep reverses it 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 variance of the noise, for free.
Speed after a change of command, from rest and after a reversal, with the lesson's three readings marked00.40.81.21.6-1-0.500.51seconds after the command changesspeed, share of the new speed3 time constantsthe threereadingsfrom rest: 95%reversed: 90%reversal startsat minus full speed
With a lag of 0.25 s, three time constants leave a step from rest 5% short and a reversal 10% short. The cell holds each command for 12 ticks first, so by its readings even a reversal has reached 98%.
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.
Measured speed against command, with a line fitted across the dead band and a line fitted from command 20 up02040608010005101520commandspeed, cm/sdeadbandmeasured at 5 and 10: nothingfitted from 20 to 100:slope 0.193, intercept -0.57fitted across the dead band:slope 0.205, intercept -1.46
The two red points inside the dead band drag the dashed line down and tilt it: it now says command 20 gives 2.6 cm/s where the robot measured 3.0. Fitting only from 20 up keeps the line on the points, and 20 to 100 is the range it is valid over.

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.

Speed against time at the start and the end of an open loop drive, with the lost and gained travel shaded051015speed, cm/slost: v times tau= 3.2 cmstartgained coasting= 3.2 cmstop, 6.17 s later...dashed: the speed the model predicts, 12.96 cm/s at command 70
Driving 80 cm at command 70 means 6.17 s at the fitted 12.96 cm/s. The robot is short by about v times tau while the motors come up and makes the same back coasting after the stop, so the two shaded areas cancel.

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.