Equation error vs simulation error

Our problem formulation inevitably begins with the data. In practice, if we have access to a physical system, instrumented using digital electronics, then we have a system in which we can apply input commands, $\bu_n$, at some discrete rate, and measure the outputs, $\by_n$ of the system at some discrete rate. We normally assume these rates our fixed, and often attempt to fit a state-space model of the form \begin{equation} \bx[n+1] = f_\balpha(\bx[n], \bu[n]), \qquad \by[n] = g_\balpha(\bx[n],\bu[n]), \label{eq:ss_model}\end{equation} where I have used $\balpha$ again here to indicate a vector of parameters. In this setting, a natural formulation is to minimize a least-squares estimation objective: $$\min_{\alpha,\bx[0]} \sum_{n=0}^{N-1} \| \by[n] - \by_n \|^2_2, \qquad \subjto \, (\ref{eq:ss_model}).$$ I have written purely deterministic models to start, but in general we expect both the state evolution and the measurements to have randomness. Sometimes, as I have written, we fit a deterministic model to the data and rely on our least-squares objective to capture the errors; more generally we will look at fitting stochastic models to the data.

We often separate the identification procedure into two parts, in which we first estimate the state $\hat{\bx}_n$ given the input-output data $\bu_n, \by_n$, and then focus on estimating the state-evolution dynamics in a second step. The dynamics estimation algorithms fall into two main categories:

• Equation error minimizes only the one-step prediction error: $$\min_{\alpha} \sum_{n=0}^{N-2} \| f_\balpha(\hat\bx_n, \bu_n) - \hat{\bx}_{n+1} \|^2_2.$$
• Simulation error captures the long-term prediction error: $$\min_{\alpha} \sum_{n=1}^{N-1} \| \bx[n] - \hat{\bx}_n \|^2_2, \qquad \subjto \quad \bx[n+1] = f_\balpha(\bx[n], \bu_n),\, \bx[0] = \hat\bx_0,$$

Within this framework of model learning, we will face all of the standard questions from supervised learning about sample efficiency and generalization.

Online optimization

Theoretical machine learning has brought a number of fresh perspectives to this old problem of system identification. Perhaps most relevant is the perspective of online optimization Rakhlin14+Hazan16+Rakhlin22. In addition to changing the formulation to considering streaming online data (as opposed to offline, batch analysis), this field has (re)popularized the notion of an old idea from decision theory: using regret as the loss function for identification.

In the regret formulation, we acknowledge that the class of parameterized models that we are searching over can likely not perfectly represent the input-output data. Rather than score ourselves directly on the prediction error, we can score ourselves relative to the predictions that would be made by best parameters in our model class. For equation-error formulations, the regret $R$ at step $N$ would take the form $$R[N] = \sum_{n=0}^{N-2} \| f_\balpha(\hat\bx_n, \bu_n) - \hat{\bx}_{n+1} \|^2_2 - \min_{\alpha^*} \sum_{n=0}^{N-2} \| f_{\alpha^*}(\hat\bx_n, \bu_n) - \hat{\bx}_{n+1} \|^2_2.$$ The goal is to produce an online learning algorithm to select $\alpha$ at each step in order to drive the regret to zero. This turns out to be an important distinction, as we will see below.

Learning models for control

So far we've described the problem purely as capturing the input-output behavior of the system. If our goal is to use this model for control, then that might not be quite the right objective for a number of reasons.

First of, predicting the observations might be more than we need for making optimal decisions; Zhang20a tells that story nicely. Imagine you are trying to swing-up and balance a cart-pole system, where the only sensor is a camera. The background images are irrelevant; we know that only the state of the cart-pole should matter. If there was a movie playing on a screen in the background, you don't want to spend all of the representational power of your model predicting the next frame of the movie, at the cost of not predicting as well the next state of the cart-pole. These challenges have led to a nice active area of research on learning task-relevant models/representations; I'll devote a section to it below after we cover the basics.

There are other considerations, as well. To be useful for control, we'll need the state of our model to be observable (from online observations) and controllable. There may also be trade-offs in model complexity. If we learn a just linear model of the dynamics then we have very powerful control design tools available, but we might not capture the rich dynamics of the world. If we learn too complex of a model, then our online planning and control techniques might become a bottleneck. Watter15 makes that point nicely. At the time of this writing, I think it's broadly fair to say that deep learning models are able to describe incredibly rich dynamics, but that our control tools for making decisions with these models are still painfully weak.

Kinematic parameters and calibration

We can separate the parameters in the multibody equations again into kinematic parameters and dynamic parameters. The kinematic parameters, like link lengths, describe the coordinate transformation from one joint to another joint in the kinematic tree. It is certainly possible to write an optimization procedure to calibrate these parameters; you can find a fairly thorough discussion in e.g. Chapter 11 of Khalil04. But I guess I'm generally of the opinion that if you don't have accurate estimates of your link lengths, then you should probably invest in a tape measure before you invest in nonlinear optimization.

