Implementation of the Bicycle Model in Matlab

Bicycle model

The Bicycle Model is a simplified vehicle handling model, where the left and right tyres are combined into one at each axle. Despite the simplicity of the model, the output of the model is accurate enough to introduce the concept of the understeering/oversteering vehicles. The model is linearised – actual vehicles are non-linear systems, however, they can be accurately approximated for a limited range of their operations.

The 2DOF of the model are the lateral velocity \( v \) and yaw velocity \(r\). The steer angle \(\delta\) is the input, and the forward velocity \(u\) is specified.
In the model, no longitudinal and lateral load transfers are assumed. Body roll and pitching are ignored. Tyres are assumed to be in the linear range.
Therefore, the model is accurate for low lateral acceleration \(a_y\).

(1) Model equations

\[ \left\{ \begin{array}{l} \mathrm{mV}(\mathrm{r}+\dot{\beta}) &=\mathrm{Y}_{\beta} \beta+\mathrm{Y}_{\mathrm{r}} \mathrm{r}+\mathrm{Y}_{\delta} \delta &= (C_F + C_R)\beta + \dfrac{1}{V}(aC_F \; – bC_R) r + (-C_F) \delta \\ \mathrm{I}_{\mathrm{z}} \dot{\mathrm{r}} &=\mathrm{N}_{\beta} \beta+\mathrm{N}_{\mathrm{r}} \mathrm{r}+\mathrm{N}_{\delta} \delta &= (a C_F \; – b C_R)\beta + \dfrac{1}{V}(a^2C_F + b^2C_R) r + (-a C_F)\delta \end{array} \right. \]

State-space model

Linear system models can be represented in state-space. \[ \left\{ \begin{array}{l} \begin{align} \dot{\textbf{x}} &= A \textbf{x} + B \textbf{u} \\ \textbf{y} &= C \textbf{x} + D \textbf{u} \\ \end{align} \end{array} \right. \] where \(\textbf{x}\) is the state vector, \(\textbf{u}\) the input vector and \(\textbf{y}\) the output vector.

Bicycle model in state-space

For the bicycle model, choose the system states \[ \textbf{x}= \begin{bmatrix} x_1 \\ x_2 \\ \end{bmatrix} = \begin{bmatrix} \beta\\ r \\ \end{bmatrix} \] where \(\beta\) is the vehicle sideslip angle in radians, and \(r \) the yaw rate in rad/s.
The input is the steer angle \(\delta\) in radians \[ \textbf{u} = \begin{bmatrix} u_1 \end{bmatrix} = \begin{bmatrix} \delta \end{bmatrix} \] Setup the outputs \[ \textbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} v \\ r \\ a_y \end{bmatrix} \] where \(v\) is the sideslip velocity in m/s, \(r\) the yaw rate in rad/s and \(a_y\) the lateral acceleration in m/s^2

From the bicycle model equations, the model in state-space form is as follows \[ \begin{array}{l} \begin{align} \begin{bmatrix} \dot{\beta} \\ \dot{r} \end{bmatrix} &= \begin{bmatrix} \dfrac{Y_\beta}{mV} & \dfrac{Y_r}{mV} -1 \\ \dfrac{N_\beta}{I_z} & \dfrac{N_r}{I_z} \end{bmatrix} \begin{bmatrix} \beta\\ r \\ \end{bmatrix} + \begin{bmatrix} \dfrac{Y_\delta}{mV}\\ \dfrac{N_\delta}{I_z} \\ \end{bmatrix} \delta \\ \begin{bmatrix} v \\ r \\ a_y \end{bmatrix} &= \begin{bmatrix} V & 0\\ 0 & 1 \\ \dfrac{Y_\beta}{m} & \dfrac{Y_r}{m} \end{bmatrix} \begin{bmatrix} \beta\\ r \\ \end{bmatrix} + \begin{bmatrix} 0\\ 0\\ \dfrac{Y_\delta}{m} \end{bmatrix} \delta \end{align} \end{array} \] where the sideslip velocity \(v = V \tan{\beta} \approx V\beta \; \) for a small angle \(\beta\)
and the lateral acceleration \(a_y = V(r + \dot{\beta}) \)

Example model parameter

(To be written …)

Equivalent damping coefficient \(c\) and spring stiffness \(k\)

The eigenvalues \(\lambda\) of the \(A\) matrix are the roots \(s\) of the characteristic equation.
The characteristic equation is \[ \begin{align} \lvert \textbf{A} – \lambda \textbf{I} \rvert &= \begin{vmatrix} \dfrac{Y_\beta}{mV} – \lambda & \dfrac{Y_r}{mV} -1 \\ \dfrac{N_\beta}{I_z} & \dfrac{N_r}{I_z} – \lambda \end{vmatrix} \\ &= \lambda^2 + (-\dfrac{N_r}{I_z} – \dfrac{Y_\beta}{mV}) \lambda + (\dfrac{N_\beta}{I_z} +\dfrac{Y_\beta N_r – Y_r N_\beta}{mV I_{z}}) = 0 \\ \end{align} \] \[ \therefore I_z \lambda^2 + (-{N_r} – \dfrac{I_z Y_\beta}{mV}) \lambda + (N_\beta +\dfrac{Y_\beta N_r – Y_r N_\beta}{mV}) = 0 \] This is analogous to the characteristic equation of the mass-spring-damper system \[ ms^2+cs+k = m(s^2 + 2\zeta\omega_{n}s + \omega_{n}^2) =0 \] i.e.

For the bicycle model, the undamped natural frequency \(\omega_{n}\), damped natural frequency \(\omega_{d}\) and damping ratio \(\zeta\) can be calculated from the equivalent damping coefficient \(c_{eq}\) and the spring stiffness \(k_{eq}\) \[ \begin{array}{l} \omega_{n} &= \sqrt{\dfrac{k_{eq}}{I_z}} = \sqrt{\dfrac{N_\beta}{I_z} +\dfrac{Y_\beta N_r – Y_r N_\beta}{mVI_z}}\\ \zeta &= \dfrac{c_{eq}}{2I_z\omega_n} = \dfrac{c_{eq}}{2 \sqrt{k_{eq}I_z}} = \dfrac{-{N_r} \; – \dfrac{I_z Y_\beta}{mV}}{2\sqrt{(N_\beta + \dfrac{Y_\beta N_r – Y_r N_\beta}{mV} ) I_z}} \\ \omega_{d} &= \omega_{n}\sqrt{1-\zeta^2}\\ \end{array} \]

(2) Matlab implementation

The built-in state-space block in Simulink is useful for the numerical time integration of the system. Below is a simple Simulink model with the state-space block.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

clear all
% System states
% x1 = beta (vehicle slip angle [rad]) 
% x2 = r (yaw rate [rad/s])
 
% Input
% u1 = delta (steer angle [rad])
delta_deg = 10;
 
% Outputs
% y1 = v  (sideslip velocity [m/s])
% y2 = r  (yaw rate [rad/s])
% y3 = ay (lateral acceleration [m/s^2]) 
 
% Boundary condition
V = 10;    % forward velocity [m/s]
 
% Initial system states
x1_0 = 0; % beta at t = 0 [rad]
x2_0 = 0; % r    at t = 0 [rad/s]
 
 
% Vehicle parameters
Cf = -1020; % front tyre cornering stiffness per tyre [N/deg]
Cr = -760;  % rear tyre cornering stiffness per tyre [N/deg] 
mf = 620;  % mass at front axle [kg]
mr = 430;  % mass at rear axle [kg]
Iz = 1560; % yaw moment of inertia [kg m^2]
l  = 2.4;  % wheelbase [m]
 
% Conversion of units
Cf2 = Cf*2*180/pi; % front cornering siffness per axle [N/rad]
Cr2 = Cr*2*180/pi; % front cornering siffness per axle [N/rad]
m = mf + mr;       % total mass of vehicle [kg]
a = l * mr/m;      % distance between c.g. and front axle [m]
b = l * mf/m;      % distance between c.g. and rear axle [m]
 
% Bicycle model derivatives
Ybeta  = Cf2 + Cr2;
Yr     = (a*Cf2 - b*Cr2)/V;
Ydelta = -Cf2;
Nbeta  = a*Cf2 - b*Cr2;
Nr     = (a^2*Cf2 + b^2*Cr2)/V;
Ndelta = -a*Cf2;
 
% State-space matrices
A = [Ybeta/(m*V), (Yr/(m*V) - 1)
    Nbeta/Iz, Nr/Iz];
 
B = [Ydelta/(m*V)
    Ndelta/Iz];
 
C = [V, 0
    0, 1
    Ybeta/m, Yr/m];
 
D = [0
    0
    Ydelta/m];

(Continued …)

Reference

[1] W. Milliken and D. Milliken, Race Car Vehicle Dynamics.