# Bicycle Balance Control Simulation

The model presented here is based on the model described here:

https://plot.ly/ipython-notebooks/bicycle-control-design/

## eval_bicycle_rhs

The file eval_bicycle_rhs.m encodes the 1st order ordinary differential equations that govern the steer-roll dynamics of a bicycle. The steer and steer rate inputs are set to constants.

function xdot = eval_bicycle_rhs(t, x, r, p)

% unpack states
omega = x(1);
theta = x(2);
psi = x(3);

% unpack inputs
beta = r(1); % steer rate
delta = r(2); % steer angle

% unpack parameters
m = p.m;
g = p.g;
h = p.h;
a = p.a;
b = p.b;
v = p.v;
I = p.I;

% calculate the state derivatives
omegadot = (m*g*h*theta - m*h/b*(a*v*beta + v^2*delta))/(I + m*h^2);
psidot = v/b*tan(delta);

% pack up the results

end


## eval_bicycle_rhs2

The file eval_bicycle_rhs2.m encodes the 1st order ordinary differential equations that govern the steer-roll dynamics of a bicycle. The steer and steer rate inputs are calculated by an input function.

function xdot = eval_bicycle_rhs2(t, x, w, p)

% unpack states
omega = x(1);
theta = x(2);
psi = x(3);

% calculate r with the passed in function w(t, x, p)
r = w(t, x, p);

% unpack inputs
beta = r(1); % steer rate
delta = r(2); % steer angle

% unpack parameters
m = p.m;
g = p.g;
h = p.h;
a = p.a;
b = p.b;
v = p.v;
I = p.I;

% calculate the state derivatives
omegadot = (m*g*h*theta - m*h/b*(a*v*beta + v^2*delta))/(I + m*h^2);
psidot = v/b*tan(delta);

% pack up the results

end


## eval_steer_control_input.m

The file eval_steer_control_input.m encodes a PD controller that generates the steer and steer rate inputs.

function r = eval_steer_control_input(t, x, p)

% unpack the states
omega = x(1);
theta = x(2);

% unpack the control gains
kp = p.kp;

% calculate the inputs
beta = kp*omega;
delta = kp*theta;

% pack up the inputs
r = [beta; delta];

end


## bicycle_main.m

The file bicycle_main.m runs the simulation with the constant inputs and the inputs from the controller.

% this is the simulation we did in class of the bicycle balance model

% Set the initial values of the states (initial conditions).
% Given an initial 5 deg roll angle and 5 deg steer angle, simulate the system.
x0 = [0;  % roll rate [rad/s]

% Generate time values for a 4 second duration.
ts = linspace(0, 4, 200);  % time values

% Set the inputs to constant values.
r = [0;  % steer rate [rad/s]

% Define all of the numerical values for the constant parameters in a structure.
p.m = 87;  % kg
p.g = 9.81;  % m/s^2
p.h = 1.0;  % m
p.a = 0.5;  % m
p.b = 1.0;  % m
p.v = 5.0;  % m/s
p.I = 3.28;  % kg m^2

% Setup the anonymous function for ode45.
eval_bicycle_rhs_anon = @(t, x) eval_bicycle_rhs(t, x, r, p);

% Integrate the ODEs (simulate) with ode45.
[ts, xs] = ode45(eval_bicycle_rhs_anon, ts, x0);

% plot the results
figure(1)
plot(ts, xs*180/pi)
xlabel('Time [s]')
legend('\omega [deg/s]', '\theta [deg]', '\psi [deg]')

% this is a second simulation of the same model but it has an input function
% with a controller

% add a parameters for the proportional control input
p.kp = 2.5; % proportional gain

% define a new anonymous function that makes use of the input function
eval_bicycle_rhs_anon2 = ...
@(t, x) eval_bicycle_rhs2(t, x, @eval_steer_control_input, p);

% simulate
[ts, xs] = ode45(eval_bicycle_rhs_anon2, ts, x0);

% plot the state trajectories
figure(2)
plot(ts, xs*180/pi)
xlabel('Time [s]')
legend('\omega [deg/s]', '\theta [deg]', '\psi [deg]')

% plot the input trajectories
rs = zeros(length(ts), 2);
for i = 1:length(ts)
rs(i, :) = eval_steer_control_input(ts(i), xs(i, :), p);
end
figure(3)
plot(ts, rs*180/pi)
xlabel('Time [s]')
legend('\beta [deg/s]', '\delta [deg]')