How to code algebric loop in MATLAB script file?

Hello. I have some state space equations which should be solved with ode45 to plot x vs t. But there is an algebric loop inside equations and two equations are dependent to each other. I mean inside d1hat there is p and inside pdot (which should be integrated and I have defined it as the 5th equation of the states), there is d1hat. So whenever I write d1hat first p is not defined before that and vice versa. The code is as below

This is the road profile

function zr = fcn5(t)
%   times
t1 = t - 3.5;
t2 = t - 6.5;
t3 = t - 8.5;
t4 = t - 11.5;
%   d(t)
d = 0.002*sin(2*pi*t) + 0.002*sin(7.5*pi*t);
%   eq 45
if (t >= 3.5 && t < 5) 
    zr = -0.0592*t1^3 + 0.1332*t1^2 + d;
elseif (t >= 5 && t < 6.5)
    zr = 0.0592*t2^3 + 0.1332*t2^2 + d;
elseif (t >= 8.5 && t < 10)
    zr = 0.0592*t3^3 - 0.1332*t3^2 + d;
elseif (t >= 10 && t < 11.5)
    zr = -0.0592*t4^3 - 0.1332*t4^2 + d;
else
    zr = d;
end
end
%----------
%this is derivative of it 
function zr_dot = roadproder(t)
t1 = t - 3.5;
t2 = t - 6.5;
t3 = t - 8.5;
t4 = t - 11.5;
%  derivative of d(t)
d = 2*pi*0.002*cos(2*pi*t) + 7.5*pi*0.002*cos(7.5*pi*t);
%   eq 45
if (t >= 3.5 && t < 5) 
    zr_dot = 3*(-0.0592*t1^2) + 2*(0.1332*t1) + d;
elseif (t >= 5 && t < 6.5)
    zr_dot = 3*(0.0592*t2^2) + 2*(0.1332*t2) + d;
elseif (t >= 8.5 && t < 10)
    zr_dot = 3*(0.0592*t3^2) - 2*(0.1332*t3)+ d;
elseif (t >= 10 && t < 11.5)
    zr_dot = 3*(-0.0592*t4^2)- 2*(0.1332*t4)+ d;
else
    zr_dot = d;
end
end
%-----
%and this is the main function which will be solved with ode 45
function [dx,N,v] = States (t,x)
% figure(3)
% plot(t,dx(2,1))
% figure(4)
% plot(t,v)
%% parameters
ksl = 14500;
ksnl = 160000;
k = 1.8;
bl = -3;
m = k;
br = 2;
mus = 59;
ms = 390;
ms0 = 350;
w = ms0*900;
bs0 = 1250;
ks0 = 13500;
k1 = 20;
k0 = 1;
S = 2;
bsl = 1385;
bsnl = 524;
kt = 190000;
bt = 14500;
g = 9.81;

fs = ksl*(x(1)-x(3))+ksnl*(x(1)-x(3))^3;
fd = bsl*(x(2)-x(4))+bsnl*(x(2)-x(4))^2;
zr = fcn5(t);
fst = kt*(x(3)-zr);
zr_dot = roadproder(t);
fdt = bt*(x(4)-zr_dot);
ft = 0;
if (kt*(x(3) - zr) + bt*(x(4) - zr_dot)) < (ms + mus)*g
    ft = fst + fdt;
end
%% sliding surface
Sigma = S*x(1)+x(2);
%% MY PROBLEM IS IN THESE TWO LINES
d1hat = p + w*Sigma; 
pdot = -w*(S*x(2)-1/ms0*(ks0*x(1)+bs0*x(2))+d1hat+k0/ms0*v); 
veq = -ms0/k0*(S*x(2)+k1*Sigma-1/ms0*(ks0*x(1))+bs0*x(2)); 
%%

dx = zeros(5,1);
dx(1,1) = x(2);
dx(3,1) = x(4);
d1hat = dx(5,1)+ w*Sigma;
dx(5,1) = pdot; 
vn = -ms0/k0*d1hat;
v = veq + vn;
%% nonideal actuator
if v>= br
    dd = -m*br;
elseif v>bl && v<br
    dd = -m*v;
elseif v<= bl
    dd = -m*bl;
end
N = m*v + dd;
dx(2,1) = 1/ms*(-fs-fd+N);
dx(4,1) = 1/mus*(fs + fd - ft - N);
end

[Dyuman Joshi]

What are the equations you are trying to solve? Please provide the relationship of the equations in mathematical form.

And how do you call the main function via ode45?

[Faezeh Yousefpour]

Sorry for the late reply.I am trying to implement this article : https://www.sciencedirect.com/science/article/abs/pii/S0022460X18307491

"Linear disturbance based sliding mode control active suspension systems with non idealactuator". And the equations are number 1 to 4 in the article. I am calling the main function like this: ( @States5, [0 15], [0 0 0 0 0])

Answer 1

[Sam Chak]

Hi @Faezeh Yousefpour

