Out of memory when using ode45(@t,x)

function xdot=r2dof(t,x,ths,spec,Kpid)
xdot=zeros(8,1);
%% set-points
th1s=ths(1);
th2s=ths(2);
%% Robot Specifications
M1=spec(3);
M2=spec(4);
L1=spec(1);
L2=spec(2);
g=9.8;
%% Inertia Matrix
b11=(M1+M2)*L1^2+M2*L2^2+2*M2*L1*L2*cos(x(4));
b12=M2*L2^2+M2*L1*L2*cos(x(4));
b21=M2*L2^2+M2*L1*L2*cos(x(4));
b22=M2*L2^2;
Bq=[b11 b12;b21 b22];
%% C Matrix
c1=-M2*L1*L2*sin(x(4))*(2*x(5)*x(6)+x(6)^2);
c2=-M2*L1*L2*sin(x(4))*x(5)*x(6);
Cq=[c1;c2];
%% Gravity Matrix
g1=-(M1+M2)*g*L1*sin(x(3))-M2*g*L2*sin(x(3)+x(4));
g2=-M2*g*L2*sin(x(3)+x(4));
Gq=[g1;g2];
%% PID Control
% PID parameters for theta 1
Kp1=Kpid(1);
Kd1=Kpid(2);
Ki1=Kpid(3);
% PID parameters for theta 2
Kp2=Kpid(4);
Kd2=Kpid(5);
Ki2=Kpid(6);
%decoupled control input
f1=Kp1*(th1s-x(3))-Kd1*x(5)+Ki1*(x(1));
f2=Kp2*(th2s-x(4))-Kd2*x(6)+Ki2*(x(2));
Fhat=[f1;f2];
F=Bq*Fhat; % actual input to the system
%% System states
xdot(1)=(th1s-x(3)); %dummy state of theta1 integration
xdot(2)=(th2s-x(4)); %dummy state of theta2 integration
xdot(3)=x(5); %theta1-dot
xdot(4)=x(6); %theta2-dot
q2dot=inv(Bq)*(-Cq-Gq+F);
xdot(5)=q2dot(1); %theta1-2dot
xdot(6)=q2dot(2); %theta1-2dot
%control input function output to outside computer program
xdot(7)=F(1);
xdot(8)=F(2);

2) System Solution and Simulation (“r2dof_cntrl.m”)

close all
clear all
clc
%% Initilization
th_int=[-pi/2 pi/2]; %initial positions
ths=[pi/2 -pi/2]; %set-points
x0=[0 0 th_int 0 0 0 0]; %states initial values
Ts=[0 20]; %time span
%% Robot Specifications
L1=1; %link 1
L2=1; %link 2
M1=1; %mass 1
M2=1; %mass 2
spec=[L1 L2 M1 M2];
%% PID Parameters
% PID parameters for theta 1
Kp1=15;
Kd1=7;
Ki1=10;
% PID parameters for theta 2
Kp2=15;
Kd2=10;
Ki2=10;
Kpid=[Kp1 Kd1 Ki1 Kp2 Kd2 Ki2];
%% ODE solving
% opt1=odeset('RelTol',1e-10,'AbsTol',1e-20,'NormControl','off');
[T,X] = ode45(@(t,x) r2dof(t,x,ths,spec,Kpid),Ts,x0);

//////////////////////////////////////////////////////

ERROR

Out of memory. The likely cause is an infinite recursion within the program.

Error in r2dof (line 75)

[T,X] = ode45(@(t,x) r2dof(1, [0 0 th_int 0 0 0 0], [pi/2 -pi/2], [1 1 1 1],[1 1 1 1 1 1]),[0 20],1);

[Torsten]

The command causing the error is not included in the code you provided:

[T,X] = ode45(@(t,x) r2dof(1, [0 0 th_int 0 0 0 0], [pi/2 -pi/2], [1 1 1 1],[1 1 1 1 1 1]),[0 20],1);

And r2dof called like this is independent of t and x which makes no sense.

[mohammed tifr]

first thank you for your answer then

the modified code

[T,X] = ode45(@(t,x) r2dof(t,x, [pi/2 -pi/2], [1 1 1 1],[1 1 1 1 1 1]),[0 20],1);

Errors(see the capture)

[Torsten]

Works for me.

But your call to ode45 is wrong.

[T,X] = ode45(@(t,x)r2dof(t,x,[pi/2 -pi/2],[1 1 1 1],[1 1 1 1 1 1 ]),tspan,x0)

would be correct.

r2dof(1,[0 0 [-pi/2 pi/2] 0 0 0 0 ],[pi/2 -pi/2],[1 1 1 1],[1 1 1 1 1 1 ])

ans = 8×1

3.1416 -3.1416 0 0 -6.6584 6.6584 6.2832 0

