The code is quite messy, with numerous unformatted local functions scattered throughout, making it visually challenging to troubleshoot. Additionally, there are some typographical errors in the equations.
However, I spent several hours rectifying the code's core structure, resulting in the ode45 solver running without errors. Nonetheless, you still need to properly rewrite the nested functions responsible for calculating the forces on the wheels and tires, both front and rear, on both the right and left sides.
Best of luck!
tspan = [0 10];
s0 = zeros(1, 16);
[t, s] = ode45(@ODE_MO, tspan, s0);
% Plotting state variables, s1 to s16
figure(1)
for j = 1:length(s0)
subplot(4,4,j)
plot(t, s(:,j)), grid on
str = sprintf('s_{%1d}', j);
title(str, 'FontSize', 9);
end
%% 8-DOF Car Model
function dsdt = ODE_MO(t, s)
% Model Parameters
Meo = 0.85;
a = 1.203;
Bf = 1.33;
Br = 1.33;
b = 1.217;
Cphi = 3511.6;
Kphi = 66185.8;
Ccfr = 30000;
Ccfl = 30000;
Ccrr = 30000;
Ccrl = 30000;
Cxfr = 50000;
Cxfl = 50000;
Cxrr = 50000;
Cxrl = 50000;
e = 0.2;
g = 10;
h = 0.5;
Ixx = 750;
Ixz = 0;
Iw = 2.1;
Iz = 2500;
l = 2.42;
m = 1280;
ms = 1160;
rw = 0.3;
% Inputs Section
Ster_Ang_Fr = 30;
Ster_Ang_Fl = 30;
Ster_Ang_rr = 0;
Ster_Ang_rl = 0;
T_FR = 0;
T_FL = 0;
T_RR = 150;
T_RL = 150; % <-- fixed
% Definitions of states
x = s(1);
xdot = s(2); % Vx
y = s(3);
ydot = s(4); % Vy
psi = s(5);
psidot = s(6); % Psi
phi = s(7);
phidot = s(8); % R
% Decoupling ψ" & ϕ" by getting the inverse of Inertia matrix
f1 = a*(Fyfr + Fyfl) - b*(Fyrr + Fyrl) + Bf/2*(Fxfl - Fxfr) + Br/2*(Fxrl - Fxrr);
f2 = - ms*e*(Fyfr + Fyfl + Fyrr + Fyrl - m*xdot*psidot)/m - ms*e*xdot*psidot + ms*g*e*sin(phi) - Kphi*phi - Cphi*phidot;
MI = [Iz -Ixx; -Ixz Ixx-(ms*e*ms*e/m)]; % Iz*psiddot - Ixx*phiddot = f1;
F = [f1; f2]; % - Ixz*psiddot + Ixx*phiddot - (ms*e*ms*e/m)*phiddot = f2;
fcn = MI\F;
psiddot = fcn(1); % ψ"
phiddot = fcn(2); % ϕ"
xddot = (Fxfr + Fxfl + Fxrr + Fxrl - Faero - Froll + m*ydot*psidot)/m;
yddot = (Fyfr + Fyfl + Fyrr + Fyrl - ms*e*phiddot - m*xdot*psidot)/m;
%% Equations of Motion
dsdt = zeros(16, 1); % initialization
%% Car dynamics
% Longitudinal motion
dsdt(1) = s(2); % x'
dsdt(2) = xddot;
% Lateral motion
dsdt(3) = s(4); % y'
dsdt(4) = yddot;
% Yaw motion
dsdt(5) = s(6); % psi'
dsdt(6) = psiddot;
% Roll motion
dsdt(7) = s(8); % phi'
dsdt(8) = phiddot;
%% Wheel dynamics
% Front right wheel
dsdt(9) = s(10);
dsdt(10) = (T_FR - rw*Ftfr)/Iw;
% Front left wheel
dsdt(11) = s(12);
dsdt(12) = (T_FL - rw*Ftfl)/Iw;
% Rear right wheel
dsdt(13) = s(14);
dsdt(14) = (T_RR - rw*Ftrr)/Iw;
% Rear left wheel
dsdt(15) = s(16);
dsdt(16) = (T_RL - rw*Ftrl)/Iw;
%% Nested Functions
function out = Fxfr()
out = 1;
end
function out = Fxfl()
out = 1;
end
function out = Fxrr()
out = 1;
end
function out = Fxrl()
out = 1;
end
function out = Faero()
out = 20;
end
function out = Froll()
out = 1;
end
function out = Fyfr()
out = 1;
end
function out = Fyfl()
out = 1;
end
function out = Fyrr()
out = 1;
end
function out = Fyrl()
out = 1;
end
function out = Ftfr()
out = 1;
end
function out = Ftfl()
out = 1;
end
function out = Ftrr()
out = 1;
end
function out = Ftrl()
out = 1;
end
end
