How to code algebric loop in MATLAB script file?

1 次查看（过去 30 天）
Faezeh Yousefpour
Faezeh Yousefpour 2023-9-20
回答： Sam Chak 2023-9-30
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
  2 个评论
Dyuman Joshi
Dyuman Joshi 2023-9-20
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
Faezeh Yousefpour 2023-9-27

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])

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Sam Chak
Sam Chak 2023-9-30
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by