function xdot=r2dof(t,x,ths,spec,Kpid)
xdot=zeros(8,1);
%% set-points
th1s=ths(1);
th2s=ths(2);
%% Robot Specifications
M1=spec(3);
M2=spec(4);
L1=spec(1);
L2=spec(2);
g=9.8;
%% Inertia Matrix
b11=(M1+M2)*L1^2+M2*L2^2+2*M2*L1*L2*cos(x(4));
b12=M2*L2^2+M2*L1*L2*cos(x(4));
b21=M2*L2^2+M2*L1*L2*cos(x(4));
b22=M2*L2^2;
Bq=[b11 b12;b21 b22];
%% C Matrix
c1=-M2*L1*L2*sin(x(4))*(2*x(5)*x(6)+x(6)^2);
c2=-M2*L1*L2*sin(x(4))*x(5)*x(6);
Cq=[c1;c2];
%% Gravity Matrix
g1=-(M1+M2)*g*L1*sin(x(3))-M2*g*L2*sin(x(3)+x(4));
g2=-M2*g*L2*sin(x(3)+x(4));
Gq=[g1;g2];
%% PID Control
% PID parameters for theta 1
Kp1=Kpid(1);
Kd1=Kpid(2);
Ki1=Kpid(3);
% PID parameters for theta 2
Kp2=Kpid(4);
Kd2=Kpid(5);
Ki2=Kpid(6);
%decoupled control input
f1=Kp1*(th1s-x(3))-Kd1*x(5)+Ki1*(x(1));
f2=Kp2*(th2s-x(4))-Kd2*x(6)+Ki2*(x(2));
Fhat=[f1;f2];
F=Bq*Fhat; % actual input to the system
%% System states
xdot(1)=(th1s-x(3)); %dummy state of theta1 integration
xdot(2)=(th2s-x(4)); %dummy state of theta2 integration
xdot(3)=x(5); %theta1-dot
xdot(4)=x(6); %theta2-dot
q2dot=inv(Bq)*(-Cq-Gq+F);
xdot(5)=q2dot(1); %theta1-2dot
xdot(6)=q2dot(2); %theta1-2dot
%control input function output to outside computer program
xdot(7)=F(1);
xdot(8)=F(2);
end

Answer 1

[Jan]

I've copied your code without any changes (but adding en end at the bottom) and it is working:

%% Initilization
th_int=[-pi/2 pi/2]; %initial positions
ths=[pi/2 -pi/2]; %set-points
x0=[0 0 th_int 0 0 0 0]; %states initial values
Ts=[0 20]; %time span
%% Robot Specifications
L1=1; %link 1
L2=1; %link 2
M1=1; %mass 1
M2=1; %mass 2
spec=[L1 L2 M1 M2];
%% PID Parameters
% PID parameters for theta 1
Kp1=15;
Kd1=7;
Ki1=10;
% PID parameters for theta 2
Kp2=15;
Kd2=10;
Ki2=10;
Kpid=[Kp1 Kd1 Ki1 Kp2 Kd2 Ki2];
%% ODE solving
% opt1=odeset('RelTol',1e-10,'AbsTol',1e-20,'NormControl','off');
[T,X] = ode45(@(t,x) r2dof(t,x,ths,spec,Kpid),Ts,x0);
plot(T, X)

function xdot=r2dof(t,x,ths,spec,Kpid)
xdot=zeros(8,1);
%% set-points
th1s=ths(1);
th2s=ths(2);
%% Robot Specifications
M1=spec(3);
M2=spec(4);
L1=spec(1);
L2=spec(2);
g=9.8;
%% Inertia Matrix
b11=(M1+M2)*L1^2+M2*L2^2+2*M2*L1*L2*cos(x(4));
b12=M2*L2^2+M2*L1*L2*cos(x(4));
b21=M2*L2^2+M2*L1*L2*cos(x(4));
b22=M2*L2^2;
Bq=[b11 b12;b21 b22];
%% C Matrix
c1=-M2*L1*L2*sin(x(4))*(2*x(5)*x(6)+x(6)^2);
c2=-M2*L1*L2*sin(x(4))*x(5)*x(6);
Cq=[c1;c2];
%% Gravity Matrix
g1=-(M1+M2)*g*L1*sin(x(3))-M2*g*L2*sin(x(3)+x(4));
g2=-M2*g*L2*sin(x(3)+x(4));
Gq=[g1;g2];
%% PID Control
% PID parameters for theta 1
Kp1=Kpid(1);
Kd1=Kpid(2);
Ki1=Kpid(3);
% PID parameters for theta 2
Kp2=Kpid(4);
Kd2=Kpid(5);
Ki2=Kpid(6);
%decoupled control input
f1=Kp1*(th1s-x(3))-Kd1*x(5)+Ki1*(x(1));
f2=Kp2*(th2s-x(4))-Kd2*x(6)+Ki2*(x(2));
Fhat=[f1;f2];
F=Bq*Fhat; % actual input to the system
%% System states
xdot(1)=(th1s-x(3)); %dummy state of theta1 integration
xdot(2)=(th2s-x(4)); %dummy state of theta2 integration
xdot(3)=x(5); %theta1-dot
xdot(4)=x(6); %theta2-dot
q2dot=inv(Bq)*(-Cq-Gq+F);
xdot(5)=q2dot(1); %theta1-2dot
xdot(6)=q2dot(2); %theta1-2dot
%control input function output to outside computer program
xdot(7)=F(1);
xdot(8)=F(2);
end