One notable exception to this is calibration with respect to joint offsets. This one can be a real nuisance in practice. Joint sensors can slip, and some robots even use relative rotary encoders, and rely on driving the joint to some known hard joint limit each time the robot is powered on in order to obtain the offset. I've worked on one advanced robot that had a quite elaborate and painful kinematic calibration procedure which involve fitting additional hardware over the joints to ensure they were in a known location and then running a script. Having an expensive and/or unreliable calibration procedure can put a damper on any robotics project.

The general approach to estimating joint offsets from data is to write the equations with the joint offset as a parameter, e.g. for the double pendulum we would write the forward kinematics as: $${\bf p}_1 =l_1\begin{bmatrix} \sin(\theta_1 + \bar\theta_1) \\ - \cos(\theta_1 + \bar\theta_1) \end{bmatrix},\quad {\bf p}_2 = {\bf p}_1 + l_2\begin{bmatrix} \sin(\theta_1 + \bar\theta_1 + \theta_2 + \bar\theta_2) \\ - \cos(\theta_1 + \bar\theta_1 + \theta_2 + \bar\theta_2) \end{bmatrix}.$$ We try to obtain independent measurements of the end-effector position (e.g. from motion capture, from perhaps some robot-mounted cameras, or from some mechanical calibration rig) with their corresponding joint measurements, to obtain data points of the form $ \langle {\bf p}_2, \theta_1, \theta_2 \rangle$. Then we can solve a small nonlinear optimization problem to estimate the joint offsets to minimize a least-squares residual.

If independent measurements of the kinematics are not available, it it possible to estimate the offsets along with the dynamic parameters, using the trigonometric identities, e.g. $s_{\theta + \bar\theta} = s_\theta c_\bar\theta + c_\theta s_\bar\theta,$ and then including the $s_\bar\theta, c_\bar\theta$ terms (separately) in the "lumped parameters" we discuss below.

Estimating inertial parameters (and friction)

Now let's think about estimating the inertial parameters of multibody system. We've been writing the manipulation equations in the form: \begin{equation}\bM({\bq})\ddot{\bq} + \bC(\bq,\dot{\bq})\dot\bq = \btau_g(\bq) + \bB\bu + \text{friction, etc.}\end{equation} Each of the terms in this equation can depend on the parameters $\balpha$ that we're trying to estimate. But the parameters enter the multibody equations in a particular structured way: the inverse dynamics equations are affine in the inertial parameters Atkeson86. Notably, the forward dynamics (solving for $\ddot{\bq}$) involves a ${\bf M}^{-1}(q)$; these dynamics are not affine in the parameters. Generic system identification approaches which treat the multibody dynamics as $\bx[n+1] = f(\bx[n], \bu[n])$ will not be leveraging this convexity directly.

To proceed with the inverse dynamics formulation, let us first choose a parameterization for the inertias. Wensing17a shows that the natural parameterization of the inertial terms for a particular rigid body, $B$, are in terms of the density-weighted moment matrices: \begin{gather*} m_B &=& \iiint_{V_B} \rho(\bx)dV \in \Re \\ {\bf h}_B &=& \iiint_{V_B} \bx \rho(\bx) dV \in \Re^3 \\ {\bf \Sigma}_B &=& \iiint_{V_B} \bx \bx^T \rho(\bx)dV \in S_3, \end{gather*} where $V_B$ is the volume of the geometry in body $B$, the integral notation is for the 3D volume integral, and $S_3$ is the group of 3x3 (real) symmetric matrices. $m_B$ is the total mass of the body, ${\bf h}_B$ is simply the mass times the position of the center of mass, and ${\bf \Sigma}_B$ is closely related to the more standard form of the inertia matrix. With the integral using $x$ represented in the body frame, $B$, we can extract ${}^B{\bf p}^{B_{com}} = \frac{1}{m}{\bf h}_B$ and ${\bf I}^{B/Bo} = \tr({\bf \Sigma}_B)\mathbb{I}_{3 \times 3} - {\bf \Sigma}_B.$ Since ${\bf \Sigma}_B$ is symmetric, we can summarize all of the inertial parameters for body $B$ in 10 real numbers: $$\balpha_B = \begin{bmatrix} m_B, {\bf h}_B^T, \mathrm{vec}_{\mathrm{sym}}({\bf \Sigma}_B) \end{bmatrix} \in \Re^{10},$$ where $\mathrm{vec}_{\mathrm{sym}}$ takes the 6 upper triangular elements of the 3x3 symmetric matrices into a vector. Naturally, if we aim to estimate the parameters of many bodies simultaneously, we can compose these into a single parameter vector $\balpha.$

