Borrar filtros
Borrar filtros

Solve an ODE with the parameters defined in the function changing in a for loop in the main script

23 visualizaciones (últimos 30 días)
I have to cycle trough different values of stifness kp in the main script with a for loop and solve an Ode. How can I modify the kp value in the function defined used as imput to an ode45 in the main script?
Thanks in advance
function xdot = Time_domain_system(t,x)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
kp = 50;
%equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 - m*g*L/2;
xdot(1) = x(2);
xdot(2) = kp*theta_ref/mx -(cx/mx)*x(2) - ((kx + kp)/mx)*x(1);
xdot = xdot';
end

Respuesta aceptada

Sam Chak
Sam Chak el 19 de Mayo de 2024 a las 7:34
Editada: Sam Chak el 19 de Mayo de 2024 a las 9:06
Let me know if the looping approach I provided is helpful. By the way, your original proportional-only control law for xdot(2) could not drive the angular position to the desired θ reference of 1. So, I made a small modification to achieve that goal wonderfully. You can compare the code to study what changes were made.
kp = [1 2 4 8 16];
for j = 1:length(kp)
[t, x] = ode45(@(t, x) Time_domain_system(t, x, kp(j)), [0 20], [0; 0]);
plot(t, x(:,1)), hold on
end
hold off, grid, xlabel t
function xdot = Time_domain_system(t, x, kp)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
% kp = 50;
% equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 - m*g*L/2;
xdot(1) = x(2);
% xdot(2) = kp*theta_ref/mx - (cx/mx)*x(2) - ((kx + kp )/mx)*x(1); % original
xdot(2) = kp*theta_ref - (cx/mx)*x(2) - ((kx + (mx*kp - kx))/mx)*x(1); % my proposal
xdot = xdot';
end
  3 comentarios
Sam Chak
Sam Chak el 19 de Mayo de 2024 a las 14:24
You are welcome, @Paolo Trenta. If you find both the for-loop and proportional control law helpful, please consider clicking 'Accept' ✔️ on the answer.
Sam Chak
Sam Chak el 19 de Mayo de 2024 a las 14:36
I found that if you add a single term to your original equation in xdot(2), then the angular position will be regulated to 1. Of course, the proportional gain also needs to be scaled accordingly.
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
mx = J+m*L^2/4;
kp = mx*[1 2 4 8 16];
for j = 1:length(kp)
[t, x] = ode45(@(t, x) Time_domain_system(t, x, kp(j)), [0 20], [0; 0]);
plot(t, x(:,1)), hold on
end
hold off, grid, xlabel t
function xdot = Time_domain_system(t, x, kp)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
% kp = 50;
% equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 - m*g*L/2;
xdot(1) = x(2);
% xdot(2) = kp*theta_ref/mx - (cx/mx)*x(2) - ((kx + kp )/mx)*x(1); % original
xdot(2) = kp*theta_ref/mx - (cx/mx)*x(2) - ((kx + (kp - kx))/mx)*x(1); % 2nd proposal
xdot = xdot';
end

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Assembly en Help Center y File Exchange.

Productos


Versión

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by