Clear Filters
Clear Filters

Matlab simulation code error

1 view (last 30 days)
hai
hai on 28 Oct 2023
Answered: Sulaymon Eshkabilov on 28 Oct 2023
I have a mechanism to simulate but it get an errors code
How to solve it?
I try to make a simulation from book exercise(Mechanisms and Robots Analysis with MATLAB®) but it get an error, please help?
clear all; clc; close all
% Input data
AB=0.22; %(m)
AC=0.08; %(m)
CD=0.20; %(m)
DE=0.60; %(m)
M = moviein(12);
incr = 0;
VB23=zeros(1,360);
VB3=zeros(1,360);
VBD=zeros(1,360);
VBE=zeros(1,360);
AB1=zeros(1,360);
AB23=zeros(1,360);
AB3=zeros(1,360);
ABD=zeros(1,360);
ABE=zeros(1,360);
TANPHI=zeros(1,360);
i=1;
for phi =0 : pi/90:6*pi
%phi = 4*pi/3 ; %(rad)
fprintf('\nphi =%g (rad)\n', phi)
n = 500*2 * pi/60;
omega1 = [0 0 n ]; % (rad/s)
fprintf('omega1 = omega2 = [ %g, %g, %g ] (rad/s)\n', omega1)
% Position of joint A
fprintf('\nPoint A \n')
xA=0;yA=0;rA = [xA yA 0];
xC=0.08; yC = 0; rC = [xC yC 0];
fprintf('rA = [ %g, %g, %g ] (m)\n', rA)
vA=[0 0 0];%(m/s) % velocity of A (fixed)
fprintf('vA = [ %g, %g, %g ] (m/s)\n', vA)
fprintf('|vA|= %g (m/s)\n',norm(vA))
% Position of joint B1
fprintf('\nPoint B \n')
xB=AB*cos(phi); yB = AB*sin(phi); rB = [xB yB 0];
fprintf('rB = [ %g, %g, %g ] (m)\n', rB)
vB1 = vA + cross(omega1,rB); % velocity of B1
fprintf('vB1 = vB2 = [ %g, %g, %g ] (m/s)\n', vB1)
fprintf('|vB1| = |vB2| = %g (m/s)\n', norm(vB1))
vB2x=vB1(1);
vB2y=vB1(2);
omega11=norm(omega1);
aB1 = [-(omega11)^2*xB -(omega11)^2*yB 0];
fprintf('aB1 = aB2 = [ %g, %g, %g ] (m/s2)\n', aB1)
%%%%%%
% tan phi3
tanphi3 = (yB-yC)/(xB-xC);
phi3=atan(tanphi3);
fprintf('phi 3 = %g (rad/s)\n', phi3)
cosphi3= cos(phi3);
fprintf('cosphi3 = %g (rad/s)\n', cosphi3)
sinphi3= sin(phi3);
fprintf('sinphi3 = %g (rad/s)\n', sinphi3)
%%%
eqnB3x = 'vB2x + vB23x = vB3x ';
eqnB23x = '-vB3x*((xC-xB)/(yC-yB))=vB2y + vB23x*((yB-yC)/(xB-xC)) ';
solB = solve(eqnB3x, eqnB23x, 'vB3x','vB23x');
vB3x= eval(solB.vB3x);
vB23x= eval(solB.vB23x);
eqnB3y = 'vB2y + vB23y = vB3y ';
eqnB23y = '-vB3y*((yC-yB)/(xC-xB))=vB2y + vB23y*((xB-xC)/(yB-yC)) ';
solB = solve(eqnB3y, eqnB23y, 'vB3y','vB23y');
vB3y= eval(solB.vB3y);
vB23y= eval(solB.vB23y);
vB3 = [vB3x vB3y 0];
vB23 = [vB23x vB23y 0];
omega3 = [0 0 vB3y/(xB-xC)];
fprintf('omega3 = [ %g, %g, %g ] (rad/s)\n', omega3)
fprintf('vB3 = [ %g, %g, %g ] \n', vB3)
fprintf('|vB3| = %g (m/s)\n', norm(vB3))
fprintf('vB23 = [ %g, %g, %g ] \n', vB23)
eqnalpha3 = '-alpha3*(yC-yB) = aB1(1) + aB23*cos(phi3) -vB23(1)*2*sin(90*pi/360-phi3) ';
eqnaB23 = 'alpha3*(xC-xB) = aB1(2) + aB23*sin(phi3) +vB23(2)*2*cos(90*pi/360-phi3) ';
solalpha = solve(eqnalpha3, eqnaB23,'alpha3',' aB23');%giai pt
alpha3pos = eval(solalpha.alpha3);
aB23pos = eval(solalpha.aB23);
fprintf('alpha3 = %g \n', alpha3pos )
fprintf('aB32pos = %g \n', aB23pos )
aB23p = [aB23pos*sinphi3 aB23pos*cosphi3 0];
fprintf('aB32 = [ %g, %g, %g ] \n', aB23p)
aB3 = alpha3pos*(rB-rC);
fprintf('aB3 = [ %g, %g, %g ] \n', aB3)
% Position of joint C
fprintf('\nPoint C \n')
xC=0.08; yC = 0; rC = [xC yC 0];
fprintf('rC = [ %g, %g, %g ] (m)\n', rC)
vC = [0 0 0];
fprintf('vC = 0 (m/s)\n')
fprintf('|vC|= 0\n')
% Position of joint D
fprintf('\nPoint D \n')
eqnD1 = '(xDsol - xC)^2 + (yDsol - yC)^2 = CD^2 ';
eqnD2 = '(yB/(xB-xC))*(yDsol/(xDsol - xC)) = -1';
solD = solve(eqnD1, eqnD2, 'xDsol',' yDsol');%giai pt
xDpositions = eval(solD.xDsol);
yDpositions = eval(solD.yDsol);
xD1= xDpositions(1);
xD2= xDpositions(2);
yD1= yDpositions(1);
yD2= yDpositions(2);
xD = xD1;
yD=yD1;
rD = [xD yD 0];
fprintf('rD = [ %g, %g, %g ] (m)\n', rD)
vD = vC + cross(omega3, rD-rC);
fprintf('vD = vD5 = [ %g, %g, %g ] (m/s)\n', vD)
fprintf('|vD| = %g (m/s)\n', norm(vD))
vDx=vD(1);
vDy=vD(2);
aD= alpha3pos*(rD-rC) - norm(omega3)^2*(rD-rC);
fprintf('aD = [ %g, %g, %g ] (m/s2)\n', aD)
% Position of joint E
eqnE = '(xEsol - xD)^2 + (yD)^2 = DE^2 ';
solE = solve(eqnE, 'xEsol');
xEpositions = eval(solE);
xE1= xEpositions(1);
xE2= xEpositions(2);
if xE1> xC
xE = xE1;
else
xE = xE2;
end
yE=0;
rE = [xE yE 0];
fprintf('\nPoint E \n')
fprintf('rE = [ %g, %g, %g ] (m)\n', rE)
eqnE1 = 'vDx - omega4*(yE-yD) = vEx ';
eqnE2 = 'vDy + omega4*(xE-xD) = 0';
solD = solve(eqnE1, eqnE2, 'vEx',' omega4');%giai pt
vEx = eval(solD.vEx);
omega4x = eval(solD.omega4);
vE = [vEx 0 0];
omega4 = [0 0 omega4x];
fprintf('vE = [ %g, %g, %g ] (m/s)\n', vE)
fprintf('|vE| = %g (m/s)\n', norm(vE))
fprintf('omega4 = [ %g, %g, %g ] \n', omega4)
aDx=aD(1);
aDy=aD(2);
eqnaEx = 'aEx = aDx + alpha4*(yE-yD) - (norm(omaga4))^2*( yE-yD)';
eqnalpha4 = 'aDy+alpha4*(xE-xD) - (norm(omega4))^2*(xE-xD)=0';
solaE = solve(eqnaEx, eqnalpha4, 'aEx',' alpha4');%giai pt
aExpos = eval(solaE.aEx);
alpha4pos = eval(solaE.alpha4);
fprintf('alpha4 = %g \n', alpha4pos )
fprintf('aE = [%g 0 0] \n', aExpos )
% Graphic of the mechanism
figure(1)
plot([xA,xB],[yA,yB],'k-*',...
[xB,xC],[yB,yC],'b-*',...
[xC,xD],[yC,yD],'g-*',...
[xD,xE],[yD,yE],'r-*','LineWidth',1.5)
% adds major grid lines to the current axes
%grid on,...
xlabel('x (m)'), ylabel('y (m)'),...
title('positions for \phi = 120 (deg)'),...
text(xA,yA,'\leftarrow A = ground',...
'HorizontalAlignment','right'),...
text(xB,yB,'--B'),...
text(xC,yC,'--C'),...
text(xD,yD,'--D'),...
text(xE,yE,'--E'),grid;
%axis([-0.3 0.8 -0.3 0.45])
xlim([-0.3 1.3]);
ylim([-0.3 0.9]);
% end of program
%incr = incr + 1;
%M( : , incr) = getframe; % record the movie
VB23(i)=norm(vB23);
VB3(i)=norm(vB3);
VD(i)=norm(vD);
VE(i)=norm(vE);
AB1(i)=norm(aB1);
AB23(i)=norm(aB23p);
AB3(i)=norm(aB3);
AD(i)=norm(aD);
AE(i)=norm(aExpos);
i=i+1;
end % end for
figure(4)
subplot ( 2, 2,1)
plot(VB23,'r');
title('|vB23|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,2)
plot(VB3,'g');
title('|vB3|')
subplot ( 2, 2,3)
plot(VD,'k');
title('|vD|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,4)
plot(VE,'b');
legend('y = sin(x)','y = cos(x)');
title('|vE|')
xlabel('0 < phi < 6*pi') % x-axis label
figure(5)
subplot ( 2, 2,1)
plot(AB23);
title('|aB23|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,2)
plot(AB3);
title('|aB3|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,3)
plot(AD);
title('|aD|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,4)
plot(AE);
title('|aE|')
xlabel('0 < phi < 6*pi') % x-axis label
%movie2avi(M,'exam7_9_velocity.avi');
Incorrect number or types of inputs or outputs for function 'solve'.
Error in me (line 60)
solB = solve(eqnB3x, eqnB23x, 'vB3x','vB23x');

Sign in to answer this question.

Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 28 Oct 2023
Here is the corrected code.
Note the errors with:
(1) symbolic math expression and solving the equations: syms, solve()
(2) singularity problem with the initial angle at phi=0 degree (0.01*pi taken to avoid this issue). To speed up simulation, take a larger step for phi, e.g.: phi= =0.01*pi:pi/10:6*pi
(3) Suggestion: instead of displaying all values with fprintf(), try to store them.
Have fun!
clear all; clc; close all
% Input data
AB=0.22; %(m)
AC=0.08; %(m)
CD=0.20; %(m)
DE=0.60; %(m)
M = moviein(12);
incr = 0;
VB23=zeros(1,360);
VB3=zeros(1,360);
VBD=zeros(1,360);
VBE=zeros(1,360);
AB1=zeros(1,360);
AB23=zeros(1,360);
AB3=zeros(1,360);
ABD=zeros(1,360);
ABE=zeros(1,360);
TANPHI=zeros(1,360);
i=1;
for phi =0.01*pi:pi/90:6*pi
%phi = 4*pi/3 ; %(rad)
fprintf('\nphi =%g (rad)\n', phi)
n = 500*2 * pi/60;
omega1 = [0 0 n ]; % (rad/s)
fprintf('omega1 = omega2 = [ %g, %g, %g ] (rad/s)\n', omega1)
% Position of joint A
fprintf('\nPoint A \n')
xA=0;yA=0;rA = [xA yA 0];
xC=0.08; yC = 0; rC = [xC yC 0];
fprintf('rA = [ %g, %g, %g ] (m)\n', rA)
vA=[0 0 0];%(m/s) % velocity of A (fixed)
fprintf('vA = [ %g, %g, %g ] (m/s)\n', vA)
fprintf('|vA|= %g (m/s)\n',norm(vA))
% Position of joint B1
fprintf('\nPoint B \n')
xB=AB*cos(phi); yB = AB*sin(phi); rB = [xB yB 0];
fprintf('rB = [ %g, %g, %g ] (m)\n', rB)
vB1 = vA + cross(omega1,rB); % velocity of B1
fprintf('vB1 = vB2 = [ %g, %g, %g ] (m/s)\n', vB1)
fprintf('|vB1| = |vB2| = %g (m/s)\n', norm(vB1))
vB2x=vB1(1);
vB2y=vB1(2);
omega11=norm(omega1);
aB1 = [-(omega11)^2*xB -(omega11)^2*yB 0];
fprintf('aB1 = aB2 = [ %g, %g, %g ] (m/s2)\n', aB1)
%%%%%%
% tan phi3
tanphi3 = (yB-yC)/(xB-xC);
phi3=atan(tanphi3);
fprintf('phi 3 = %g (rad/s)\n', phi3)
cosphi3= cos(phi3);
fprintf('cosphi3 = %g (rad/s)\n', cosphi3)
sinphi3= sin(phi3);
fprintf('sinphi3 = %g (rad/s)\n', sinphi3)
%%%
syms vB23x vB3x
eqnB3x = vB2x + vB23x == vB3x ;
eqnB23x = -vB3x*((xC-xB)/(yC-yB))==vB2y + vB23x*((yB-yC)/(xB-xC));
solB = solve(eqnB3x, eqnB23x, vB3x,vB23x);
vB3x= eval(solB.vB3x);
vB23x= eval(solB.vB23x);
syms vB3y vB23y
eqnB3y = vB2y + vB23y == vB3y ;
eqnB23y = -vB3y*((yC-yB)/(xC-xB))==vB2y + vB23y*((xB-xC)/(yB-yC));
solB = solve(eqnB3y, eqnB23y, vB3y,vB23y);
vB3y= eval(solB.vB3y);
vB23y= eval(solB.vB23y);
vB3 = [vB3x vB3y 0];
vB23 = [vB23x vB23y 0];
omega3 = [0 0 vB3y/(xB-xC)];
fprintf('omega3 = [ %g, %g, %g ] (rad/s)\n', omega3)
fprintf('vB3 = [ %g, %g, %g ] \n', vB3)
fprintf('|vB3| = %g (m/s)\n', norm(vB3))
fprintf('vB23 = [ %g, %g, %g ] \n', vB23)
syms alpha3 aB23
eqnalpha3 = -alpha3*(yC-yB) == aB1(1) + aB23*cos(phi3) -vB23(1)*2*sin(90*pi/360-phi3);
eqnaB23 = alpha3*(xC-xB) == aB1(2) + aB23*sin(phi3) +vB23(2)*2*cos(90*pi/360-phi3);
solalpha = solve(eqnalpha3, eqnaB23, alpha3, aB23);%giai pt
alpha3pos = eval(solalpha.alpha3);
aB23pos = eval(solalpha.aB23);
fprintf('alpha3 = %g \n', alpha3pos )
fprintf('aB32pos = %g \n', aB23pos )
aB23p = [aB23pos*sinphi3 aB23pos*cosphi3 0];
fprintf('aB32 = [ %g, %g, %g ] \n', aB23p)
aB3 = alpha3pos*(rB-rC);
fprintf('aB3 = [ %g, %g, %g ] \n', aB3)
% Position of joint C
fprintf('\nPoint C \n')
xC=0.08; yC = 0; rC = [xC yC 0];
fprintf('rC = [ %g, %g, %g ] (m)\n', rC)
vC = [0 0 0];
fprintf('vC = 0 (m/s)\n')
fprintf('|vC|= 0\n')
% Position of joint D
fprintf('\nPoint D \n')
syms xDsol yDsol
eqnD1 = (xDsol - xC)^2 + (yDsol - yC)^2 == CD^2;
eqnD2 = (yB/(xB-xC))*(yDsol/(xDsol - xC)) == -1;
solD = solve(eqnD1, eqnD2, xDsol, yDsol);%giai pt
xDpositions = eval(solD.xDsol);
yDpositions = eval(solD.yDsol);
xD1= xDpositions(1);
xD2= xDpositions(2);
yD1= yDpositions(1);
yD2= yDpositions(2);
xD = xD1;
yD=yD1;
rD = [xD yD 0];
fprintf('rD = [ %g, %g, %g ] (m)\n', rD)
vD = vC + cross(omega3, rD-rC);
fprintf('vD = vD5 = [ %g, %g, %g ] (m/s)\n', vD)
fprintf('|vD| = %g (m/s)\n', norm(vD))
vDx=vD(1);
vDy=vD(2);
aD= alpha3pos*(rD-rC) - norm(omega3)^2*(rD-rC);
fprintf('aD = [ %g, %g, %g ] (m/s2)\n', aD)
% Position of joint E
syms xEsol
eqnE = (xEsol - xD)^2 + (yD)^2 == DE^2 ;
solE = solve(eqnE, xEsol);
xEpositions = eval(solE);
xE1= xEpositions(1);
xE2= xEpositions(2);
if xE1> xC
xE = xE1;
else
xE = xE2;
end
yE=0;
rE = [xE yE 0];
fprintf('\nPoint E \n')
fprintf('rE = [ %g, %g, %g ] (m)\n', rE)
syms vEx omega4
eqnE1 = vDx - omega4*(yE-yD) == vEx ;
eqnE2 = vDy + omega4*(xE-xD) == 0;
solD = solve(eqnE1, eqnE2, vEx, omega4);%giai pt
vEx = eval(solD.vEx);
omega4x = eval(solD.omega4);
vE = [vEx 0 0];
omega4 = [0 0 omega4x];
fprintf('vE = [ %g, %g, %g ] (m/s)\n', vE)
fprintf('|vE| = %g (m/s)\n', norm(vE))
fprintf('omega4 = [ %g, %g, %g ] \n', omega4)
aDx=aD(1);
aDy=aD(2);
syms aEx alpha4
eqnaEx = aEx == aDx + alpha4*(yE-yD) - (norm(omega4))^2*( yE-yD);
eqnalpha4 = aDy+alpha4*(xE-xD) - (norm(omega4))^2*(xE-xD)==0;
solaE = solve(eqnaEx, eqnalpha4, aEx, alpha4);%giai pt
aExpos = eval(solaE.aEx);
alpha4pos = eval(solaE.alpha4);
fprintf('alpha4 = %g \n', alpha4pos )
fprintf('aE = [%g 0 0] \n', aExpos )
% Graphic of the mechanism
figure(1)
plot([xA,xB],[yA,yB],'k-*',...
[xB,xC],[yB,yC],'b-*',...
[xC,xD],[yC,yD],'g-*',...
[xD,xE],[yD,yE],'r-*','LineWidth',1.5)
% adds major grid lines to the current axes
%grid on,...
xlabel('x (m)'), ylabel('y (m)'),...
title('positions for \phi = 120 (deg)'),...
text(xA,yA,'\leftarrow A = ground',...
'HorizontalAlignment','right'),...
text(xB,yB,'--B'),...
text(xC,yC,'--C'),...
text(xD,yD,'--D'),...
text(xE,yE,'--E'),grid;
%axis([-0.3 0.8 -0.3 0.45])
xlim([-0.3 1.3]);
ylim([-0.3 0.9]);
% end of program
%incr = incr + 1;
%M( : , incr) = getframe; % record the movie
VB23(i)=norm(vB23);
VB3(i)=norm(vB3);
VD(i)=norm(vD);
VE(i)=norm(vE);
AB1(i)=norm(aB1);
AB23(i)=norm(aB23p);
AB3(i)=norm(aB3);
AD(i)=norm(aD);
AE(i)=norm(aExpos);
i=i+1;
end % end for
figure(4)
subplot ( 2, 2,1)
plot(VB23,'r');
title('|vB23|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,2)
plot(VB3,'g');
title('|vB3|')
subplot ( 2, 2,3)
plot(VD,'k');
title('|vD|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,4)
plot(VE,'b');
legend('y = sin(x)','y = cos(x)');
title('|vE|')
xlabel('0 < phi < 6*pi') % x-axis label
figure(5)
subplot ( 2, 2,1)
plot(AB23);
title('|aB23|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,2)
plot(AB3);
title('|aB3|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,3)
plot(AD);
title('|aD|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,4)
plot(AE);
title('|aE|')
xlabel('0 < phi < 6*pi') % x-axis label
%movie2avi(M,'exam7_9_velocity.avi');

Categories

Find more on Simulink in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!