Remarkably, the inverse dynamics are affine in these parameters. This means that the manipulator equations above can be refactored into the form \begin{equation}{\bf W}(\bq_n,\dot{\bq}_n, \ddot{\bq}_n, \bu_n) \balpha_B + \bw_0(\bq_n,\dot{\bq}_n, \ddot{\bq}_n, \bu_n) = 0. \label{eq:least_squares}\end{equation} We sometimes refer to the ${\bf W}$ which is assembled by vertically stacking these equations evaluated at all time samples, $n$, as the "data matrix". Given data trajectories of positions, $\bq_n,$ and torque measurements, $\bu_n,$ we can estimate $\dot{\bq}_n$ and $\ddot{\bq}_n$ with (acausal) filtering. Then we can estimate the parameters using least-squares estimation. As such, we can leverage all of the strong results and variants from linear estimation. For instance, we can add terms to regularize the estimate (e.g. to stay close to an initial guess), and we can write efficient recursive estimators for optimal online estimation of the parameters using recursive least-squares. My favorite recursive least-squares algorithm uses incremental QR factorizationKaess08.

Not all values of $\alpha_B$ describe physically-valid inertias. For instance, we know that $m_B$ has to be positive. There are constraints on the inertia matrix itself. Wensing17a showed that these constraints are exactly transcribed as positive-definiteness of the "pseudo-inertia": $${\bf J}(\balpha_B) := \begin{bmatrix} {\bf \Sigma}_B & {\bf h}_B \\ {\bf h}_B^T & m_B \end{bmatrix} \succ 0.$$ If we know the geometry of the body, $B$, that we are estimating, then we can optionally write additional convex constraints on these parameters. For instance, we know that the center of mass must lie inside the convex hull of the geometry: $${}^B{\bf p}^{B_{com}} \in \mathrm{conv}(V_B) \Rightarrow {\bf h}_B \in m_B \, \mathrm{conv}(V_B).$$ This non-negative scaling of a convex set is still a convex set (known as the "perspective" of the set). Wensing17a also gives a convex constraint on the second-moment matrix based on e.g. the minimum-volume Loewner-John ellipsoid containing the known geometry.

Do these additional constraints matter? In the limit of infinite data, then one might hope that the least-squares estimator would obtain the correct values without the constraints. But their are two foils to this ideal. First, we often wish to operate in a relatively low-data regime. Second, it is important to understand that not all of the inertial parameters will be identifiable. For instance, if you imagine a robotic arm bolted to a table where the first link is a rotation about the $z$ axis. Then no amount of data obtained from this robot will given any information about the moments of inertia of the first link about the $x$ and $y$ axes. Naively including these parameters in the estimation problem could not only lead to poor estimates of the unidentifiable parameters, but can even render the physical-validity constraints on the inertial parameters vacuous.

There are two general approaches to dealing with the unidentifiable parameters. First, we can exclude them from the optimization. There are a number of algorithms for symbolically extracting the identifiable parameters, but it is more common (and likely just as effective) to simply evaluate the null-space of a sufficiently rich data matrix, $W$, to find the unidentifiable parameters by using, e.g., $\bR_1\alpha_l$ from the QR factorization: $${\bf W} = \begin{bmatrix} \bQ_1 & \bQ_2 \end{bmatrix} \begin{bmatrix} \bR_1 \\ 0 \end{bmatrix}, \quad \Rightarrow \quad {\bf W}\balpha_B = \bQ_1 (\bR_1 \balpha_B) = \bQ_1 \bar{\balpha}_{B},$$ where $\bar{\balpha}_B$ are the identifiable parameters. To implement the physical-validity constraints, we would need to use estimates of the unidentifiable parameters. Alternatively, one can formulate the optimization directly on the full parameters and add a regularization term which (ideally) operates in the nullspace of the primary objective to drive the unidentifiable parameters towards a prior estimate. Lee19 gives a thorough treatment of this approach and provides a convex formulation of a geometrically consistent regularizer, which they call the entropic divergence.

The least-squares inverse-dynamics formulation also allows one to simultaneously fit a number of other relevant dynamic parameters which enter linearly into the equations. These include reflected inertial (as a diagonal matrix times $\ddot{\bq}$), joint friction (as a constant torque multiplied by the sign of $\dot{\bq}$), and damping (as a linear multiple of $\dot{\bq}$).

Allow me to mention one subtle point about numerical scaling. All of the equations in the least-squares inverse-dynamics formulation (\ref{eq:least_squares}) have units of torque. But the torques used on the base link of a robotic arm might be numerically much larger than the torques corresponding to the more distal joints. Lee19 addresses this by solving the unconstrained least-squares problem once, then taking the error covariance of this solution into the objective function for the constrained least-squares optimization. This serves to make the least-squares objective coordinate independent.

