# Identification of the benchmark canonical form

This is a description of an identification procedure for the bicycle based around the benchmark canonical form.

## Model structure

The identification of the linear dynamics of the bicycle can be formulated with respect to the benchmark canonical form realized in [Meijaard2007]. $$ M \ddot{q} + v C_1 \dot{q} + [g K_0 + v^2 K_2] q = T $$ where the time varying states roll and steer are collected in the vector \( q = [\phi \quad \delta]^T \) and the time varying inputs roll torque and steer torque are collected in the vector \(T = [T_\phi \quad T_\delta]^T\). I extend the equations to properly account for the lateral perturbation force, \( F \), which was the actual input we delivered during the experiments. It contributes to both the roll torque and steer torque equations. $$ M \ddot{q} + v C_1 \dot{q} + [g K_0 + v^2 K_2] q = T + H F $$ where \(H = [H_{\phi F} \quad H_{\delta F}]^T\) is a vector describing the linear contribution of the lateral force to the roll and steer equations. \(H_{\phi F}\) is approximately the distance from the ground to the force application point. \(H_{\delta F}\) is a distance that is a function of the bicycle geometry (trail, wheelbase) and the longitudinal location of the force application point. \(H_{\delta F} H = [H_{\phi F}\). I compute \(H\) for each rider/bicycle from the state space form of the linear equations of motion. $$ \dot{x} = A x + B u $$ where \( x = [\phi \quad \delta \quad \dot{\phi} \quad \dot{\delta}]^T \) and \( u = [F \quad T_\phi \quad T_\delta]^T \). The state and input matrices can be sectioned. $$ A = \begin{bmatrix} 0 & I \\ A_l & A_r \end{bmatrix} $$ $$ B = \begin{bmatrix} 0 & 0\\ B_F & B_T \end{bmatrix} $$ where \(A_l\) and \(A_r\) are the 2 x 2 submatrices corresponding to the states and their derivatives, respectively. \(B_F\) and \(B_T\) are the 2 x 1 and 2 x 2 submatrices corresponding to the lateral force and the torques, respectively. The benchmark canonical form can now be written as $$ B_T^{-1} [ \ddot{q} - A_r \dot{q} - A_l q] = T + B_T^{-1} B_F F $$ where $$ M = B_T^{-1} $$ $$ vC_1 = -B_T^{-1} A_r $$ $$ [gK_0 + v^2K_2] = -B_T^{-1} A_l $$ $$ H = B_T^{-1} B_F $$

## Data processing

The roll and steer angle were measured with potentiometers. The roll and steer rates were measured with a combination of rate gyros. The steer torque was measured with a torque sensor and the lateral force with a load cell. The wheel speed was measured with a DC tachometer. See my draft chapter for more details of the measurements.

All of the signals were filtered with a second order low pass Butterworth filter at 15 Hz. The roll and steer accelerations were computed by numerically differentiating the roll and steer rate signals. The mean was subtracted from all the signals except the lateral force.

I explored complementary filtering various signals but found little improvement in the signal quality. The iPython notebook worksheet details what I tried.

## Identification

A simple analytic identification problem can be formulated. If we have good measurements of \(q\), their first and second derivatives, forward speed \(v\) and the inputs \(T_\delta\) \(F\), the entries in the \(M\), \(C1\), \(K0\), \(K2\) and \(H\) matrices can be identified by forming two simple regressions, one for each equation in the canonical form. I use the instantaneous speed at each time step rather than the mean over a run to improve accuracy with respect to the speed parameter. The roll and steer equation each can be put into a simple linear form $$ \Gamma \Theta = Y $$ where \( \Theta \) are the unknown parameters and \( \Gamma \) and \( \Theta \) are made up of the inputs and outputs measured during a run. \(Theta\) can generally be all or a subset of the entries in the canonical matrices. If there are \(N\) samples in a run and we desire to find \(M\) entries in the equation, then \(\Gamma\) is an \( N \times M \) matrix and \(Y\) is an \(N \times 1\) vector. The Moore-Penrose puesdo inverse can be employed to solve for \(\Theta\) analytically. The estimate of the unknown parameters is then $$ \hat{\Theta} = [\Gamma^T\Gamma]^{-1}\Gamma^T Y $$ For example if we fix the mass terms in the steer torque equation and let the rest be free the linear equation is $$ \begin{bmatrix} v(1) \dot{\phi}(1) & v(1) \dot{\delta}(1) & g \phi(1) & g \delta(1) & v(1)^2 \phi(1) & v(1)^2 \delta(1) & - F(1)\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ v(N) \dot{\phi}(N) & v(N) \dot{\delta}(N) & g \phi(N) & g \delta(N) & v(N)^2 \phi(N) & v(N)^2 \delta(N) & - F(N)\\ \end{bmatrix} \begin{bmatrix} C_{1\delta\phi}\\ C_{1\delta\delta}\\ K_{0\delta\phi}\\ K_{0\delta\delta}\\ K_{2\delta\phi}\\ K_{2\delta\delta}\\ H_{\delta F} \end{bmatrix} = \begin{bmatrix} T_\delta(1) - M_{\delta\phi} \ddot{\phi}(1) - M_{\delta\delta} \ddot{\delta}(1)\\ \vdots\\ T_\delta(N) - M_{\delta\phi} \ddot{\phi}(N) - M_{\delta\delta} \ddot{\delta}(N) \end{bmatrix} $$ The error in the fit is $$ \hat{Y} - Y = \Gamma \hat{\Theta} - Y $$ I'm not quite sure what the "noise" model is in this model structure. It doesn't seem to be strictly an output error model, as we are not necessarily finding only the error in the forcing predictions, as \(\hat{Y}\) can be a combination of kinematic and force measurements.

This equation can be solved for each run individually, a portion of a run or if there are \(P\) runs then they can be all collected in \(\Gamma\) for best fit over multiple runs. If all the runs have the same number of time steps, \(\Gamma\) then becomes an \( NP \times M \) matrix and \(Y\) an \(NP\times 1\) vector if each run is the same length.

Secondly, all of the parameters in the canonical matrices need not be estimated. The analytical benchmark bicycle model [Meijaard2007] gives us a good idea of which entries in the matrices we may be more certain about from our physical parameters measurements. I went through the benchmark formulation and eliminated free parameters based on these rules:

- If the parameter is greatly affected by trail, leave it free.
- If the parameter is greatly affected by the front assembly moments and products of inertia, leave it free.
- If the parameter is equal or near to zero, fix it.

For the roll equation this leaves \(M_{\phi\delta}\), \(C_{1\phi\delta}\), and \(K_{0\phi\delta}\) as free parameters. And for the steer equations this leaves \(M_{\delta\phi}\), \(M_{\delta\delta}\), \(C_{1\delta\phi}\), \(C_{1\delta\delta}\), \(K_{0\delta\phi}\), \(K_{0\delta\delta}\), \(K_{2\delta\delta}\), and \(H_{\delta F}\) as free. Previous attempts at identifying all of the parameters and identification from a state space perspective confirmed that this choice of fixed parameters was reasonable for the roll equation and that little is known in the steer equation.

I start by identifying the roll equation. I have much more certainty in the roll equation estimates from first principles. I then use the assumption that \(M_{\phi\delta} = M_{\delta\phi}\) and \( K_{0\phi\delta} = K_{0\delta\phi} \) to fix these values in the steer equation to the ones identified in the roll equation, leaving less free parameters in the steer equation. This matrix symmetry is likely enforced in reality due to the simple coupling of the front and rear frames. It is also possible to solve the roll and steer equations simultaneously and enforce the symmetry but this was the easier solution for the time being. Finally, I identify the remaining steer equation coefficients.

The code for this analysis can be found in https://github.com/moorepants/BicycleSystemID/tree/master/scripts.

## Results

I have data for three riders on the same bicycle, performing two maneuvers, on two different environments. I don't suspect that the dynamics of the system should vary much with respect to different maneuvers, but there is potentially a small variation across riders and there may be variation across environments because of the differences in the wheel to floor interaction. If computing the best fit model across a series of runs this leaves these four combinations:

- All riders in both environments (n=1)
- All riders in each environment (n=2)
- Each rider in both environments (n=3)
- Each rider in each environment (n=6)

### All riders in both environments

Here I'll make the assumption that the model doesn't vary much across riders or environments and give the results for the last item in the list. I run the fit over 374 runs giving about 142 minutes of data sampled at 200 Hz.

Matrix | First Principles (Whipple) |
Identified | Percent Difference |
---|---|---|---|

\(M\) | \(\begin{bmatrix} 129.86277082 & 2.29077079\\ 2.29077079 & 0.22039034 \end{bmatrix}\) | \(\begin{bmatrix} x & 2.28375768 \pm 7.05777612e-05\\ 2.28375768 & 0.20003257 \pm 1.94915567e-03 \end{bmatrix}\) | \(\begin{bmatrix} 0 & 0.30614621\\ 0.30614621 & 9.23713925 \end{bmatrix}\) |

\(C_1\) | \(\begin{bmatrix} 0. & 42.07257347\\ -0.3150094 & 1.3924634 \end{bmatrix}\) | \(\begin{bmatrix} x & 23.94008131\pm 3.15769856e-04\\ -0.88845094 \pm 4.33887933e-06 & 1.73368327 \pm 5.66893205e-07 \end{bmatrix}\) | \(\begin{bmatrix} 0 & 43.0981294\\ -182.03950055 & 24.50476423 \end{bmatrix}\) |

\(K_0\) | \(\begin{bmatrix} -114.59019659 & -2.36844326\\ -2.36844326 & -0.73938563 \end{bmatrix}\) | \(\begin{bmatrix} x & -3.91797216 \pm 1.94915567e-03\\ -3.91797216 & -0.66780731 \pm 2.30169695e-06 \end{bmatrix}\) | \(\begin{bmatrix} 0 & -65.42393996\\ -65.42393996 & -9.68078349 \end{bmatrix}\) |

\(K_2\) | \(\begin{bmatrix} 0. & 102.9447477\\ 0. & 2.19688326 \end{bmatrix}\) | \(\begin{bmatrix} x & x\\ x & 2.33973324 \pm 1.19898796e-06 \end{bmatrix}\) | \(\begin{bmatrix} 0 & 0\\ 0 & 6.50239301]] \end{bmatrix}\) |

