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

1 visualizzazione (ultimi 30 giorni)
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

Risposta accettata

Sam Chak
Sam Chak il 19 Mag 2024
Modificato: Sam Chak il 19 Mag 2024
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 Commenti
Sam Chak
Sam Chak il 19 Mag 2024
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 il 19 Mag 2024
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

Accedi per commentare.

Più risposte (0)

Prodotti


Release

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by