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

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

Answer 1

[Sam Chak]

Hi @Paolo Trenta

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

[Paolo Trenta]

Thank you Sam!

[Sam Chak]

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]

Hi @Paolo Trenta

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