\(H\) | \(\begin{bmatrix} 0.91545833\\ 0.01101888 \end{bmatrix}\) | \(\begin{bmatrix} x\\ 0.0086155 \pm 1.62592752e-09 \end{bmatrix}\) | \(\begin{bmatrix} 0\\ 21.81146898 \end{bmatrix}\) |

Notice that the identified model seems to have a similar eigenvalue locus as the Whipple model except for the higher weave speed and the faster caster mode. Notice the arm model has an unstable capsize mode for all speeds shown. There is a stable oscillatory mode for all speeds shown, with a similar frequency as the Whipple model for speeds above 2 m/s. In the Bode plots, the Whipple model seems to better predict the low speed behavior and the arm model better predicts the high speed behavior.

The question remains as to how well the identified model can predict the measured motion given measured inputs. This is ultimately our measurement of model quality. The following graphs give some example comparisons for various riders and various speeds of the simulation results of the model identified from all runs (all riders and both environments) along with the Whipple and Arm model.

These graphs give evidence that the identified model from all rider and both environments is better at predicting the measure motion of the bicycle given the measured inputs. The variance in the prediction is on the order of 50-80% for the identified model. I will have to run these fits on all of the runs to get an idea of the quality of the model, but this look very promising.

### All riders in each environment

The following results make the assumption that the differences in the riders' inertia are negligible, but that the treadmill and pavilion floor may give different results.