Solving the least-squares optimizations described here subject to convex (semidefinite) constraints requires solving small semidefinite programs. There has been some nice work on parameterizations which lead to unconstrained, smooth, but nonconvex, optimization problems which can be solved with e.g. gradient descent Rucker22. In some cases, one can still guarantee that all local minima are global minima. But my impression from talking to the authors is that numerics of the semidefinite programming formulation can be notable better, and that it should be preferred if SDP solvers are available.

Simultaneous kinematic and inertial identification via lumped parameters.

Although it is highly recommended to separate kinematic an inertial estimation, some amount of joint estimation is possible when needed. We know that the inverse dynamics are affine in the inertial parameters. But it turns out the can also be written as affine in the "lumped parameters": nonlinear combinations of inertial and kinematic which do not depend on the data. I feel this is best explained with a few simple examples.

It is worth taking a moment to reflect on this factorization. First of all, it does represent a somewhat special statement about the multibody equations: the nonlinearities enter only in a particular way. For instance, if I had terms in the equations of the form, $\sin(m \theta)$, then I would not be able to produce an affine decomposition separating $m$ from $\theta$. Fortunately, that doesn't happen in our mechanical systems Khalil04. Again, this structure is particular to the inverse dynamics, as we have written here. If you were to write the forward dynamics, multiplying by $\bM^{-1}$ in order to solve for $\ddot{\bq}$, then once again you would destroy this affine structure.

Symbolic dynamics:
(0.10000000000000001 * v(0) - u(0) + (pow(v(1), 2) * mp * l * sin(q(1))) + (vd(0) * mc) + (vd(0) * mp) - (vd(1) * mp * l * cos(q(1))))
(0.10000000000000001 * v(1) - (vd(0) * mp * l * cos(q(1))) + (vd(1) * mp * pow(l, 2)) + 9.8100000000000005 * (mp * l * sin(q(1))))
        

The existence of the lumped-parameter decomposition reveals that the equation error for lumped-parameter estimation, with the error taken in the torque coordinates, can again be solved using least squares. Importantly, because we have reduced this to a least-squares problem, we can also understand when it will not work. Once again, not all of the lumped-parameters are identifiable, so we will typically need a second step to operate on the identifiable lumped parameters or to add a regularization term.

Estimated Lumped Parameters:
(mc + mp).  URDF: 11.0,  Estimated: 10.905425349337081
(mp * l).  URDF: 0.5,  Estimated: 0.5945797067753813
(mp * pow(l, 2)).  URDF: 0.25,  Estimated: 0.302915745122919

Should we be happy with only estimating the (identifiable) lumped parameters? Isn't it the true original parameters that we are after? The linear algebra decomposition of the data matrix (assuming we apply it to a sufficiently rich set of data), is actually revealing something fundamental for us about our system dynamics. Rather than feel disappointed that we cannot accurately estimate some of the parameters, we should embrace that we don't need to estimate those parameters for any of the dynamic reasoning we will do about our equations (simulation, verification, control design, etc). The identifiable lumped parameters are precisely the subset of the lumped parameters that we need.

Identification using energy instead of inverse dynamics.

In addition to leveraging tools from linear algebra, there are a number of other refinements to the basic recipe that leverage our understanding of mechanics. One important example is the "energy formulations" of parameter estimation Gautier97.

We have already observed that evaluating the equation error in torque space (inverse dynamics) is likely better than evaluating it in state space space (forward dynamics), because we can avoid inverting the mass matrix and stay linear in the lumped parameters. The total energy of the system (kinetic + potential) is also affine in the lumped parameters. We can use the relation $$\dot{E}(\bq,\dot\bq,\ddot\bq) = \dot\bq^T (\bB\bu + \text{ friction, ...}).$$

Why might we prefer to work in energy coordinates rather than torque? The differences are apparent in the detailed numerics. In the torque formulation, we find ourselves using $\ddot\bq$ directly. Conventional wisdom holds that joint sensors can reliably report joint positions and velocities, but that joint accelerations, often obtained by numerically differentiating twice, tend to be much more noisy. In some cases, it might be numerically better to apply finite differences to the total energy of the system rather than to the individual joints ^†. ^† Sometimes these methods are written as the numerical integration of the power input, rather than the differentiation of the total energy, but the numerics should be no different. Gautier96 provides a study of various filtering formulations which leads to a recommendation in Gautier97 that the energy formulation tends to be numerically superior. This distinction is probably more important in online parameter estimation, where the filters must be causal, as opposed to offline parameter estimation where acausal filters can be used.

Residual physics models with linear function approximators

