CALSPAN¶

The earliest comprehensive bicycle model validation began at CALSPAN in the late 60’s. This included several revolutionary studies, in one of which they made use of a rocket to apply know step torques to an uncontrolled riderless bicycle. In another, simulations of slalom maneuvers were visually compared with video footage [RM71].

Eaton¶

David Eaton’s work ([Eat73b], [Eat73a], [ES73]) may be the closest example to the work presented in this chapter. He did his PhD work at the University of Michigan under the Highway Safety Research Institute. His dissertation focused on the experimental validation of the motorcycle modeling work of [Sha71] and the human controller modeling work of [Wei72]. He did this with two sets of experiments: 1) identification of the uncontrolled dynamics of the motorcycle under perturbations, and 2) identification of the rider controller during roll stabilization tasks. The latter of which will be discussed in the next section.

His initial experiments were aimed at validating and identifying the passive motorcycle system. During these experiments, his subjects road a motorcycle with their bodies rigidly braced to the frame and hands-free at speeds of 15, 30, and 45 mph (6.7, 13.4, and 20.1 m/s) along side a pace car which recorded the output from roll angle, roll rate, and steer angle sensors. The brace and open loop response allowed rigid rider modeling assumptions to be used. Weights were dropped from one side of the motorcycle to induce a step roll torque and the rider used a single pulse in steering torque to the handlebars to right the motorcycle in roll after the drop. These experiments were impressively dangerous and would be hard pressed for approval by the Institutional Review Board if done today, but well designed for the typical modeling assumptions. The resulting time histories of the measured system outputs were compared to simulations of Sharp’s model [Sha71] augmented with a variety of tire models of Eaton’s design. He found good agreement between the experiments and the models for higher speeds, but felt that a more complex tire model was needed to predict the wobble mode in slower speed runs.

The second set of experiments were more tame. The three riders simply balanced the motorcycle on a straight path at two speeds, 15 mph and 30 mph, for a total of 38 runs. He added a steer torque transducer bar above the handlebars. The rider controlled the motorcycle with one hand and the rider applied torque was recorded along with the other signals. No perturbations were necessary, as the rider’s natural control actions excited the system in a wide enough bandwidth. From this data he was able to identify the motorcycle steer torque to roll angle transfer function through the spectral densities of the measured signals (by dividing the cross spectrum of the roll angle and steer torque signal by the power spectrum of the steer torque). The identified transfer functions show good agreement with the augmented Sharp motorcycle model at the 30 mph speeds, but less so for the 15 mph runs.

His generated frequency responses from the second experiments provided an empirical model, while the simulation comparisons from the first experiments were validation rather than identification.

Weir, Zellner, Teper¶

Weir, Zellner, and Teper performed an extensive experimental study on motorcycle handling qualities for the U.S. National Highway Traffic Safety Administration in the late 70’s, [WZT79]. This was a follow up to both the CALSPAN studies and [Tag75], both under or related to the same Administration. There is little to no explicit system identification in the study but some important elements are there. In terms of the passive model identification, they present steady state comparisons of their experimental data to their models with varying degrees of qualitative agreement and generally good ability to predict the conditions at which sign reversals in torque are needed to maintain a steady turn. They also compare single lane change simulations of a controlled vehicle to their measured data by visual inspection. They unfortunately admit that adjusting the first principles models to better fit their measured data was outside the scope of the project. But this gives some early examples of model evaluation with respect to good quality data.

James¶

Stephen James published a study in 2002 [Jam02] in which he attempted to identify the linear dynamics of an off-road motorcycle. He measured steering torque, steer angle, speed, roll rate, and yaw rate with his subjects manually exciting the vehicle through steer torque during runs at various speeds on a straight single lane road. He made use of black box ARX SIMO identification routines of 6th and 7th order (his and others motorcycles models are usually 10th+ order) to tease out the weave and wobble eigenvalues. He compares the identified eigenvalues, eigenvectors and frequency responses to his motorcycle model and claims good fits based on visual interpretation of the plots. The agreement is questionable due to the lack of statistics in the model comparisons and little validation of his first principles model which assumes a rigid rider. The study does show that there is the possibility of identification of multiple modes of motion with simple manual excitation of the handlebars. He also used these techniques to identify the same motorbike with a single wheel trailer in [Jam05].

Biral et al.¶

[BBCL03] performed a nice study to identify motorcycle dynamics under an oscillatory steer torque input. They measured steer torque, roll rate, steer angle, and yaw rate with an instrumented motorcycle. They performed slalom maneuvers at speeds from 2 to 30 m/s at three sets of cone spacings in the slalom course. The resulting time histories were close to ideal sinusoids. They used curve fitting to find amplitude and phase relationships among the measured signals. The results were plotted on Bode plots for comparison to the frequency response of several first principles models. The models predict the experimental data and their motorcycle model is shown to do a better job than other models from literature. This claim is only based on visual inspection. I would say this technique and others like it are more of an ad hoc method of system identification of the vehicle dynamics because they rely heavily on very specific input and output characteristics, but nevertheless seems to be effective. Making use of formal system identification techniques could potentially give more reliable results and the ability to better characterize the uncertainty in the predictions.

Kooijman¶

Jodi Kooijman has worked on experimental validation of the benchmark bicycle [MPRS07] linear equations of motion for a riderless bicycle [Koo06], [KSM08], [KS09]. His instrumented bicycle measured the steer angle, forward speed, roll rate, and yaw rate. Due to the fact that the bicycle can be stable at certain speeds he was able to launch the bicycle in and around the stable speed range and perturb the bicycle with a lateral unmeasured impulse and record the stable decay in the steer, roll, and yaw rates. The post perturbation time histories of the measured signals provided nice decaying oscillations and curves could be fit to find both the time decay constant and frequency of oscillation. These were then compared to the predicted weave response based on the first principle model numerically populated with measured physical parameters of the bicycle. He found good prediction abilities of the weave mode between 4 and 6 m/s. The “goodness” of fit were gaged by visual inspection with no uncertainty estimates in the models or the results from the dynamic measurements. The method was not able to predict the heavily damped caster mode nor the capsize mode. He also demonstrated that the measured dynamics were the same when the experiments were performed on a treadmill.

In [KMP+11], Jodi constructed a bicycle with very unusual physical characteristics including negative trail and canceled angular momentum of the wheels. He performed similar experiments to his Master’s thesis work. They show the comparison of a single stable experiment in which yaw and roll rates were measured and compared to the predictions of the benchmark bicycle.

[Ste09] and [ER11] both perform experiments similar to Kooijman’s with similar results, although Steven’s results vary in the ability of the model to predict the data for various configurations of his adjustable bicycle.

These also fall into the ad hoc system identification techniques that take advantage of the stability at certain speeds and very specific output characteristics. The variability in reproducibility in the studies from other researchers should be noted.

Chen and Doa¶

[CD10] develop a first principles non-linear bicycle model with a fuzzy controller and use it to generate stable simulations for various speeds. They then do an output error grey box identification on the resulting data with respect to the non-zero and non-unity entries of the state, input and output matrices (i.e. just the entries of the acceleration equations). The identification is done for a discrete number of speeds in the range 1 to 15 m/s. The eigenvalues are calculated of the resulting identified, speed-dependent A matrices and the root locus plotted versus speed.