I've got the code to run correctly from a mathematical perspective. However, it's essential to examine how the observer can estimate the lumped disturbance. Varying the parameter values can significantly affect control performance. While I'm not an expert in suspension systems, it's crucial to grasp the physics behind selecting these values rather than relying solely on optimization or reinforcement learning algorithms to determine them for you. Otherwise, you may find it challenging to provide a technical explanation or justification for your choices in journal papers, or answering questions in conferences, or even defending your dissertation.

tspan  = [0 15];
x0     = [0 0 0 0 0];
[t, x] = ode15s(@odefcn, tspan, x0);
plot(t, x(:,1)), grid on
xlabel('t'), ylabel('x_{1}')

%% Quarter-Car Suspension System
function [dxdt, N, v] = odefcn(t, x)
    %% parameters
    ksl    = 14500;
    ksnl   = 160000;
    k      = 1.8;
    bl     = -3;
    m      = k;
    br     = 2;
    mus    = 59;
    ms     = 390;
    ms0    = 350;
    w      = ms0*900;
    bs0    = 1250;
    ks0    = 13500;
    k1     = 20;
    k0     = 1;
    S      = 2;
    bsl    = 1385;
    bsnl   = 524;
    kt     = 190000;
    bt     = 14500;
    g      = 9.81;
    
    %% Forces produced by the spring and damper
    fs     = ksl*(x(1) - x(3)) + ksnl*(x(1) - x(3))^3;
    fd     = bsl*(x(2) - x(4)) + bsnl*(x(2) - x(4))^2;
    
    %% Force produced by the tire
    zr     = fcn5(t);
    zr_dot = roadproder(t);
    fst    = kt*(x(3) - zr);
    fdt    = bt*(x(4) - zr_dot);
    if fst + fdt < (ms + mus)*g
        ft = fst + fdt;
    else
        ft = 0;
    end
    
    %% Sliding surface
    sigma   = S*x(1) + x(2);
    
    %% Estimation of disturbance
    d1hat   = x(5) + w*sigma;

    %% State-feedback controller with disturbance canceller
    veq     = - ms0/k0*(S*x(2) + k1*sigma - 1/ms0*(ks0*x(1) + bs0*x(2)));   % state-feedback 
    vn      = - ms0/k0*d1hat;                                               % disturbance canceller
    v       = veq + vn;

    %% Non-ideal actuator
    if     v >= br
        dd  = - m*br;
    elseif v >  bl && v < br
        dd  = - m*v;
    elseif v <= bl
        dd  = - m*bl;
    end
    N       = m*v + dd;
    
    %% Differential equations
    dxdt    = zeros(5,1);
    %  Equation of motion
    dxdt(1) = x(2);
    dxdt(2) = - 1/ms*(fs + fd - N);
    dxdt(3) = x(4);
    dxdt(4) =  1/mus*(fs + fd - ft - N);
    %  Disturbance Observer
    dxdt(5) = - w*(S*x(2) - 1/ms0*(ks0*x(1) + bs0*x(2)) + d1hat + k0/ms0*v);
end

%% Local functions
%  -------------------
%  Road profile, zr(t)
function zr = fcn5(t)
    %   times
    t1 = t - 3.5;
    t2 = t - 6.5;
    t3 = t - 8.5;
    t4 = t - 11.5;
    %   d(t)
    d  = 0.002*sin(2*pi*t) + 0.002*sin(7.5*pi*t);
    %   eq 45
    if (t >= 3.5 && t < 5) 
        zr = -0.0592*t1^3 + 0.1332*t1^2 + d;
    elseif (t >= 5 && t < 6.5)
        zr = 0.0592*t2^3 + 0.1332*t2^2 + d;
    elseif (t >= 8.5 && t < 10)
        zr = 0.0592*t3^3 - 0.1332*t3^2 + d;
    elseif (t >= 10 && t < 11.5)
        zr = -0.0592*t4^3 - 0.1332*t4^2 + d;
    else
        zr = d;
    end
end

%  Derivative of zr
function zr_dot = roadproder(t)
    t1 = t - 3.5;
    t2 = t - 6.5;
    t3 = t - 8.5;
    t4 = t - 11.5;
    %  derivative of d(t)
    d = 2*pi*0.002*cos(2*pi*t) + 7.5*pi*0.002*cos(7.5*pi*t);
    %   eq 45
    if (t >= 3.5 && t < 5) 
        zr_dot = 3*(-0.0592*t1^2) + 2*(0.1332*t1) + d;
    elseif (t >= 5 && t < 6.5)
        zr_dot = 3*(0.0592*t2^2) + 2*(0.1332*t2) + d;
    elseif (t >= 8.5 && t < 10)
        zr_dot = 3*(0.0592*t3^2) - 2*(0.1332*t3)+ d;
    elseif (t >= 10 && t < 11.5)
        zr_dot = 3*(-0.0592*t4^2)- 2*(0.1332*t4)+ d;
    else
        zr_dot = d;
    end
end