The term "residual physics" has become quite popular recently (e.g. Zeng20) as people are looking for ways to combine the modeling power of deep neural networks with our best tools from mechanics. But actually this idea is quite old, and there is a natural class of residual models that fit very nicely into our least-squares lumped-parameter estimation. Specifically, we can consider models of the form: \begin{equation}\bM({\bq})\ddot{\bq} + \bC(\bq,\dot{\bq})\dot\bq = \btau_g(\bq) + \balpha_r {\bPhi}(\bq,\dot\bq) + \bB\bu + \text{friction, etc.},\end{equation} with $\balpha_r$ the additional parameters of the residual and $\bPhi$ a set of (fixed) nonlinear basis functions. The hope is that these residual models can capture any "slop terms" in the dynamics that are predictable, but which we did not include in our original parametric model. Nonlinear friction and aerodynamic drag are commonly cited examples.

Common choices for these basis functions, $\bPhi(\bq,\dot\bq)$, for use with the manipulator equations include radial basis functionsSanner91 or wavelets Sanner98. Although allowing the basis functions to depend on $\ddot\bq$ or $\bu$ would be fine for the parameter estimation, we tend to restrict them to $\bq$ and $\dot\bq$ to maintain some of the other nice properties of the manipulator equations (e.g. control-affine).

Due to the maturity of least-squares estimation, it is also possible to use least-squares to effectively determine a subset of basis functions that effectively describe the dynamics of your system. For example, in Hoburg09a, we applied least-squares to a wide range of physically-inspired basis functions in order to make a better ODE approximation of the post-stall aerodynamics during perching, and ultimately discarded all but the small number of basis functions that best described the data. Nowadays, we could apply algorithms like LASSO for least-squares regression with an $\ell_1$-regularization, or Brunton16a uses an alternative based on sequential thresholded least-squares.

Experiment design as a trajectory optimization

One key assumption for any claims about our parameter estimation algorithms recovering the true identifiable lumped parameters is that the data set was sufficiently rich; that the trajectories were "parametrically exciting". Basically we need to assume that the trajectories produced motion so that the data matrix, ${\bf W}$ contains information about all of the identifiable lumped parameters. Thanks to our linear parameterization, we can evaluate this via numerical linear algebra on ${\bf W}$.

Moreover, if we have an opportunity to change the trajectories that the robot executes when collecting the data for parameter estimation, then we can design trajectories which try to maximally excite the parameters, and produce a numerically-well conditioned least squares problem. One natural choice is to minimize the condition number of the data matrix, ${\bf W}$. The condition number of a matrix is the ratio of the largest to smallest singular values $\frac{\sigma_{max}({\bf W})}{\sigma_{min}({\bf W})}$. The condition number is always greater than one, by definition, and the lower the value the better the condition. The condition number of the data matrix is a nonlinear function of the data taken over the entire trajectory, but it can still be optimized in a nonlinear trajectory optimization (Khalil04, § 12.3.4).

Online estimation and adaptive control

The field of adaptive control is a huge and rich literature; many books have been written on the topic (e.g Åström13). Allow me to make a few quick references to that literature here.

I have mostly discussed parameter estimation so far as an offline procedure, but one important approach to adaptive control is to perform online estimation of the parameters as the controller executes. Since our estimation objective can be linear in the (lumped) parameters, this often amounts to the recursive least-squares estimation. To properly analyze a system like this, we can think of the system as having an augmented state space, $\bar\bx = \begin{bmatrix} \bx \\ \balpha \end{bmatrix},$ and study the closed-loop dynamics of the state and parameter evolution jointly. Some of the strongest results in adaptive control for robotic manipulators are confined to fully-actuated manipulators, but for instance Moore14 gives a nice example of analyzing an adaptive controller for underactuated systems using many of the tools that we have been developing in these notes.

As I said, adaptive control is a rich subject. One of the biggest lessons from that field, however, is that one may not need to achieve convergence to the true (lumped) parameters in order to achieve a task. Many of the classic results in adaptive control make guarantees about task execution but explicitly do not require/guarantee convergence in the parameters.

Identification with contact

Can we apply these same techniques to e.g. walking robots that are making and breaking contact with the environment?

There is certainly a version of this problem that works immediately: if we know the contact Jacobians and have measurements of the contact forces, then we can add these terms directly into the manipulator equations and continued with the least-squares estimation of the lumped parameters, even including frictional parameters.

One can also study cases where the contact forces are not measured directly. For instance, Fazeli17a studies the extreme case of identifiability of the inertial parameters of a passive object with and without explicit contact force measurements.

From state observations

