Learning Objectives

After completing this lab you will be able to:

• formulate the explicit first order ordinary differential equations for the roll and steer dynamics of a bicycle
• translate ordinary differential equations into a computer function that evaluates the equations at any given point in time
• develop a function that evaluates state dependent input functions
• develop a function that evaluates state and input dependent output variables
• numerically integrate ordinary differential equations with Octave/Matlab's ode45
• create complete and legible plots of the resulting input, state, and output trajectories
• create a report with textual explanations and plots of the simulation
• evaluate a bicycle model for stability at different speeds
• design a feedback controller that stabilizes a bicycle at different speeds
• design a bicycle that is self-stable at low speeds

Introduction

In this lab, you will investigate the dynamics of a bicycle with and without a controller to assist in balancing. Bicycles and motorcycles share the same core dynamic behavior, so this lab will introduce you to phenomena that apply to both vehicles. Below is a video that demonstrates a bicycle's uncanny ability to be self-stable, i.e. return to an upright configuration even when perturbed from its equilibrium. This is, of course, not always possible as we know that a bicycle will simply fall over if it has no forward speed and we also know that bicycles require some control from the rider in general. In this lab, you will simulate the Whipple-Carvallo bicycle model, investigate its dynamics, apply a controller to stabilize it at low speeds, and design a bicycle that can self-stabilize at lower speeds.

System Description

The Whipple-Carvallo bicycle model is essentially the simplest model that behaves like a bicycle and can exhibit self-stability. "Self-stability" is defined as open loop stability or stability without any control inputs. The model was developed by both Whipple and Carvallo around the turn of the 20th century. In 2007, Meijaard, Papdopoulos, Ruina, and Schwab validated the model and presented a benchmark formulation [Meijaard2007]. You'll be using this formulation for the lab. This model assumes no-slip non-deforming tires, a flat level surface, frictionless revolute joints, and a human rider that is rigidly fixed to the rear bicycle frame. See the cited paper for more details about the history of the model, its derivation, and in-depth analysis. Figure 1 shows the basic geometry and the four rigid bodies: [B] rear frame and human body, [H] handlebar and fork, [R] rear wheel, and [F] front wheel.

Figure 1: Schematic of the benchmark bicycle model.

function [M, C1, K0, K2] = compute_benchmark_bicycle_matrices(p)
% COMPUTE_BENCHMARK_BICYCLE_MATRICES - Returns the canonical matrices of the
% benchmark Whipple-Carvallo bicycle model linearized about the upright
% constant velocity configuration. It uses the parameter definitions from:
%
% J.P Meijaard, Jim M Papadopoulos, Andy Ruina and A.L Schwab, 2007,
% Linearized dynamics equations for the balance and steer of a bicycle: a
% benchmark and review, Proc. R. Soc. A. 4631955–1982
% http://doi.org/10.1098/rspa.2007.1857
%
% Syntax: [M, C1, K0, K2] = compute_benchmark_bicycle_matrices(p)
%
% Inputs:
% p - structure
%     A dictionary of the benchmark bicycle parameters. Make sure your units
%     are correct and match the benchmark paper's units.
%
% Outputs:
% M - matrix, size 2x2
%     The mass matrix.
% C1 - matrix, size 2x2
%     The part of the damping matrix that is proportional to the speed, v.
% K0 - matrix, size 2x2
%     The part of the stiffness matrix proportional to gravity, g.
% K2 - matrix, size 2x2
%     The part of the stiffness matrix proportional to the speed squared,
%     v^2.

mT = p.mR + p.mB + p.mH + p.mF;
xT = (p.xB*p.mB + p.xH*p.mH + p.w*p.mF)/mT;
zT = (-p.rR*p.mR + p.zB*p.mB + p.zH*p.mH - p.rF*p.mF)/mT;

ITxx = (p.IRxx + p.IBxx + p.IHxx + p.IFxx + p.mR*p.rR^2 + p.mB*p.zB^2 + ...
        p.mH*p.zH^2 + p.mF*p.rF^2);
