필터 지우기
필터 지우기

Matlab simulation code error

조회 수: 1 (최근 30일)
hai
hai 2023년 10월 28일
답변: Sulaymon Eshkabilov 2023년 10월 28일
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');

이 질문에 답변하려면 로그인하십시오.

답변 (1개)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023년 10월 28일
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');

카테고리

Help Center 및 File Exchange에서 Simulink에 대해 자세히 알아보기

태그

제품


릴리스

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by