Let's start our treatment with the easy case: fitting a linear model from direct (potentially noisy) measurements of the state. Fitting a discrete-time model, $\bx[n+1] = \bA\bx[n] + \bB\bu[n] + \bw[n]$, to sampled data $(\bu_n,\bx_n = \bx[n]+\bv[n])$ using the equation-error objective is just another linear least-squares problem. Typically, we form some data matrices and write our least-squares problem as: \begin{gather*} {\bf X}' \approx \begin{bmatrix} \bA & \bB \end{bmatrix} \begin{bmatrix} {\bf X} \\ {\bf U} \end{bmatrix}, \qquad \text{where} \\ {\bf X} = \begin{bmatrix} \mid & \mid & & \mid \\ \bx_0 & \bx_1 & \cdots & \bx_{N-2} \\ \mid & \mid & & \mid \end{bmatrix}, \quad {\bf X}' = \begin{bmatrix} \mid & \mid & & \mid \\ \bx_1 & \bx_2 & \cdots & \bx_{N-1} \\ \mid & \mid & & \mid \end{bmatrix}, \quad {\bf U} = \begin{bmatrix} \mid & \mid & & \mid \\ \bu_0 & \bu_1 & \cdots & \bu_{N-2} \\ \mid & \mid & & \mid \end{bmatrix}.\end{gather*} By the virtues of linear least squares, this estimator is unbiased with respect to the uncorrelated process and/or measurement noise, $\bw[n]$ and $\bv[n]$.

Model-based Iterative Learning Control (ILC)

Local time-varying linear model along a trajectory + iLQR. Bregmann ADMM

The Hydrodynamic Cart-Pole

One of my favorite examples of model-based ILC was a series of experiments in which we explored the dynamics of a "hydrodynamic cart-pole" system. Think of it as a cross between the classic cart-pole system and a fighter jet (perhaps a little closer to the cartpole)!

Here we've replaced the pole with an airfoil (hydrofoil), turned the entire system on its side, and dropped it into a water tunnel. Rather than swing-up and balance the pole against gravity, the task is to balance the foil in its unstable configuration against the hydrodynamic forces. These forces are the result of unsteady fluid-body interactions; unlike the classic cart-pole system, this time we do not have an tractable parameterized ODE model for the system. It's a perfect problem for system identification and ILC.

In a series of experiments, first we attempted to stabilize the system using an LQR controller derived with an approximate model (using flat-plate theory). This controller didn't perform well, but was just good enough to collect relevant data in the vicinity of the unstable fixed point. Then we fit a linear model, recomputed the LQR controller using the model, and got notably better performance.

To investigate more aggressive maneuvers we considered making a rapid step change in the desired position of the cart (like a fighter jet rapidly changing altitude). Using only the time-invariant LQR balancing controller with a shifted set point, we naturally observed a very slow step-response. Using trajectory optimization on the time-invariant linear model used in balancing, we could do much better. But we achieved considerably better performance by iteratively fitting a time-varying linear model in the neighborhood of this trajectory and performing model-based ILC.

These experiments were quite beautiful and very visual. They went on from here to consider the effect of stabilizing against incoming vortex disturbances using real-time perception of the oncoming fluid. If you're at all interested, I would encourage you to check out John Robert's thesisRoberts12 and/or even the video of his thesis defense.

Consider adding a video or two here? John just shared his thesis slides with me in dropbox, and I found some of the original videos in movies/RobotLocomotion/WaterTunnel. But the videos have big black borders, so could really use some cleanup. How reasonable (locally) are LQG models through contact?

From input-output data (the state-realization problem)

In the more general form, we would like to estimate a model of the form \begin{gather*} \bx[n+1] = \bA \bx[n] + \bB \bu[n] + \bw[n]\\ \by[n] = \bC \bx[n] + {\bf D} \bu[n] + \bv[n]. \end{gather*} Once again, we will apply least-squares estimation, but combine this with the famous "Ho-Kalman" algorithm (also known as the "Eigen System Realization" (ERA) algorithm) Ho66. My favorite presentation of this algorithm is Oymak19.