ITxz = (p.IBxz + p.IHxz - p.mB*p.xB*p.zB - p.mH*p.xH*p.zH + p.mF*p.w*p.rF);
p.IRzz = p.IRxx;
p.IFzz = p.IFxx;
ITzz = (p.IRzz + p.IBzz + p.IHzz + p.IFzz + p.mB*p.xB^2 + p.mH*p.xH^2 + ...
        p.mF*p.w^2);

mA = p.mH + p.mF;
xA = (p.xH*p.mH + p.w*p.mF)/mA;
zA = (p.zH*p.mH - p.rF*p.mF)/mA;

IAxx = (p.IHxx + p.IFxx + p.mH*(p.zH - zA)^2 + p.mF*(p.rF + zA)^2);
IAxz = (p.IHxz - p.mH*(p.xH - xA)*(p.zH - zA) + p.mF*(p.w - xA)*(p.rF + ...
        zA));
IAzz = (p.IHzz + p.IFzz + p.mH*(p.xH - xA)^2 + p.mF*(p.w - xA)^2);
uA = (xA - p.w - p.c)*cos(p.lambda) - zA*sin(p.lambda);
IAll = (mA*uA^2 + IAxx*sin(p.lambda)^2 + 2*IAxz*sin(p.lambda)*cos(p.lambda) + ...
        IAzz*cos(p.lambda)^2);
IAlx = (-mA*uA*zA + IAxx*sin(p.lambda) + IAxz*cos(p.lambda));
IAlz = (mA*uA*xA + IAxz*sin(p.lambda) + IAzz*cos(p.lambda));

mu = p.c/p.w*cos(p.lambda);

SR = p.IRyy/p.rR;
SF = p.IFyy/p.rF;
ST = SR + SF;
SA = mA*uA + mu*mT*xT;

Mpp = ITxx;
Mpd = IAlx + mu * ITxz;
Mdp = Mpd;
Mdd = IAll + 2*mu*IAlz + mu^2*ITzz;
M = [Mpp, Mpd; Mdp, Mdd];

K0pp = mT*zT;
K0pd = -SA;
K0dp = K0pd;
K0dd = -SA*sin(p.lambda);
K0 = [K0pp, K0pd; K0dp, K0dd];

K2pp = 0.0;
K2pd = (ST - mT*zT)/p.w*cos(p.lambda);
K2dp = 0.0;
K2dd = (SA + SF*sin(p.lambda))/p.w*cos(p.lambda);
K2 = [K2pp, K2pd; K2dp, K2dd];

C1pp = 0.0;
C1pd = (mu*ST + SF*cos(p.lambda) + ITxz/p.w*cos(p.lambda) - mu*mT*zT);
C1dp = -(mu*ST + SF*cos(p.lambda));
C1dd = (IAlz/p.w*cos(p.lambda) + mu*(SA + ITzz / p.w*cos(p.lambda)));
C1 = [C1pp, C1pd; C1dp, C1dd];

end

Initial Conditions

For the simulations, set the initial conditions as:

This will start the vehicle off with a small roll and steer angle. See Integrating the Equations for how to set up the initial condition vector. Make sure that your initial conditions are arranged in the same order as your state variables.

Time Steps

Simulate the system for 5 seconds with time steps of 1/100th of a second.

Eigenvalues and Stability

Matlab and Octave can easily calculate the roots of the characteristic equation, i.e. the eigenvalues, given a state dynamics matrix \(\mathbf{A}\). The eig() function calculates the eigenvalues. It is used like so:

evals = eig(A);

The variable evals will then contain a vector of eigenvalues. The functions real() and imag() are also helpful to extract the real and imaginary parts of the complex eigenvalues.

See the Matlab documentation for more about these functions:

• https://www.mathworks.com/help/matlab/ref/eig.html
• https://www.mathworks.com/help/matlab/ref/real.html
• https://www.mathworks.com/help/matlab/ref/imag.html

plot(speeds, real_part_of_evals, '.')

Using a Root Locus To Choose a Feedback Gain