The resulting eigenvalues seem to behave like the benchmark bicycle but the capsize mode is shown to go unstable briefly at a speed lower than the stable speed range. They did not attempt to characterize or identify the process noise even though they generated the data with a known model with known input noise. Also their non-linear bicycle equations of motion [CD06] were never validated against any other accepted models. Both of these could potentially explain the discrepancies in their identification. Their identification procedure does show that it may be possible to get good estimates of a linear model of the vehicle alone from noisy data regardless of the controller which stabilizes the vehicle.

Doria¶

In [DFT12] experiments are performed where a motorcycle rider excites the steering with a pulse and lets the motorcycle oscillate while the rider keeps his hands on the handlebars (as opposed to Eaton’s hands-free experiments). The resulting dynamical measurements are nice decaying sinusoidal-like motions and the authors fit optimal curves to the data. They identify the time constants, and frequency and phase information to construct the eigenvalues and eigenvectors of the excited mode. The empirically derived eigenvectors show some resemblance to the model’s predictions.

van Lunteren and Stassen¶

At Delft University of Technology in the Man-Machines research group, Drs. van Lunteren and Stassen began work in 1962 to identify the human controller for a normal population of subjects and report on their work into the early 70’s ([VLS67], [VLS69], [Sta69], [VLS70a], [VLS70b], [VLS70c], [VLS73], [SVLB+73]). They chose a bicycle simulator as the plant because it was a common task that average people could do and their studies could focus on a wider population of individuals as compared to most previous work based around trained pilots. The bicycle simulator did not capture all of the essential dynamics of a real bicycle as it’s operation was based on only the simplified roll dynamics of Whipple’s model, but nonetheless offered a similarly complex roll stabilization control task as a normal bicycle would. The simulator was controlled by both the steering angle and the rider’s lean angle, both of which are questionable inputs as have been pointed out as early as [RL72].

They assumed the rider’s control actions can be described by a PID controller with time delays on each feedback variable and mention that this controller was chosen instead of a McRuer style controller primarily due to limitations of their computational equipment. The error in the roll angle is fed into two PID controllers each with a time delay: one to output the corrective steer angle and the other to output the corrective lean angle. They introduce a remnant term for each control action and the external disturbances to the bicycle model.

The identification goal was to find the six gains and two time delays in which the controller performed as a human would. The preferred method was a real time estimation routine due to the speed of computations and reasonable agreement their correlation method. The results indicated that the subjects used no integral control (i.e. only position and rate feedback). They could identify within a bandwidth of about 2 Hz and noticed that when the system was undisturbed there was a 0.5 Hz dominant frequency in the rider’s control actions. The rate feedback was more dominant in generating the lean control input than it was for the steer control input. Also, they found the time delay for lean to be larger than the steer time delay and postulate that the steer action is a result of cerebral activity while the lean is more of a reflexive pattern. Another finding resulting from analysis of Nyquist plots of different riders’ identified control actions showed that riders chose different control actions. They attribute this to the roll stabilization being a sub-critical task (i.e. a more difficult task may force different riders to adopt similar control behavior). They also investigated the effects of drugs, such as alcohol, on the riders control behavior. They found correlations from drug dose to time delays and the error in the control actions. Their later studies introduced better identification methods and they found discrepancies in the identified time delays of the later work as compared to the previous work. For example, the steer control time delay was originally found to be around 1.5 seconds and the improved methods found the delay to be around 0.7 seconds. The discrepancy was attributed to the bias due to remnant in their early work. They also introduced a visual tracking task into the simulator but had difficulties in getting reliable transfer function identification as compared to the roll stabilization transfer functions which improved in quality due to longer trials of 35 minutes.

The methods developed in their studies are excellent and thorough examples of early parameter identification in human control tasks. The simpler plant dynamics were most likely beneficial in reducing the uncertainty in the identified parameters, but the choice of angles as inputs instead forces of torques may not be a realistic enough model of the human’s actuation.

Eaton¶

After feeling confident in his motorcycle identification results, Eaton made use of the Wingrove-Edwards ([WE68], [WEC69], [Win71], [Edw72]) method in tandem with an impulse identification to identify the human controller. The remnant element was large with respect to the torque that was linearly correlated with the roll angle, but the human control element was identified with a simple gain and time delay for most of the high speed runs. The time delay identification of about 0.3 seconds was very repeatable across all runs. Furthermore, he demonstrated that the crossover model was evident in the resulting closed loop rider-motorcycle transfer functions.

Eaton is one of very few who have identified the rider controller during actual single track vehicle tests with confidence in the underlying passive rider-vehicle model. This study has influenced the work in this chapter in many ways.

Doyle¶

A recently uncovered study on the manual control of a bicycle from a psychologist’s perspective has some very non-traditional techniques and outlooks for the understanding of the control system employed while balancing a bicycle [Doy87]. Anthony Doyle’s paper [Doy88] on his thesis topic opens with “The old saw says that once learned it is never forgotten, but what exactly is learned has been by no means clear.” This reflection points to the great complexity behind balancing a bicycle, such an easily gained skill. He chooses to study the bicycle over a simpler task partially due to the fact that the rider has little freedom in effective control strategies and partially because it is a skill many people can do.

His goal was to determine how much of the rider’s control actions can be accounted for without involving higher cerebral functions. He mentions the Weir and Zellner work and the fact that its focus is on motorcycles at high speed, and questions whether the control employed for their system is simply a different version of the one employed on a bicycle at low speed or whether they are different control methodologies altogether.

He was aware of the inherent stability that bicycles can provide and constructs an instrumented bicycle where the head angle, trail, and front wheel gyro effects are eliminated so that “all steer movements are a result of the human’s control”. He also mentions, but doesn’t use, a body brace to eliminate unnecessary body movements and he blindfolds his subjects so that their sensory information is limited to proprioception and vestibular cues. He mentions the arm and upper body movements and how it is difficult to tease out the deliberate movements versus the passive dynamics of the body. With the instrumented bicycle he conducts low speed steady turn and balancing tasks and measures speed, roll rate, and steer angle.

Along with the experimental data, he developed a bicycle and rider model with accompanying controller. The derivation of the bicycle model is questionable due to the non-traditional methods, but he does end up with a model which behaves like a bicycle including speed dependent stability. He is aware of the need to roll the bicycle frame in the direction of the desired turn for directional control and how counter-steering plays a roll in this. This concept leads to the primary inner loop being chosen as roll control and his control structure resembles that of Weir’s work in terms of sequential loops. He cites the crossover model and is aware that humans can adjust their gains as needed for good performance. The controller is traditional in most senses and follows the patterns by McRuer, Weir, and Eaton, but he adds in the ability to add discrete pulses to the roll angle. He feeds back roll acceleration and integrates it to get roll angular velocity. This is basically a continuous PD control on roll rate. But his non-continuous addition to the controller is based on a fuzzy logic-like rule “Make a pulse against the lean whenever it gets bigger that 1.6 degrees.”

It seems like he gets somewhat close matches to the experimental traces from his control model simulations without the discrete pulses, but then adds in pulses (single or multiple) to the steering so that the traces matches more closely. His identification technique and criterion is focused around a detailed examination of the patterns in the time histories in a very qualitative way.

His results focus on the evidence for intermittent control and finds the traditional gains to be inversely proportional to speed. He claims the balancing part of the control system is done primarily in the lower cortex.

To me, Doyle’s work emphasizes the need for close collaboration between psychologists and control engineers to formalize the theory for human balance. His intermittent control theory may be valid, but due to the unusual model development, simulation and analysis techniques it is hard to gauge whether the need for intermittent control was simply an artifact of poor vehicle modeling. His insight into the human control theory is very enlightening and his ways of wording bring the theory outside the traditional control framework for an expansion in understanding.