First, observe that \begin{align*} \by[0] =& \bC\bx[0] + {\bf D}\bu[0] + \bv[0], \\ \by[1] =& \bC(\bA\bx[0] + \bB\bu[0] + \bw[0]) + {\bf D}\bu[1] + \bv[1], \\ \by[n] =& \bC\bA^n\bx[0] + {\bf D}\bu[n] + \bv[n] + \sum_{k=0}^{n-1}\bC\bA^{n-k-1}(\bB\bu[k] + \bw[k]). \end{align*} For the purposes of identification, let's write $\by[n]$ as a function of the most recent $N+1$ inputs (for $k \ge N$): \begin{align*}\by[n] =& \begin{bmatrix} \bC\bA^{N-1}\bB & \bC\bA^{N-2}\bB & \cdots & \bC\bB & {\bf D}\end{bmatrix} \begin{bmatrix} \bu[n-N] \\ \bu[n-N+1] \\ \vdots \\ \bu[n] \end{bmatrix} + {\bf \delta}[n] \\ =& {\bf G}[n]\bar\bu[n] + {\bf \delta}[n] \end{align*} where ${\bf \delta}[n]$ captures the remaining terms from initial conditions, noise and control inputs before the (sliding) window. ${\bf G}[n] = \begin{bmatrix} \bC\bA^{N-1}\bB & \bC\bA^{N-2}\bB & \cdots & \bC\bB & {\bf D}\end{bmatrix},$ and $\bar\bu[n]$ represents the concatenated $\bu[n]$'s from time $n-N$ up to $n$. Importantly we have that ${\bf \delta}[n]$ is uncorrelated with $\bar\bu[n]$: $\forall k > n$ we have $E\left[\bu_k {\bf \delta}_n \right] = 0.$ This is sometimes known as Neyman orthogonality, and it implies that we can estimate ${\bf G}$ using simple least-squares $\hat{\bf G} = \argmin_{\bf G} \sum_{n \ge N} \| \by_n - {\bf G}\bar\bu_n \|^2.$ Oymak19 gives bounds on the norm of the estimation error as a function of the number of samples and the variance of the noise.

How should we pick the window size $N$? By Neyman orthogonality, we know that our estimates will be unbiased for any choice of $N \ge 0$. But if we choose $N$ too small, then the term $\delta[n]$ will be large, leading to a potentially large variance in our estimate. For stable systems, the $\delta$ terms will get smaller as we increase $N$. In practice, we choose $N$ based on the characteristic settling time in the data (roughly until the impulse response becomes sufficiently small).

If you've studied linear systems, ${\bf G}$ will look familiar; it is precisely this (multi-input, multi-output) matrix impulse response, also known as the "Markov parameters". In fact, estimating $\hat{\bf G}$ may even be sufficient for control design, as it is closely-related to the parameterization used in disturbance-based feedback for partially-observable systems Sadraddini20+Simchowitz20. But the Ho-Kalman algorithm can be used to extract good estimates $\hat\bA, \hat\bB, \hat\bC, \hat{\bf D}$ with state-dimension $\dim(\bx) = n_x$ from $\hat{\bf G},$ assuming that the true system is observable and controllable with order at least $n_x$ and the data matrices we form below are sufficiently richOymak19.

It is important to realize that many system matrices, $\bA, \bB, \bC,$ can describe the same input-output data. In particular, for any invertible matrix (aka similarity transform), ${\bf T}$, the system matrices $\bA, \bB, \bC$ and $\bA', \bB', \bC'$ with, $$\bA' = {\bf T}^{-1} \bA {\bf T},\quad \bB' = {\bf T}^{-1} \bB,\quad \bC' = \bC{\bf T},$$ describe the same input-output behavior. The Ho-Kalman algorithm returns a balanced realization. A balanced realization is one in which we have effectively applied a similarity transform, ${\bf T},$ which makes the controllability and observability Grammians equal and diagonalBrunton19 and which orders the states in terms of diminishing effect on the input/output behavior. This ordering is relevant for determining the system order and for model reduction.

Note that $\hat{\bf D}$ is the last block in $\hat{\bf G}$ so is extracted trivially. The Ho-Kalman algorithm tells us how to extract $\hat\bA, \hat\bB, \hat\bC$, with another application of the SVD on suitably constructed data matrices (see e.g. Oymak19, Juang01, or Brunton19).

Ho-Kalman identification of the cart-pole from keypoints

Let's repeat the cart-pole example. But this time, instead of getting observations of the joint position and velocities directly, we will consider the problem of identifying the dynamics from a camera rendering the scene, but let us proceed thoughtfully.

The output of a standard camera is an RGB image, consisting of 3 x width x height real values. We could certainly columnate these into a vector and use this as the outputs/observations, $y(t)$. But I don't recommend it. This "pixel-space" is not a nice space. For example, you could easily find a pixel that has the color of the cart in one frame, but after an incremental change in the position of the cart, now (discontinuously) takes on the color of the background. Deep learning has given us fantastic new tools for transforming raw RGB images into better "feature spaces", that will be robust enough to deploy on real systems but can make our modeling efforts much more tractable. My group has made heavy use of "keypoint networks" Manuelli19 and self-supervised "dense correspondences" Florence18a+Florence20+Manuelli20a to convert RGB outputs into a more consumable form.