*All of the following graphs follow the same color and line style scheme as the ones above. I also haven't match frequencies in the phase plots.*

#### Horse Treadmill

#### Pavilion Floor

The environment doesn't seem to change the predicted dynamics too much.

### Each Rider on Both Floors

The following results assume that the environment has little difference on the models, but that the rider's do.

#### Charlie

#### Jason

#### Luke

Charlie's eigenvalue plot looks very different from the other two riders, with an unstable capsize mode similar to the arm model. It is interesting that the Arm model is different for Charlie than Jason and Luke. Luke and Jason's models seem similar to each other. The frequency response of Charlie doesn't seem to vary as much with respect to the other riders.

#### Each rider in each environment

I won't bore you with the 18 graphs, but they are all in this folder.

## Discussion

This method of estimating the bicycle models is the most promising that I've done so far. I am able to generate a single models for large sets of data. The resulting models act and look like bicycles, but I have enforced a lot of structure on the model by fixing a great deal of the parameters. The 4th order model still seems to be able to capture the essential dynamics and a modified steer torque equation which includes some kind of tire related torques may be able to account for the differences in the identification and the first principle models. These results give less credence to the Arm model and in the process I realized that I haven't properly linearized the Arm model in most of my previously report work. When I try to reduce my correct 19 x 19 A matrix to it's essential 4 x 4 size, I'm not properly accounting for the holonomic constraints. I hope to have that worked out so that I can compare the numerical values of the matrix entries.

The model predictions from the identified model simulation look very promising. My main question at this point, is how this model addresses the process and measurement noise and whether the parameter identification is very accurate.