Lange¶

Peter de Lange’s recent Master thesis work [DL11] focused on identifying the rider controller from the data that he helped us collect while interning at our lab. He used the Whipple bicycle model, a simplified second order representation of the human’s neuromuscular dynamics (natural frequency 2.17 rad/s and damping ratio of 1.414[2]) and a PID like controller with a 0.03 second time delay. The controller structure had gains proportional to the integral of the angle, the angle, the angular rate and the angular acceleration for roll and steer. The control task was defined as simple roll stabilization (i.e. track a roll angle of zero degrees), even though the data was collected during heading and roll tracking tasks.

He used a four step process for identifying the rider controller 1) he “removed” the human remnant by averaging the time histories over several individual perturbations, 2) he identified a high order finite impulse response model (only a function of previous inputs) for the lateral force to steer angle SISO pair (lateral perturbation force input and steer angle as output) 3) low pass filtered the resulting responses, and 4) he identified the rider controller parameters with a grey box model using the filtered FIR simulation results as the base data. The grey box model was parameterized with eight gains and a time delay. He was able to identify the gains, but the time delay always gave a resulting unstable model, so he dropped it. Furthermore, all of the gains were not necessary for good model predictions so he eliminated the unnecessary gains systematically to find the critical feedback elements. These turned out to be the gains for roll angle, roll rate, steer rate, and the integral of the steer angle. The first three are as one may expect and he concludes that the steer angle integral could be equated to yaw angle feedback since they are proportional in the linear sense.

Peter’s approach hinges on the averaging process in step one. The human remnant is large relative to the measurements and averaging potentially removes data that isn’t necessarily noise. This averaging is atypical, as process noise models are usually employed to account for these variations in the data. Using a model such as ARMAX instead of the two step averaging and FIR model would potentially allow one to identify the underlying linear model without removing potentially valid data in the time history averaging process. Or all of the steps could be combined into a state space grey box formulation with a process noise model, for a more direct route to identifying the free parameters. But these methods have their difficulties and will be described later in the Chapter.

Treadmill¶

Dr. James Jones at the veterinary school at here at Davis graciously let us use their horse treadmill (Graber Ag Kagra Mustang 2200) during their downtime, Figure 11.1. The treadmill is 1 meter wider and 5 meters long and has a speed range from 0.5 m/s to 17 m/s. This was only a third of the width treadmill at Vrije Universiteit in Amsterdam, but after some practice runs we felt that narrow lane changes and the lateral perturbations could be successfully performed. We used the treadmill because the environment was very controllable, in particular with regard to fixed constant speeds, and it offered the ability to do very long run durations within a broad speed range. Potentially both the side railings and the belt side curbs added to rider’s lack of lateral movement space and changed the riders’ control strategy.

Sideview of the horse treadmill while Luke was riding the bicycle.

Pavilion¶

The bicycle was designed in such a way that all of the data collection equipment was on board and was suitable for data collection in a free environment. After lengthy bureaucratic negotiations, we were able to make use of the UCD pavilion floor for the experiments, Figure 11.2. The floor was made of a stiff rubber[3] and provided a rectangular, wind-free space of about 100’ by 180’ (30 m by 55 m). We rode around the perimeter to build up speed and did our maneuvers on a straight section about 100 feet (30 m) long. We were not able to travel at speeds higher than about 7 m/s as the tires would slip in the final turn into the test section (this seemed to be due to the dust on the floor). This indoor environment provided a wind free area which was more akin to the environment bicyclists normally ride in than did the treadmill.

Overhead view of the pavilion floor during a perturbation run.

Method¶

For this analysis, I limit the inputs to steer torque and lateral force and the outputs to roll angle, steer angle, roll rate, and steer rate. The ideal fourth order system can be described with the following continuous state space description

where \(\eta(t)\) are the outputs and \(\mathbf{H}\) is the identity matrix.

Assuming that this model structure can adequately capture the dynamics of interest of the bicycle-rider system, our goal is to accurately identify the unknown parameters \(\theta\) which are made up of the unspecified entries in the \(\mathbf{F}\) and \(\mathbf{G}\) matrices. To do this one needs to recognize that this continuous formulation is not compatible with noisy discrete data. The following difference equation can be assumed if we sample the continuous system at \(t=kT\), \(k=1,2,\dots\), with \(T\) being the sample period and making the assumption that the variables are constant over the sample period (i.e. zero order hold).

The additional terms \(w\) and \(v\) represent the process and measurement noise vectors, respectively, which are assumed to be sequences of white Gaussian noise with zero mean and some covariance. By making use of the Kalman filter, Equation (2) can be transformed such that the optimal estimate of the states, \(\hat{x}\), with respect to the process and measurement noise covariance are utilized, see [Lju98].

where \(\mathbf{K}\) is the Kalman gain matrix. \(\mathbf{K}\) is a function of \(\mathbf{A}(\theta)\), \(\mathbf{C}(\theta)\) and the covariance and cross covariance of the process and measurement noises, but it can also be directly parameterized by \(\theta\). With that, this equation is called the directly parameterized innovations form [Lju98] and the entries of the four matrices in equation (3) can be estimated directly.

The \(\mathbf{A}\) and \(\mathbf{B}\) matrices are related to \(\mathbf{F}\) and \(\mathbf{G}\) by

where \(T\) is the sampling period. With a linear assumption \(\mathbf{A}\) and \(\mathbf{B}\) can even be directly estimated in discrete form by

The one step ahead predictor for the innovations form is

where \(q\) is the forward shift operator (\(q u(t) = u(t+1)\)) [Lju98]. The predictor is a vector of length \(p\) where each entry is a ratio of polynomials in \(q\). These are transfer functions in \(q\) from the previous inputs and outputs to the current output. In general, the coefficients of \(q\) are non-linear functions of the parameters \(\theta\).

We can now construct the cost function, which will enable the computation of the parameters which give the best fit using optimization methods. We’d like to minimize the error in the predicted output with respect to the measured output at each time step. First form \(Y_N\) which is a \(pN \times 1\) vector containing all of the current outputs at time \(kT\).

where \(p\) is the number of outputs and \(N\) is the number of samples. Then organize the predictor vector, \(\hat{Y}_N(\theta)\), the one step ahead prediction of \(Y_N\) given \(y(s)\) and \(u(s)\) where \(s \leq t - 1\)

The cost function is then the norm of the difference of \(Y_N\) and \(\hat{Y}_N(\theta)\) for all \(k\).

The value of \(\theta\) which minimizes the cost function is the best prediction

where \(Z^N\) is the set of all the measured inputs and outputs.

In general, the minimization problem is not trivial and may be susceptible to many of the issues associated with optimization including local minima. The number of unknown parameters in the \(\mathbf{K}\) matrix is a function of the number of states and the number of outputs, in our case in \(\mathbf{R}^{4\times4}\) which more than doubles the number of unknowns present in the \(\mathbf{A}\) and \(\mathbf{B}\) matrices. It is thus critical to reduce the number of unknown parameters to have a higher chance of finding the global minimum of the cost function. The accuracy of the system parameters depends on the ability to estimate the \(\mathbf{K}\) matrix along with the other parameters.

Before identification I further processed all of the signals that were generally symmetric about zero by subtracting the means over time. For some of the pavilion runs, this may have actually introduced a small bias, as the short duration runs with unbalanced perturbations may not have zero mean.