A root locus is a parametric plot of the complex plane that shows how the roots of the characteristic equation change with respect to changes in a single parameter in a linear model. Any parameter can be chosen but the most common parameter is a proportional control gain of a closed loop negative feedback controller. It turns out that the Whipple-Carvallo bicycle model can be stabilized by "steering into the fall". One way to formalize this statement is to apply steer torque in the direction of the roll angular rate, i.e. if you roll towards the right, steer towards the right. You'll first need to determine the transfer function that relates the steer torque input \(T_\delta\) to the roll angular rate \(\omega\):

You can use Cramer's rule to find this transfer function, as in the book, or there are various Matlab/Octave functions that can help, for example ss2tf() can convert a state space system to equivalent transfer functions. If you use Cramer's rule you'll need to use the tf() function to construct the transfer function in Matlab/Octave. There isn't any reason to do any algebra by hand, use Matlab/Octave to do the number crunching.

Once you have a transfer function you can use the rlocus() function to make a standard root locus plot as a function of the negative feedback gain \(k\). You'll need to examine this plot to determine what gain value ensures that all roots are in the left half of the complex plane.

• https://www.mathworks.com/help/matlab/ref/ss2tf.html
• https://www.mathworks.com/help/matlab/ref/tf.html
• https://www.mathworks.com/help/control/ref/tf.rlocus.html

Checking your work

You can use the BicycleParameters web app to visualize the bicycle geometry and to double check your eigenvalue plots. Visit the app here:

http://bicycleparameters.herokuapp.com/

This app was recently developed by a UCD physics undergraduate, Noah Sanders, and is powered by similar calculations that you are doing in this lab. It is still in beta and there may be some bugs, but you should be able to edit the parameters in the table to see changes in the two figures.

Deliverables

In your lab report, show your work for creating and evaluating the simulation model. Include any calculations you had to do, for example those for state equations, initial conditions, input equations, time parameters, and any other parameters. Additionally, provide the indicated plots and answer the questions below. Append a copy of your Matlab/Octave code to the end of the report. The report should follow the report template and guidelines.

Submit a report as a single PDF file to Canvas by the due date that addresses the following items:

1. Create a function defined in an m-file that evaluates the right hand side of the ODEs, i.e. evaluates the state derivatives. See Defining the State Derivative Function for an explanation.

2. Create two functions defined each in an m-file that calculates the two requested inputs: no control and with control. See Time Varying Inputs for an explanation.

3. Create a function defined in an m-file that calculates the requested outputs. See Outputs Other Than the States and Outputs Involving State Derivatives for an explanation.

4. Create a script in an m-file that utilizes the above functions to simulate system for the scenarios described below. This should setup the constants, integrate the dynamics equations, and plot each state, and output versus time. See Integrating the Equations for an explanation.

5. Make a plot where the Y axis is the real part of the eigenvalues (determines growth or decay and the rate) and the X axis is the bicycle speed. The speeds should vary from 0 to 10 m/s (0 to ~22 mph). Discuss how the bicycle's self-stability relates to the speed and discuss any speed values that stand out and what the eigenvalues at this speed tell you about the motion. Use simulation (states vs time: 0 to 5 seconds) plots to back up your explanations at the particular speeds of interest. Be sure to explain your simulation plots in terms of growth/decay rates, oscillation frequencies, and phase differences in the states. Limit your Y axis to +/- 90 degs for unstable simulations.

6. Develop two roll feedback controllers using a root locus to select an appropriate gain to stabilize the bicycle at 2 m/s and 8 m/s. At 2 m/s feedback roll rate. At 8 m/s feedback roll angle. Both controllers should control steer torque. Choose gains that give realistic decay rates and oscillation frequencies (a human's arm control bandwidth has an upper bound of about 10 rad/s). Simulate the system at these speeds with the controller defined in an input function to show that the system is stable. Discuss the sign of the gain values.

7. Design a bicycle by changing any number of the parameters values so that it is self-stable (no control) at a very low range of speeds. Try not to choose parameter values that would be impossible to actually build, e.g. a wheelbase that is 100 meters long or rear frame/human mass of 1 kg. There is no single solution to this. Once you chose a bike, simulate it at a very low speed and discuss the simulation results and the advantages/disadvantages of the bicycle design. A practical use of a bicycle that is self-stable at low speeds could be to make it easy for people to learn to ride a bicycle. The videos below show bicycles designed for low-speed stability:

Assessment Rubric