The existence of these tools makes it reasonable to assume that we have observations representing the 2D positions of some number of "key points" that are rigidly fixed to the cart and to the pole. The location of any keypoint $K$ on a pole, for instance, is given by the cart-pole kinematics: $${}^{W}{\bf p}^{K} = \begin{bmatrix} x \\ 0 \end{bmatrix} + \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} {}^{P}{\bf p}^{K},$$ where $x$ is the position of the cart, and $\theta$ is the angle of the pole. I've used Monogram notation to denote that ${}^{P}{\bf p}^{K}$ is the Cartesian position of a point $A$ in the pole's frame, $P$, and ${}^{W}{\bf p}^{K}$ is that same position in the camera/world frame, $W$. While it is probably unreasonable to linearize the RGB outputs as a function of the cart-pole positions, $x$ and $\theta$, it is reasonable to approximate the keypoint positions by linearizing this equation for small angles.

You'll see in the example that we recover the impulse response quite nicely.

The big question is: does Ho-Kalman discover the 4-dimensional state vector that we know and like for the cart-pole? Does it find something different? We would typically choose the order by looking at the magnitude of the singular values of the Hankel data matrix. In this example, in the case of no noise, you see clearly that 2 states describe most of the behavior, and we have diminishing returns after 4 or 5:

As an exercise, try adding process and/or measurement noise to the system that generated the data, and see how they effect this result.

Adding stability constraints

Details coming soon. See, for instance Umenberger18.

Autoregressive models

Another important class of linear models predict the output directly from a history of recent inputs and outputs: \begin{align*} \by[n+1] =&\bA_0 \by[n] + \bA_1 \by[n-1] + ... + \bA_k \by[n-k] \\ & + \bB_0\bu[n] + \bB_1 \bu[n-1] + ... \bB_k \bu[n-k]\end{align*} These are the so-called "AutoRegressive models with eXogenous input (ARX)" models. The coefficients of ARX models can be fit directly from input-output data using linear least squares.

While it is certainly possible to think of these models as state-space models, where we collect the finite history $\by$ and $\bu$ (except not $\bu[n]$) into our state vector. But this is not necessarily a very efficient representation; the Ho-Kalman algorithm might be able to find a much smaller state representation. This is particularly important for partially-observable systems that might require a very long history, $k$.

Statistical analysis of learning linear models

From state observations

Following analogously to our discussion on linear dynamical systems, the first case to understand is when we have direct access to state and action trajectories: $s[\cdot], a[\cdot]$. This corresponds to identifying a Markov chain or Markov decision process (MDP). The transition matrices can be extracted directly from the statistics of the transitions: $$T_{ijk} = \frac{\text{number of transitions to } s_i\text{ from }s_j, a_k}{\text{total number of transitions from }s_j, a_k}.$$ Not surprisingly, for this estimate to converge to the true transition probabilities asymptotically, we need the identification (exploration) dynamics to be ergodic (for every state/action pair to be visited infinitely often).

For large MDPs, analogous to the modal decomposition that we described for linear dynamical systems, we can also consider factored representations of the transition matrix. [Coming soon...]

Identifying Hidden Markov Models (HMMs)

Generating training data

One extremely interesting question that arises when fitting rich nonlinear models like neural networks is the question of how to generate the training data. You might have the general sense that we would like data that provides some amount of "coverage" of the state space; this is particularly challenging for underactuated systems since we have dynamic constraints restricting our trajectories in state space.

For multibody parameter estimation we discussed using the condition number of the data matrix as an explicit objective for experiment design. This was a luxury of using linear least-squares optimization for the regression and it takes advantage of the specialized knowledge we have from the manipulator equations about the model structure. Unfortunately, we don't have an analogous concept in the more general case of nonlinear least squares with general function approximators. This approach also works for identifying linear dynamical systems. The challenge is here perhaps considered less severe since for linear models in the noise-free case even data generated from the impulse response is sufficient for identification.

This topic has received considerable attention lately in the context of model-based reinforcement learning (e.g. Agarwal20b). Broadly speaking, in the ideal case we would like the training data we use for system identification to match the distribution of data that we will encounter when executing the optimal policy. But one cannot know this distribution until we've designed our controller, which (in our current discussion) requires us to have a model. It is a classic "chicken and the egg" problem. In most cases, it speaks to the importance to interleaving system identification and control design instead of the simpler notion of performing identification once and then using the model forevermore.

From state observations

State-space models from input-output data (recurrent networks)

Input-output (autoregressive) models

You might have heard about a neural network model called the Transformercite attention is all you need. There are a massive number of online resources where you can learn about Transformers. These started in the natural language processing (NLP) community, where they rapidly replaced recurrent neural network models which were, at the time, struggling to handle very long sentences (trajectories). Transformers get around this problem by switching back to an autoregressive-style model, but they compress the very long history of previous observations into a compact representation using attention mechanisms. This is what now enables large language models, like chat-GPT, to e.g. write an entire essay, "remembering" names/concepts that it used early in the essay correctly very late in the essay.

Particle-based models

Object-centric models

Modeling stochasticity

Control design for neural network models