Out of memory when using ode45(@t,x)
4 views (last 30 days)
Show older comments
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);
3 Comments
Torsten
on 10 Dec 2022
Edited: Torsten
on 10 Dec 2022
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 ])
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
Answers (2)
Jan
on 10 Dec 2022
Edited: Jan
on 10 Dec 2022
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.
0 Comments
Sam Chak
on 11 Dec 2022
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
0 Comments
See Also
Categories
Find more on PID Controller Tuning in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!