I made use of the Matlab System Identification Toolbox for the identification of the parameters \(\theta\) in each run of this model structure. In particular, a structured idss object was built with the initial guesses of the unknown parameters based on the Whipple model and the initial guesses for the initial conditions and the Kalman gain matrix being equal to zero. All of my attempts at identifying the Kalman gain matrix were plagued by local minima.

Results¶

It turns out that finding a model that meets the criterion is not too difficult when the output error form is considered (\(\mathbf{K}=0\)). This model may be able to explain the data well, but the parameter estimation is potentially poor because the parameters in the state and input matrices are adjusted such that the results fit both the true trajectories and the noise. Global minima in the search routine are quickly found when the number of parameters is between 10 and 14. When the \(\mathbf{K}\) matrix is added the number of unknown parameters increases by 16 and the global minimum becomes more difficult to find and I was rarely able, if ever, to find the global minimum for the general problem, even when reducing the number of outputs to one.

Figure 11.7 shows typical example input and output data for a single run (#596) with both steer torque and lateral force as inputs. The plot compares the simulation response of the input to the measured response. Notice that the identified model predicts the trajectory extremely well. Similar results are found for the majority of the runs. The Whipple model predicts the trajectory directions but the magnitudes are large, meaning that for a given trajectory, the Whipple model requires less torque than that measured. The Whipple model with the arm inertial effects does a better job than the basic Whipple model, but still has some magnitude differences with respect to the identified model. It also has a harder time predicting the roll angle than both the identified model and the Whipple model.

Example results for the identification of a single run (#596). The experimentally measured steer torque and lateral force are shown in the top two graphs. All of the signals were filtered with a 2nd order 15 hertz low pass Butterworth filter. The remaining four graphs show the simulation results for the Whipple model (W), Whipple model with the arm inertia (A), and the identified model for that run (I) plotted with the measured data (M). The percentages give the percent of variance explained by each model.

The identified models are almost always unstable due to the high weave critical speed and, even though the measured inputs stabilize the true system, they will not necessarily stabilize the models. This poses an issue when gauging the model quality by the percentage variance of the output data explained by the model. A model that blows up during the simulation may not necessarily be a bad model, but it will return a very small percent variance and lose its ability to be compared by that criterion. [BBCL03] and [TJ10] both are able to run feed-forward simulations of their motorcycle models with the measured steering torque. They both are dealing with high speed motorcycles which typically only have a slightly unstable capsize mode. [TJ10] uses a controller to compensate the torque for unbounded errors so that the simulation doesn’t blow up. The method I use here is to choose short duration portions of the runs for simulation and search for the best set of initial conditions to keep the model stable during the duration. This generally works but there are ultimately some incomparable runs due to this issue.

I use this structured state space output error identification procedure for a collection of experiments (\(n=368\)) over a range of speeds between about 1 and 9 m/s. Figures 11.8 and 11.9 plot the identified coefficients of the dynamical equations of motion (i.e. the bottom two rows of the \(\mathbf{F}\) and \(\mathbf{G}\) matrices) as a function of speed for all of the experiments using box plots. Both the Whipple (green) and arm (red) model predictions are superimposed for comparison. The first notable thing is that the coefficients seem to generally have large variance, especially as the speed increases. Secondly, the roll acceleration, \(\ddot{\phi}\), equation seems to be better predicted by the two models and the data has less spread at the lower speeds, barring the \(\dot{\phi}\) coefficient which has large spread and no apparent relationship with speed for both equations. The roll equation also seems to have less spread in the experimental data. For example, the \(a_{\ddot{\phi}\delta}\) coefficient appears to be very tight and the first principles models predict it very well. The constant, linear, and quadratic trends in the coefficients are somewhat visible in the data but the variance in the coefficients clouds it. This variability in the coefficient predictions depends on many things including data quality, the ability to identify a process noise model, speed being constant during the run, choice of unknown coefficients, and more. With all of these improved, detailed regression models may be able to reveal the true trends[4]. Nonetheless, these graphs reveal several important things:

• The identified models predict the data well with most having mean predicted variance of the four outputs above 70% (but this is tightly coupled to run duration, i.e. longer runs have more room for model error).
• Some of the coefficients are well predicted by the Whipple model and can be fixed from first principles calculations, notably: \(a_{\ddot{\phi}\phi}\), \(a_{\ddot{\phi}\delta}\) and \(b_{\ddot{\delta}T_\delta}\) and maybe even \(a_{\ddot{\delta}\delta}\).
• The roll rate coefficients are highly variable with poor prediction by the models. Deficiencies in the first principles models are likely.
• Either the higher speed runs are outliers, or the behavior of the system changes more rapidly with speeds above 5 m/s or so.
• Some coefficients spread around zero giving inconsistent signs and others give opposite signs from those the first principles models expect.

State coefficients of the linear dynamical equations of motion plotted as a function of speed. Each box plot represents the distribution of that parameter for a small range of speeds, i.e. speed bin. The width of the box is proportional to the total duration of the runs in that speed bin. The green line is the Whipple model and the red line is the arm model. Only experiments with a mean fit percentage greater than zero are shown. The orange line is the model identified with the canonical method using runs done by Luke in the pavilion which is presented and discussed in the next section. Generated by src/systemidentification/coefficient_box_plot.py.

Input coefficients of the linear dynamical equations of motion plotted as a function of speed. Each box plot represents the distribution of that parameter for a small range of speeds, i.e. speed bin. The width of the box is proportional to the total duration of the runs in that speed bin. The green line is the Whipple model and the red line is the arm model. Only experiments with a mean fit percentage greater than zero are shown. The orange line is the model identified with the canonical method using runs done by Luke in the pavilion which is presented and discussed in the next section. Generated by src/systemidentification/coefficient_box_plot.py.

Figure 11.10 gives another view of the resulting data. It is a frequency response plot at the mean speed for a set of runs. The blue lines give the mean and one standard deviation bounds of the magnitude and phase of the system transfer function \(\frac{\phi}{T_\delta}(s)\) for the set of runs. Even though the spread in the identified parameters seems high in Figures 11.8 and 11.9, the Bode plot shows that the identified system frequency response is not as variable, especially in magnitude. It is also apparent that the experimental magnitude mean has a -5 to -10 dB offset across the frequency range shown with respect to the Whipple model, although the Whipple model does fall within one standard deviation of the mean. This correlates with the amplitude differences in the trajectories shown in Figure 11.7. Notice that the arm model has little to no offset between 2 and 10 rad/s, thus the better amplitude matching. The frequency response may give a better indication of the overall identified model quality.

Frequency response of steer torque to roll angle for a set of runs at \(4.0 \pm 0.3\) m/s. The solid blue line is the mean from the identified runs and is bounded by the standard deviation, the dotted blue lines. The green line is the Whipple model and the red line is for the model which accounted for the arm inertial effects.

Conclusion¶

I have shown that a fourth order structured state space model is both adequate for describing the motion of the bicycle under manual control in a speed range from approximately 1.5 m/s to 9 m/s. The fact that higher order models may not be necessary for bicycle dynamic description is an important finding. More robust models of single track vehicles are typically higher than 4th order, with degrees of freedom associated with tire slip, frame flexibilities, and rider biomechanics. These findings suggest that the more complex models may be overkill for many modeling purposes. The data subsequently also reveals that fourth order archetypal first principles models are not robust enough to fully describe the dynamics. The deficiencies are most likely due to un-modeled effects, with the knife-edge, no side-slip wheel contact assumptions being the most probable candidate. Un-modeled rider biomechanics such as passive arm stiffness and damping and head motion may play a role too. The uncertainty in the estimates of the physical parameters from Chapter Physical Parameters is not large enough to explain the difference between the identified coefficient identification and their predictions from first principles. It is likely that something as simple as a “static” tire scrub torque is needed to improve the fidelity of the first principles derivations, but that doesn’t preclude that the additional of a tire slip model might also improve the models.

Model structure¶

The identification of the linear dynamics of the bicycle can be formulated with respect to the benchmark canonical form realized in [MPRS07], Equation (11). If the time varying quantities in the equations are all known at each time step, the coefficients of the linear equations can be estimated given enough time steps.

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\). This equation assumes that the velocity is constant with respect to time as the model was linearized about a constant velocity equilibrium, but the velocity can also potentially be treated as a time varying parameter if the acceleration is negligible. 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.

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 torque 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. For our normal geometry bicycles, including the one used in the experiments, \(H_{\delta F} << H_{\phi F}\). I estimate \(H\) for each rider/bicycle from geometrical measurements and the state space form of the linear equations of motion calculated in Chapter Extensions of the Whipple Model.

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 partitioned.

where \(\mathbf{A}_l\) and \(\mathbf{A}_r\) are the 2 x 2 sub-matrices corresponding to the states and their derivatives, respectively. \(\mathbf{B}_F\) and \(\mathbf{B}_T\) are the 2 x 1 and 2 x 2 sub-matrices corresponding to the lateral force and the torques, respectively. The benchmark canonical form can now be written as

where

The location of the lateral force application point is the same for Charlie and Luke because they used the same seat height. The force was applied just below the seat, which was adjustable in height for different riders.

Data processing¶

Chapter Davis Instrumented Bicycle details how each of the signals were measured and processed. For the following analysis, 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 with a central differencing method except for the end points being handled by forward and backward differencing. The mean was subtracted from all the signals except the lateral force.

Identification¶

A simple analytic identification problem can be formulated from the canonical form. If we have good measurements of \(q\), their first and second derivatives, forward speed \(v\), and the inputs \(T_\delta\) and \(F\), the entries in \(\mathbf{M}\), \(\mathbf{C}_1\), \(\mathbf{K}_0\), \(\mathbf{K}_2\), and \(H\) can be identified by forming two simple regressions, i.e. 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 as it has some variability.

The roll and steer equation each can be put into a simple linear form

where \(\Theta\) is a vector of the unknown coefficients and \(\mathbf{\Gamma}\) and \(Y\) are made up of the inputs and outputs measured during a run. \(\Theta\) can 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 \(\mathbf{\Gamma}\) is an \(N \times M\) matrix and \(Y\) is an \(N \times 1\) vector. The Moore-Penrose pseudo inverse can be employed to solve for \(\Theta\) analytically. The estimate of the unknown parameters is then

For example, if we fix the mass terms in the steer torque equation and let the rest be free the linear equation is

The error in the fit is

The covariance of \(\Theta\), Equation (23), of the parameter estimations can be computed with respect to the error.

Equations (19), (21), (22), and (23) can be solved for each run individually, a portion of a run, or a set of runs. Secondly, all of the parameters in the canonical matrices need not be estimated. The analytical benchmark bicycle model [MPRS07] gives a good idea of which entries in the matrices we may be more certain about from our physical parameters measurements in Chapter Physical Parameters. I went through the benchmark formulation and fixed the parameters based on these rules, keeping in mind that even little numbers can have large effects to the resulting dynamics:

• 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 identically equal to zero, fix it.

The reasoning for these assumptions are:

• Trail is difficult to measure and be certain about, especially since the bicycle tire deforms and creates a variable shaped tire contact patch which depends on the tires properties, pressure, and the configuration of the bicycle. The true trail is what is typically called pneumatic trail and gives the location in the tire patch at which the resultant contact force acts.
• The front frame moments and products of inertia play a large role in the steer dynamics and I’m not as confident in the estimation of these due to the fact that our apparatus was more suited to the estimate the inertial properties of the rear frame than the front.
• The zero entries in the velocity dependent stiffness and damping matrices that correspond to the roll angle and rate are assumed to hold from first principles.

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 equation 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 parameters.

I start by identifying the three coefficients of the roll equation for the given data. This choice is due to there being more certainty in the roll equation estimates from first principles because there are fewer unknown parameters.

I then enforce the assumptions 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 fewer free parameters in the steer equation[5]. This symmetry is enforced to coincide with the theory and choice of coordinates. Finally, I identify the remaining steer equation coefficients with

Results¶

I selected data for three riders on the same bicycle, performing two maneuvers, in two different environments. There is little reason to believe the dynamics of the passive system should vary much with respect to different maneuvers, but there is potentially variation across riders due to the differences in their inertial properties and there may be variation across environments because of the differences in the wheel to floor interaction. I opted to compute the best fit model across series of runs to benefit from the large dataset. This leaves these four scenarios:

• All riders in both environments, one data set
• All riders in each environment, two data sets
• Each rider in both environments, three data sets
• Each rider in each environment, six data sets

A total of 12 different models can be derived from this perspective.

All riders in both environments¶

This section details the example results for one subset of data. Here I make the assumption that the best fit model doesn’t vary much across riders or environments. The assumption that the passive model of the bicycle and rider are similar with respect to rider can be justified by recognizing that the Whipple model predicts little difference in the dynamics with respect to the three bicycle/rider combinations. On the other hand, I have little reason to believe the environments are the same except that both floors are made of a rubber like material. I calculated the best fit over 374 runs giving about 142 minutes of data sampled at 200 Hz, \(N=1720647\).

The eigenvalues as a function of speed of the identified model can be compared to those of the Whipple and arm models, see Section Rider Arms in Chapter Extensions of the Whipple Model. Figure 11.11 shows the root locus of the three models. The oscillatory weave mode exists in all three models, with it always being stable in the arm model but being unstable at lower speeds in the other two models. The identified model’s oscillatory weave mode is unstable over most of the shown speed range. Above 3 m/s or so, the Whipple model’s oscillatory weave mode diverges from the identified model to different asymptotes. The arm model weave mode diverges somewhere in between. Note that the arm model has an unstable real mode for all speeds. Figure 11.12 gives a different view of the root locus allowing one to more easily compare the real parts of the eigenvalues. The imaginary parts of the weave mode have similar curvature with respect to speed for all the models, with the identified model having about 1 rad/s larger frequency of oscillation for all speeds. The identified model does have a stable speed range where the Whipple model under predicts the weave critical speed by almost 2 m/s. The identified caster mode is much faster than the one predicted by the Whipple model which is somewhat counterintuitive because tire scrub torques would probably tend to slow the caster mode. Although, the pneumatic trail and the rider’s arm inertia could play a larger role that expected. Furthermore, the caster mode may not be well identified because the experiments were not design to specifically excite it.

Root locus of the identified model (circle), the Whipple model (diamond), and the arm model (triangle) with respect to speed in m/s. Generated by src/systemidentification/canonical_plots.py.

Real and imaginary parts of the eigenvalues as a function of speed for model (I)dentified from all runs, the (W)hipple model and the (A)rm model. Generated by src/systemidentification/canonical_plots.py.

The identification process is based on identifying the input/output relationships among measured variables. The frequency response provides a view into these relationships. Figures 11.13 to 11.16 give a picture of how the first principles models compare to the identified model with respect to frequency response from a roll torque input. The frequency band from 1 rad/s to 12 rad/s is of most concern as it bounds a reasonable range that the human may be able to operate within.

The roll torque to roll angle response Figure 11.13 shows that at 2 m/s (solid lines) the response of the Whipple model and the identified model are practically identical across the frequency range shown. On the other hand, the arm model predicts the high frequency (> 4 rad/s) behavior really well, but the magnitude is as great as 50 dB larger and the phase off by 180 degrees at lower frequencies. At 4 m/s the Whipple model still predicts the phase well, but the magnitudes do not match below 100 rad/s with 10 dB difference at lower frequencies. At this speed the arm model matches the identified model as low as 3.5 rad/s but the low frequency phase still being 180 degrees off. The magnitude of the identified model stays quite constant among the 4, 6, and 9 m/s models. As the speed increases the Whipple and arm models are less predictive of the identified model at low frequencies, but tend to match out past 4 rad/s or so.

\(\frac{\phi}{T_\phi}\) frequency response of the three models, (I)dentified, (W)hipple, and (A)rm, at four speeds (2, 4, 6, and 9 m/s). The color indicates the model and the line type indicates the speed. Generated by src/systemidentification/canonical_plots.py.

Figure 11.14 shows the steer angle response with respect to the roll torque. Once again the Whipple model is almost identical to the identified model at the lowest speed, 2 m/s. The arm model has about a +5 dB offset at low frequencies and a -10 dB offset at high frequencies with the models only being similar just around 4 rad/s. At 9 m/s the Whipple model has similar magnitude and phase above 7 rad/s whereas the low frequency shows that the Whipple model has up to +30 dB magnitude difference with respect to the identified model. The arm model behaves in much the same way at 9 m/s.

\(\frac{\delta}{T_\phi}\) frequency response of the three models, (I)dentified, (W)hipple, and (A)rm, at four speeds (2, 4, 6, and 9 m/s). The color indicates the model and the line type indicates the speed. Generated by src/systemidentification/canonical_plots.py.

The steer torque to roll angle transfer function, Figure 11.15 may be the most important to model accurately as it is the primary method of controlling the bicycle’s direction, i.e. commanding roll allows one to command yaw. At 2 m/s the Whipple model magnitude matches at lower frequencies (< 4 rad/s) better and the arm model better at higher frequencies (> 4 rad/s). A 4 m/s and above the magnitude of identified varies little with speed, which contrasts the stronger speed dependence of the first principles models. The low frequency behavior of the identified model is not well predicted by the Whipple and arm models at the three highest speeds but about about 3 rad/s the arm model shows better magnitude matches than the Whipple model.

Frequency response of the three models, (I)dentified, (W)hipple, and (A)rm, at four speeds (2, 4, 6, and 9 m/s). The color indicates the model and the line type indicates the speed. Generated by src/systemidentification/canonical_plots.py.

The steer torque to steer angle frequency response, Figure 11.16, shows that at speeds above 2 m/s the first principle models do not predict the response well at low frequencies. The response changes more drastically with respect to speed for the first principles models than for the identified model.

Frequency response of the three models at four speeds. The color indicates the model and the line type indicates the speed. Generated by src/systemidentification/canonical_plots.py.

Comparison of identified models¶

Tables 11.3, 11.4, and 11.5 present the identified canonical coefficient values from all twelve of the chosen data subsets. The variances of the coefficient estimates with respect to how well the model fits the data in the steer equation are quite low except for the \(H_{\delta F}\) parameter. The low variance is partially due to the large datasets but also due to the quality of the resulting fits. The \(H_{\delta F}\) is highly dependent on the trail which is expected to be difficult to identify. The coefficient estimates are not all repeatable among the data subsets. In particular \(C_{1\delta\phi}\) and \(H_{\delta F}\) have the largest relative variance among the models (> 62%). \(M_{\delta\delta}\), \(C_{1\phi\delta}\), \(C_{1\delta\delta}\), and \(K_{2\delta\delta}\) are the most consistent among the models (< 18%) with the remaining coefficients falling somewhere in between. This is odd because these data sets are not independent of one another where some are subsets of the others but once again the identification procedure does not explicitly account for process noise and the coefficient estimates reflect this. A sensitivity study for each coefficient with respect to the various physical parameters could help reveal which parameters less likely fit into the Whipple model mold.

It is interesting to note that \(C_{1 \delta \phi}\) deviates quite significantly from the Whipple model prediction. This term depends on the wheel radii, wheel rotational inertia, wheelbase, steer axis tilt, and trail. All of these but trail are easily measured so it is tempting to solve for trail given \(C_{1 \delta \phi}\) and the other measured parameters described in [MPRS07].

The results for each data subset are given in Table 11.6. On average the values are extremely unrealistic when compared to the measured geometric trail of \(c=0.0599\) meters. This may imply that including the effects of pneumatic trail would not be enough to improve the predictive capabilities of the Whipple model. It also points to the potentially erroneous assumptions made in this identification effort such as the order of the system and which parameters are fixed. Even thought the model is quite adequate for describing the measured motion, it lacks the complexness to identify the precise deficients in the first principles assumptions.

The measurement errors, model structure and order dictate how well the models can predict the input-output behavior of each run or even each perturbation. The ideal goal is to select one or a small set of models that do a good job at predicting the measured behavior for each run. The previous section’s state space methods have already shown that the Whipple and arm models may not provide adequate predictions.

Figures 11.17, 11.18, and 11.19 plot the steer torque to roll angle frequency response for three speeds: 2 m/s, 5.5 m/s and 9.0 m/s for each of the models in Tables 11.3, 11.4, and 11.5. At the lowest speed, all of the models have a similar frequency response, especially in the frequency band between about 1 an 20 rad/s. At 5.5 m/s the models are similar at a higher bandwidth, 4 to 30 rad/s, and at 9.0 m/s even higher, 10 to 50 rad/s. Notice that the frequency band where the models are most similar shifts to higher frequencies at higher speeds. The model derived from all of the data (all rider and all runs), gives an average model and if this model is significantly better at predicting the measured behavior of the Whipple and arm models, it may be a good general candidate model for this bicycle.

Steer torque to roll angle frequency responses at 2.0 m/s for all the identified models in Tables 11.3, 11.4, and 11.5.

Steer torque to roll angle frequency responses at 5.5 m/s for all the identified models in Tables 11.3, 11.4, and 11.5.

Steer torque to roll angle frequency responses at 9.0 m/s for all the identified models in Tables 11.3, 11.4, and 11.5.

The predictive capability and quality of a given model can be quantified by an assortment of criterion and methods. I’ve made use of two criterion to judge the quality of these models with respect to given data. The first is to simulate the system given the measured inputs. This method works well when the open loop system is stable, but if it is unstable as so in the case of this bicycle, it may be difficult to simulate. Searching for initial conditions that give rise to a stable model for the duration of the run or simulating by weighting the future error less may relieve the instability issues. Another option is to see how well the inputs are predicted given the measured outputs. The canonical form of the equations lend themselves to this and only two inputs per run need be checked.

The predicted torques on the system are

and the measured torques as

\(y_p\) and \(y_m\) can be computed for each run along with the variance accounted for (VAF), by the model for both the total roll torque and the steer torque

I compute the VAF for each of the 374 runs used in the canonical identification outlined in Equation (29) using each of the 12 identified models and both the Whipple and Arm models. This percentage can be used as a criterion with which to judge the ability of one model versus another to predict the measurement. I then take the median of the VAF over each of the 12 sets of runs, Tables 11.7 and 11.8. The results give an idea of how well the various models are able to predict the data for all of the runs in a given set.

Tables 11.7 and 11.8 give the median for each set of runs in each column for each model given in the row for roll and steer respectively. The maximum VAF in the column gives a measure of the best model for predicting each individual run in that set of runs. Intuitively, I would expect that the diagonal of the upper 12 rows would be the maximum in each column due to the fact that that model was derived from that set of runs, but that is not always the case. I believe that this can be explained by the fact that there are more outlier runs in some sets. These outliers have enough effect in the resulting regressions, that the models generated from sets of runs with fewer outliers are able to predict the data better.

The models are able to predict the steer torque much better than the roll torque. The roll torque should be zero in all of the runs without disturbances but the roll equations do not predict a zero value. This is also reflected in the negative median values of all the runs with disturbances in much of Table 11.7. We fixed six of the nine parameters in the roll equation to those of the Whipple model and fixed three of the nine parameters in the steer equation. These extra degrees of freedom can partially explain why the steer predictions are better than the roll predictions. The model for the roll torque is more susceptible to the noise in the rate and angle measurements and has consequences of +/- 50 Nm variation in the predicted roll torques. These are unfortunately comparable in magnitude to the measured roll torques due to the lateral perturbations. But Table 11.7 can still be used to gage which models are better with reference to each other. The values in Table 11.7 are only generated from the runs with disturbances, as a relative measure of quality to zero is hard to make.

Tables 11.7 and 11.8 reveal:

• The arm model is poor at predicting the steer torque.
• The models derived from Charlie’s runs are poorer at predicting the inputs.
• The Whipple model is not too bad at predicting steer torque, but on average about 10% worse than the best models.
• The models identified from the pavilion runs are generally the best (with exception of Charlie’s). The ones generated from Luke and Jason’s runs are typically the best at predicting both steer torque and roll torque, with Luke’s giving better roll torque predictions.
• The roll torque is poorly predicted by all models when it is supposed to be zero. This raises implications in the validity of the roll equation and the potential need for tire slip models.

It may seem odd that a model identified from the subset of runs of one rider in one environment is the best at predicting the runs on a individual basis, but the uncertainty and error in both the data and the model structures don’t dictate that this can’t be. Keep in mind that all of the frequency responses of all 12 models shown in Figures 11.17, 11.18, and 11.19 are probably bounded in the uncertainty of the predicted responses and each can be considered a “good” model, even including the Whipple model.

The second method of evaluating the quality of the identified models is to simulate the model with the measured inputs and compare the predicted outputs with the measured outputs with a similar VAF criterion. I simulated all 14 models with the inputs from the 374 runs and computed the VAF explained by the model for each output. Since the models are typically unstable at all of the speeds we tested, I searched for the set of initial conditions which minimizes the VAF for all outputs. For the majority of runs and models, this is sufficient to find a stable simulation for the duration of the run. However, this is not always the case. For long duration runs I select a random 20 second section of the data to simulate, reducing the likelihood that the simulation blows up due to the model’s instability. Finally, I ignore any outputs VAFs that are less than -100 percent as they are most likely due to unstable simulations. Table 11.9 presents the median percent variance accounted for across all runs for each model and each output. The best model seems to be the one generated from the data with Luke on the Pavilion Floor once again, but these results differ from the previous otherwise.

• For all outputs other than roll angle, the arm model is better than the Whipple model.
• The models from Charlie’s data fare much better than the input comparisons and are better than some of Jason’s.
• The model identified from the data with Jason on the Pavilion floor is very poor in roll angle prediction as opposed to being a good choice from the input comparison results.
• All of the identified models are better predictors than the first principles models.

The mean percent variance across the outputs can be computed and the models ranked by the mean, Table 11.10. The best model seems to be LP and the AA is also a pretty good predictor. Notice that the Whipple model is poorer than the arm model.

The orange lines in Figures 11.8 and 11.9 correspond to the L-P model which allows comparison of the results of the canonical identification process with those of the state space identified models. The L-P model seems be better at fitting the data, especially in the steer acceleration equation, but the large variance in the state space coefficients is still a problem. This lends more confidence that that the L-P model is a better model choice than the Whipple or the arm model.

Discussion¶

Canonical identification

The canonical realization, Section Canonical Identification, is a good method for identifying a model for data from multiple runs. It relies on quality measurements of the coordinates, rates, and accelerations. I use numerical differentiation of the rates to get the accelerations instead of direct measurements and the angles are not perfectly related to the rates through differentiation because they were measured from different sensors. The noise in the measurements and whether the measurements are accurately the derivatives of one another have bearing on the results. It is possible to identify all of the entries in the canonical matrices, but it is likely some of those are easily predicted from first principles so I fix them to the Whipple model predictions if that is the case. This leaves the roll equation mostly known and the steer equation mostly unknown and the results reflect better fits with respect to steer than roll as a result. This formulation does not explicitly account for process noise, so it may be be susceptible to similar accuracy errors as the state space formulation is. The advantage in this method is the ability to use large sets of data for the calculations. It is extremely surprising that a model from a small subset of the data is better at predicting all of the runs on an individual basis.

Input comparison

The input prediction comparisons do not predict the roll torque well at all. It seems that the roll torque equation magnifies the noise in the coordinate, rate, and acceleration measurements such that the resulting noise in the roll torque is equivalent in magnitude to the roll perturbations. But the roll perturbations do seem to clearly be present in the predictions. This results in difficultly comparing the quality of the resulting models with respect to the roll equation. This also points to the potential deficiencies in the Whipple roll torque equation and these large magnitude roll torques may also be due to inaccurate modeling of the tire dynamics.

Output comparison

The output comparison (simulations) give more reasonable results because all four outputs generally fit well across runs given the measured inputs. It is surprising that the ranking of model prediction ability is different for the input comparisons than the output comparisons, but the fact that the model identified from Luke’s pavilion runs is the best from both comparisons, at least gives credence to its further adoption. The first principles models are dead last in the ranking and the model identified from Jason’s pavilion runs is surprisingly poor due to poor roll angle prediction.

Whipple model

The input comparisons show that the Whipple model is relatively reasonable at predicting the data but the output comparisons make it out to be much poorer. The Whipple model is clearly the worst at explaining the variance in the steer angle, roll rate, and steer rate outputs and is second to worst in roll angle. Also contrary to the Whipple model predictions, the weave mode of all identified models is unstable until at least 8-9 m/s or so. The caster mode is also typically much faster in the identified models. This implies that the real system does not benefit from open loop stability at all during most normal speed bicycling and that the rider is always responsible for stabilizing the vehicle. This may explain why none of the riders ever felt comfortable enough to try hands-free riding while strapped into the harness.

Arm model

I had hypothesized that the arm model would better predict the measured motion because it more accurately modeled the fact that we allowed the rider full use of his arms to control the steering and that the arms effectively added to the front frame inertia. This was validated with respect to the output simulation comparisons. They predict that the arm model is much better than the Whipple model for most of the output variables. But this is in contrast to the input comparison predictions which were the opposite, with the Whipple being much better than the arm model. More work is needed to verify which model is better and why the results differ at all.

Rider biomechanics may mot be modeled

The models identified from Charlie’s runs are different than those of Jason and Luke. The models from Charlie’s runs do not predict the runs very well, even including the subset of Charlie’s data. I’m not sure if the rider’s arm stiffness can affect these results or how much the different riders can effect this if we are only searching for the passive bicycle-rider model. The other potential explanation is that there are too many outliers in Charlie’s runs. This could have something to do with the runs that had time synchronization issues.

Predicting derivatives

The roll angle is more poorly predicted than the other variables. The steer angle and steer rate are better predicted than the roll angle and roll rate. In the output comparisons I enforce that the roll rate is the derivative of the roll and the same for the steer variables. The poorer prediction of roll angle is probably due to noise and error in the independent measurements of these variables. I toyed with complementary filters to try to combine the angle and rate measurements in a way that filtered and enforced the derivative relationship between the measured variables, but did not have much luck improving the results. It may be better to focus on one each of the roll and steer variables for minimization purposes. It is well known that fitting models with many fewer inputs to outputs is difficult and the fewer outputs reduces the number of noise terms to estimate.

Conclusion¶

The best candidate model for the measured system is the model identified from the data subset with rider Luke and the pavilion floor. I find no reason to use different models for each rider or environment as this model does a better job at predicting than the other models derived from matching subsets. I will use the model identified from the set of runs with Luke on the pavilion floor as the base bicycle model for identifying the controller for all the runs in the following section of this Chapter. I could use the individual bicycle identifications for each run, but using a model derived from a set of runs has the advantage that it will be less affected by the lack of identifying the process noise explicitly.

Suggestions for improving the results.

• Fit to MISO models instead of MIMO for simplification.
• Fix at least the \(a_{\phi\delta}\) coefficient and make the noise with respect to the kinematical equations equal to zero giving 11-13 total parameters to fit.
• Use model reduction techniques to combine many MISO and SISO models for a given run.
• Use better initial guess techniques and try global optimizers to get to the best solution.

Loop Shaping¶

The values of the parameters are difficult to directly compare with the ones found by the loop shaping technique described in Chapter Control due to some of the educated guesses about the rider’s internal control choices. But it is possible to evaluate the same loops’ frequency responses in an effort to understand how the rider chooses the gains.

Figure 11.25 shows the frequency response of the closed inner-most loop for a particular speed bin. The theory presented in [HMH12] and Chapter Control postulates that the rider chooses a gain such that the damping ratio of the high frequency neuromuscular peak around 10 rad/s and about gives a 10 dB peak (\(\zeta=0.15\)). Figure 11.25 shows that there may be a more heavily damped neuromuscular peak, but the large variability in the lower frequencies indicates that the rider’s choice is not so constant. This may be explained by the fact that it is more critical for the roll rate loop to exhibit the neuromuscular peaking and the differences in the three rider’s behavior. I’ve previously shown in Chapter Control that there are not necessarily unique gains in the first two loops to achieve the correct peaking in the roll rate loop.

Frequency responses of the closed steer angle loop around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid black lines give the mean gain and phase bounded by the dotted black lines which indicate one standard deviation.

Figure 11.26 shows the frequency response of the closed roll rate loop \(\frac{\dot{\phi}}{\dot{\phi}_c}(s)\). As was just pointed out, the theory is that this loop, which completes the inner control loop, must exhibit this typical neuromuscular peaking seen in human-machine control tasks. The mean magnitude plots indicates that this is the case and that the riders do choose their steer and roll rate gains such that the inner loops exhibit the typical response.

Frequency responses of the closed roll rate loop around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid black line give the mean gain and phase bounded by the dotted black lines which indicate the one sigma standard deviation.

Our theory for the selection of the remaining three gains suggests that the loops follow the dictates of the crossover model (i.e. 20 db slope around crossover) and that the crossover frequencies start at 2 and each successive loop is closed at half the previous crossover frequency such that the roll angle, heading and lateral deviation loops are closed in that order. Figures 11.27, 11.28, and 11.29 show the empirically derived frequency responses for the remaining loops.

Frequency responses of the open roll angle loop around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid blacks line give the mean gain and phase bounded by the dotted black lines which indicate one standard deviation.

Frequency responses of the open heading loop around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid blacks line give the mean gain and phase bounded by the dotted black lines which indicate one standard deviation.

Frequency responses of the open lateral deviation loop around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid black lines give the mean gain and phase bounded by the dotted black lines which indicate one standard deviation.

The median crossover frequencies for this 4 m/s speed bin are \(\omega_{\phi c}=3.66\), \(\omega_{\psi c}= 0.98\), and \(\omega_{y_q c}=1.06\). This indicates that the rider is more aggressive than our theory predicts with the higher roll crossover frequency in the roll angle loop. Also, the riders do not crossover at half the previous frequency, with the two outer loops crossing over at about the same frequency indicating that the lateral deviation tracking is more pertinent to the rider than the theory suggests. This may be partially due to the narrowness of the treadmill which constitutes 80% of the data.

Our hypothesis is that the crossover frequencies are constant with respect to speed, so the distribution of the crossover frequencies for all runs should have a tight distribution regardless of speed. Secondly, we postulated that the rider crosses over the roll angle loop around 2 rad/s. Figure 11.30 shows the distribution of the crossover frequencies for each loop where the median values of the crossover frequencies are \(w_{\phi c}=3.31\), \(w_{\psi c}= 1.20\), and \(w_{y_q c}=0.98\), all with somewhat large standard deviations.

Median crossover frequencies in rad/s for the open outer loops for all runs at all speeds.

If the medians are any indication of the true crossover frequency the rider is more aggressive than predicted and the two outer loop crossover frequencies are much closer in value than predicted. The one-half rule of thumb does seem to hold for the inner two loops.

Full System Response¶

Once the loops are all closed the system output response to the two inputs \(F_B\) and \(y_{qc}\) can be examined. Given a commanded input \(y_{qc}\), the tracking performance can be gauged by the closed loop Bode plot \(\frac{y_q}{y_{qc}}(s)\) given in Figure 11.31. Notice that in general, the response is well behaved with a 0 dB magnitude out to about 10 rad/s, all poles are stable and the phase lag increases gradually with frequency.

Frequency responses of the closed loop system lateral deviation tracking response around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid black lines give the mean gain and phase bounded by the dotted black lines which indicate one standard deviation.

Similarly, the system’s response to the lateral disturbance force Figure 11.32 is favorable for all the identified models and shows good disturbance rejection.

Frequency responses of the closed loop system disturbance rejection response around 4 m/s. The grey lines plot the response of each individual run in the speed bin while the solid black lines give the mean gain and phase bounded by the dotted black lines which indicate one standard deviation.

The reader may note that some of the Bode plots among Figures 11.25-11.29 indicate instability in that amplitude and phase characteristics synonymous with quadratic poles with small, positive damping ratios are in evidence. For example, one can see the evidence of such poles around 1-2 rad/s in Figure 11.25 and 10 rad/s in Figure 11.26. It should be emphasized that these Bode plots represent either (1) open-loop transfer functions, or (2) inner, closed loop transfer functions. Hence, they do not imply final, outer closed loop instabilities. Indeed, the Bode plots of the outer-most, closed loop transfer function of Figure 11.31 (lateral deviation) do not indicate any such quadratic poles with small positive damping ratios. The gain selection procedure in Chapter Control attempted to insure stability at least by closing of the roll angle loop, but these plots show that may not always be required for adequate performance for the complete system.