So if you get errors, you do not run the posted code.

If you still have problem, post the code you run.

Answer 2

[Sam Chak]

Hi @mohammed tifr

The desired control torques cannot be integrated directly in and , because the actual torques are delivered by the motors. This, some minor modifications are made to assume and include the motor dynamics in and , as well as joint dynamics and .

However, I didn't check if the kinematics and are true for large angles. Swinging ° to ° are relatively large motion. Logically, there should be some constraints (or maybe singularity), like the elbow and shoulder.

%% Initilization
th_int = [-pi/2  pi/2];         % initial positions
ths    = [ pi/2 -pi/2];         % setpoints
x0     = [0 0 th_int 0 0 0 0];  % states initial values
Ts     = linspace(0, 10, 1001); % time span

%% Robot Specifications
L1   = 1;       % link 1
L2   = 1;       % link 2
M1   = 1;       % mass 1
M2   = 1;       % mass 2
spec = [L1 L2 M1 M2];

%% PID Parameters
% PID parameters for theta 1
Kp1 = 15;
Kd1 =  7;
Ki1 = 10;
% PID parameters for theta 2
Kp2 = 15;
Kd2 = 10;
Ki2 = 10;

Kpid = [Kp1 Kd1 Ki1 Kp2 Kd2 Ki2];

%% ODE solving
% opt1=odeset('RelTol',1e-10,'AbsTol',1e-20,'NormControl','off');
[t, x] = ode45(@(t,x) r2dof(t, x, ths, spec, Kpid), Ts, x0);
plot(t, x(:,3:4)*180/pi), grid on, ylabel('Angle, [\circ]')
legend('\theta_1', '\theta_2', 'fontsize', 16)
yticks(-150:60:150)

plot(t, x(:,7:8)), grid on, ylabel('Torque, [Nm]')
legend('\tau_1', '\tau_2', 'fontsize', 16)

function xdot = r2dof(t, x, ths, spec, Kpid)
    xdot = zeros(8, 1);
    
    %% set-points
    th1s = ths(1);
    th2s = ths(2);
    
    %% Robot Specifications
    M1 = spec(3);
    M2 = spec(4);
    L1 = spec(1);
    L2 = spec(2);
    g  = 9.8;
    
    %% Inertia Matrix
    b11 = (M1 + M2)*L1^2 + M2*L2^2 + 2*M2*L1*L2*cos(x(4));
    b12 = M2*L2^2 + M2*L1*L2*cos(x(4));
    b21 = M2*L2^2 + M2*L1*L2*cos(x(4));
    b22 = M2*L2^2;
    Bq  = [b11 b12; b21 b22];
    
    %% C Matrix
    c1  = - M2*L1*L2*sin(x(4))*(2*x(5)*x(6) + x(6)^2);
    c2  = - M2*L1*L2*sin(x(4))*x(5)*x(6);
    Cq  = [c1; c2];
    
    %% Gravity Matrix
    g1  = - (M1 + M2)*g*L1*sin(x(3)) - M2*g*L2*sin(x(3) + x(4));
    g2  = - M2*g*L2*sin(x(3) + x(4));
    Gq  = [g1; g2];
    
    %% PID Control
    % PID parameters for theta 1
    Kp1 = Kpid(1);
    Kd1 = Kpid(2);
    Ki1 = Kpid(3);
    % PID parameters for theta 2
    Kp2 = Kpid(4);
    Kd2 = Kpid(5);
    Ki2 = Kpid(6);
    
    % decoupled control input
    f1   = - Kp1*(x(3) - th1s) - Kd1*x(5) - Ki1*(x(1));
    f2   = - Kp2*(x(4) - th2s) - Kd2*x(6) - Ki2*(x(2));
    Fhat = [f1; f2];
    F    = Bq*Fhat;         % desired torques computed by formulas to the system, x7 & x8 are the actual torques delivered by the motors to the joints
    
    %% System states
    xdot(1) = x(3) - th1s;  % dummy state of theta1 integration
    xdot(2) = x(4) - th2s;  % dummy state of theta2 integration
    xdot(3) = x(5);         % theta1-dot
    xdot(4) = x(6);         % theta2-dot
    
    q2dot   = inv(Bq)*(- Cq - Gq + [x(7); x(8)]);   % actuated torques are delivered by the motors
    xdot(5) = q2dot(1);     % theta1-2dot
    xdot(6) = q2dot(2);     % theta1-2dot
    
    % control input function output to outside computer program
    % (assume the motor is fast enough to be approximated by a 1st-order function)
    xdot(7) = - 1e2*(x(7) - F(1));  % can modify the motor constant
    xdot(8) = - 1e2*(x(8) - F(2));
end