Fit Params to a run¶
In this tutorial we find the mass, the stiffness and the damping of a robot from a log of its motion, as one does when a model has to match a real robot.
The run¶
The robot is a mass on a spring and a damper that slides along \(x\). A controller moves it with a virtual spring to a goal that wanders, which gives motions at several speeds, and we log the run. This robot is the truth, and the fit does not see its numbers:
import casadi as ca
import numpy as np
import matplotlib.pyplot as plt
import virtualmodelcontrol as vmc
def slider(mass, stiffness, damping):
line = vmc.models.JointSpace(1, unit="m")
robot = vmc.Mechanism("slider", model=line)
x = robot.joint(0)
robot.add("m", vmc.Inertance(x, mass)) # [kg]
robot.add("spring", vmc.LinearSpring(x, stiffness)) # [N/m]
robot.add("damper", vmc.LinearDamper(x, damping)) # [N·s/m]
return robot
truth = slider(2.0, 40.0, 1.5)
goal = vmc.Custom(lambda t: 0.1 * ca.sin(2 * t) + 0.05 * ca.sin(9 * t),
[vmc.Time()], dim=1, unit="m")
ctrl = vmc.Mechanism("ctrl")
ctrl.add("pull", vmc.LinearSpring(truth.joint(0) - goal, 60.0))
ctrl.add("damp", vmc.LinearDamper(truth.joint(0), 3.0))
controller = vmc.VMCController(
vmc.compile(vmc.VirtualMechanismSystem(truth, ctrl)))
plant = vmc.sim.ModelPlant(truth, max_step=1e-3)
log = vmc.sim.run(plant, controller, vmc.sim.SimClock(1 / 200), T=10.0)
The fit¶
We start from a model that has the right structure and the wrong numbers, and name the Params
to find. A Param is named after its component and its own name: the mass of the Inertance added as
m is m.inertance, the stiffness of spring is spring.stiffness (robot.params lists them all; see
Parameters). fit_params takes the robot, the names (or globs) of its Params and the runs: each a log from
vmc.sim.run, or a dict of arrays with one row per sample, t, q, v and motor_torque, and
a too if you have the accelerations. It asks
which values make the robot’s own dynamics, at the logged motion, agree with the torques that
were sent. It returns the fitted values and their standard errors std, each by name, and the
root-mean-square rms of what is left over:
from virtualmodelcontrol.identification import fit_params
guess = slider(1.0, 10.0, 0.2)
names = ["m.inertance", "spring.stiffness",
"damper.damping"] # component.Param
fit = fit_params(guess, names, [log], smoothing=11)
for name in names:
value, std = float(fit.values[name]), float(fit.std[name])
print(f"{name:18s} {value:8.4f} ± {std:.1e}")
m.inertance 2.0009 ± 3.3e-05
spring.stiffness 39.9983 ± 1.4e-03
damper.damping 1.5191 ± 2.2e-04
The true values are 2 kg, 40 N/m and 1.5 N·s/m. The log holds no accelerations, so the fit
gets them by smoothing the velocity and differentiating it, which is where its error comes from
(smoothing is the number of samples of the smoothing). A run that has an a of its own, as a
simulation can give, is used as it is. After the fit, 0.0027 N is left over
per sample. The numbers after ± are the standard errors: how far the values may be off if what is
left over is noise. Here it is not noise, it is the smoothing, so the errors are much larger than
the ± says (the damping is 1.3 % off, 85 times its
standard error). Read the ± to compare how well the run determines one Param against another,
not as a bound on the error.
What it asks of the data¶
A Param shows only in motion that depends on it. The mass shows in the accelerations, the
stiffness in the position, the damping in the velocity, so a run that only sits still fits
none of them: every Param keeps the value it started from, with a standard error of inf. A Param of a nonlinear component,
such as the saturating damper TanhDamper, is fitted the same way, by iterating from the values
the robot has, and the fit keeps every Param within its bounds.
To use a result, put its values into the robot: guess.params[name].value = fit.values[name]
for each name. The robot then simulates, plans and is controlled with the fitted numbers.
The hanging arm fits the stiffness and damping of a real kind of arm
from logged runs, and Tuning and Optimizing a virtual mechanism use
such a model.