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);
thetadot = omega;
psidot = v/b*tan(delta);
% pack up the results
xdot = [omegadot; thetadot; psidot];
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);
thetadot = omega;
psidot = v/b*tan(delta);
% pack up the results
xdot = [omegadot; thetadot; psidot];
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]
5*pi/180; % roll angle [rad]
0]; % heading angle [rad]
% 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]
5*pi/180]; % steer angle [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]')