Category: Matlab
Category Archives: Matlab
When I run this code in Matlab, it pops up “Function or variable ‘six_bar’ is not recognized”. How to solve the error?
clear;
i1=0.1;
i3=0.59;
i4=0.159;
i6=0.36;
i61=0.578;
omega1=0.87;
alpha=0;
hd=pi/180;
du=180/pi;
%调用子函数,six_bar返回牛头刨床的位移,角速度,角加速度
for n1 = 1:459
theta1(n1)=-2*pi+5.8119+(n1-1)*hd;
ii=[i1,i3,i4,i6,i61];
[theta,omega,alpha]=six_bar(theta1(n1),omega1,ii);
s3(n1)=theta(1);
theta3(n1)=theta(2);
theta4(n1)=theta(3);
sE(n1)=theta(4);
v2(n1)=omega(1);
omega3(n1)=omega(2);
omega4(n1)=omega(3);
vE(n1)=omega(4);
a2(n1)=alpha(1);
alpha3(n1)=alpha(2);
alpha4(n1)=alpha(3);
aE(n1)=alpha(4);
end
%3.位移,角速度，角加速度和牛头刨床图形输出
figure(1);
n1=1:459;
t=(n1-1)*2*pi/360;
subplot(2,2,1);%绘角位移及位移线图
plot(t,theta3*du,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角位移及位移线图’)
axis auto;
[haxes,~,~]=plotyy(t,theta4*du,t,sE); %
grid on;
hold on;
xlabel(‘时间/s’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角位移/circ’);
axes(haxes(2));
ylabel(‘位移/m’);
hold on;
grid on;
text(1.15,-0.65,’theta_3′);
text(3.4,0.27,’theta_4′);
text(2.25,-0.15,’s_E’);
subplot(2,2,2); %绘角速度及速度线图
plot(t,omega3,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角速度及速度线图’)
axis auto;
[haxes,~,~]=plotyy(t,omega4,t,vE); %
grid on;
hold on;
xlabel(‘时间/s’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角速度/radcdots∧{-1}’);
axes(haxes(2));
ylabel(‘速度/mcdots∧{-1}’);
hold on;
grid on;
text(3.1,0.35,’omega_3′);
text(2.1,0.1,’omega_4′);
text(5.5,0.45,’v_E’);
subplot(2,2,3); %绘角加速度及加速度线图
plot(t,alpha3,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角加速度及加速度线图’)
axis auto;
[haxes,~,~]=plotyy(t,alpha4,t,aE); %
grid on;
hold on;
xlabel(‘时间’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角加速度/radcdots∧{-2}’);
axes(haxes(2));
ylabel(‘加速度/mcdots∧{-2}’);
hold on;
grid on;
text(1.5,0.3,’alpha_3′);
text(3.5,0.51,’alpha_4′);
text(1.5,-0.11,’a_E’);
subplot(2,2,4); %牛头刨床机构
n1=20;
x(1)=0;y(1)=0;
x(2)=(s3(n1)*1000-50)*cos(theta3(n1));
y(2)=(s3(n1)*1000-50)*sin(theta3(n1));
x(3)=0;y(3)=i6*1000;
x(4)=i1*1000*cos(theta1(n1));
y(4)=s3(n1)*1000*sin(theta3(n1));
x(5)=(s3(n1)*1000+50)*cos(theta3(n1));
y(5)=(s3(n1)*1000+50)*sin(theta3(n1));
x(6)=i3*1000*cos(theta3(n1));
y(6)=i3*1000*sin(theta3(n1));
x(7)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1));
y(7)=i3*1000*sin(theta3(n1))+i4*1000*sin(theta4(n1));
x(8)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))-900;
y(8)=i61*1000;
x(9)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))+600;
y(9)=i61*1000;
x(10)=(s3(n1)*1000-50)*cos(theta3(n1));
y(10)=(s3(n1)*1000-50)*sin(theta3(n1));
x(11)=x(10)+25*cos(pi/2-theta3(n1));
y(11)=y(10)-25*sin(pi/2-theta3(n1));
x(12)=x(11)+100*cos(theta3(n1));
y(12)=y(11)+100*sin(theta3(n1));
x(13)=x(12)-50*cos(pi/2-theta3(n1));
y(13)=y(12)+50*sin(pi/2-theta3(n1));
x(14)=x(10)-25*cos(pi/2-theta3(n1));
y(14)=y(10)+25*sin(pi/2-theta3(n1));
x(15)=x(10);
y(15)=y(10);
x(16)=0;
y(16)=0;
x(17)=0;
y(17)=i6*1000;
k=1:2;
plot(x(k),y(k));
hold on;
k=3:4;
plot(x(k),y(k));
hold on;
k=5:9;
plot(x(k),y(k));
hold on;
k=10:15;
plot(x(k),y(k));
hold on;
k=16:17;
plot(x(k),y(k));
hold on;
grid on;
axis ([-500 600 0 650]);
title(‘牛头刨床’);
grid on;
xlabel(‘mm’);
ylabel(‘mm’);
plot(x(1),y(1),’o’);
plot(x(3),y(3),’o’);
plot(x(4),y(4),’o’);
plot(x(6),y(6),’o’);
plot(x(7),y(7),’o’);
hold on;
grid on;
xlabel(‘mm’);
ylabel(‘mm’);
axis([-400 600 0 650]);
function [theta, omega, alpha]=six_bar(theta1,omega1,ii)
i1=ii(1);
i3=ii(2);
i4=ii(3);
i6=ii(4);
i61=ii(5);
%1.计算角位移和线位移
s3 =sqrt((i1*cos(theta1))*(i1*cos(theta1))+(i6+i1*sin(theta1))*(i6+i1*sin(theta1)));
%s3表示滑块2相对于CD杆的位移
theta3 =acos((i1*cos(theta1 ))/s3 );
%theta3表示杆3转过角度
theta4 =pi-asin((i61-i3*sin(theta3 ))/i4);
%theta4表示杆4转过角度
sE =i3*cos(theta3 )+i4*cos(theta4 );
%sE表示杆5的位移
theta(1)=s3;
theta(2)=theta3;
theta(3)=theta4;
theta(4)=sE;
%2.计算角速度和线速度
A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
-cos(theta3 ),s3 *sin(theta3 ),0,0;
0,i3*sin(theta3 ),i4*sin(theta4 ),1;
0,i3*cos(theta3 ),i4*cos(theta4 ),0];
B=[i1*cos(theta1 );i1*sin(theta1 );0;0];
omega=A(omega1*B);
v2 =omega(1);
%滑块2的速度
omega3 =omega(2);
%构件3的角速度
omega4 =omega(3);
%构件4的角速度
vE =omega(4);
%构件5的速度
%3.计算角加速度和加速度
A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
%从动件位置参数矩阵
cos(theta3 ),-s3 *sin(theta3 ),0,0;
0,i3*sin(theta3 ),i4*sin(theta4 ),1;
0,i3*cos(theta3 ),i4*cos(theta4 ),0];
At=[omega3 *cos(theta3 ),(v2 *cos(theta3 )-s3*omega3 *sin(theta3 )),0,0;
-omega3 *sin(theta3 ),(-v2 *sin(theta3 )-s3 *omega3 *cos(theta3 )),0,0;
0,i3*omega3 *cos(theta3 ),i4*omega4 *cos(theta4 ),0;
0,-i3*omega3 *sin(theta3 ),-i4*omega4 *sin(theta4 ),0];
Bt=[-i1*omega1*sin(theta1 );-i1*omega1*cos(theta1 );0;0];
alpha=A(-At*omega+omega1*Bt);
%机构从动件的加速度列阵
a2 =alpha(1);
%a2表示滑块2的加速度
alpha3 =alpha(2);
%alpha3表示杆3的角加速度
alpha4 =alpha(3);
%alpha4表示杆4的角加速度
aE =alpha(4);
%构件5的加速度
endclear;
i1=0.1;
i3=0.59;
i4=0.159;
i6=0.36;
i61=0.578;
omega1=0.87;
alpha=0;
hd=pi/180;
du=180/pi;
%调用子函数,six_bar返回牛头刨床的位移,角速度,角加速度
for n1 = 1:459
theta1(n1)=-2*pi+5.8119+(n1-1)*hd;
ii=[i1,i3,i4,i6,i61];
[theta,omega,alpha]=six_bar(theta1(n1),omega1,ii);
s3(n1)=theta(1);
theta3(n1)=theta(2);
theta4(n1)=theta(3);
sE(n1)=theta(4);
v2(n1)=omega(1);
omega3(n1)=omega(2);
omega4(n1)=omega(3);
vE(n1)=omega(4);
a2(n1)=alpha(1);
alpha3(n1)=alpha(2);
alpha4(n1)=alpha(3);
aE(n1)=alpha(4);
end
%3.位移,角速度，角加速度和牛头刨床图形输出
figure(1);
n1=1:459;
t=(n1-1)*2*pi/360;
subplot(2,2,1);%绘角位移及位移线图
plot(t,theta3*du,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角位移及位移线图’)
axis auto;
[haxes,~,~]=plotyy(t,theta4*du,t,sE); %
grid on;
hold on;
xlabel(‘时间/s’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角位移/circ’);
axes(haxes(2));
ylabel(‘位移/m’);
hold on;
grid on;
text(1.15,-0.65,’theta_3′);
text(3.4,0.27,’theta_4′);
text(2.25,-0.15,’s_E’);
subplot(2,2,2); %绘角速度及速度线图
plot(t,omega3,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角速度及速度线图’)
axis auto;
[haxes,~,~]=plotyy(t,omega4,t,vE); %
grid on;
hold on;
xlabel(‘时间/s’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角速度/radcdots∧{-1}’);
axes(haxes(2));
ylabel(‘速度/mcdots∧{-1}’);
hold on;
grid on;
text(3.1,0.35,’omega_3′);
text(2.1,0.1,’omega_4′);
text(5.5,0.45,’v_E’);
subplot(2,2,3); %绘角加速度及加速度线图
plot(t,alpha3,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角加速度及加速度线图’)
axis auto;
[haxes,~,~]=plotyy(t,alpha4,t,aE); %
grid on;
hold on;
xlabel(‘时间’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角加速度/radcdots∧{-2}’);
axes(haxes(2));
ylabel(‘加速度/mcdots∧{-2}’);
hold on;
grid on;
text(1.5,0.3,’alpha_3′);
text(3.5,0.51,’alpha_4′);
text(1.5,-0.11,’a_E’);
subplot(2,2,4); %牛头刨床机构
n1=20;
x(1)=0;y(1)=0;
x(2)=(s3(n1)*1000-50)*cos(theta3(n1));
y(2)=(s3(n1)*1000-50)*sin(theta3(n1));
x(3)=0;y(3)=i6*1000;
x(4)=i1*1000*cos(theta1(n1));
y(4)=s3(n1)*1000*sin(theta3(n1));
x(5)=(s3(n1)*1000+50)*cos(theta3(n1));
y(5)=(s3(n1)*1000+50)*sin(theta3(n1));
x(6)=i3*1000*cos(theta3(n1));
y(6)=i3*1000*sin(theta3(n1));
x(7)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1));
y(7)=i3*1000*sin(theta3(n1))+i4*1000*sin(theta4(n1));
x(8)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))-900;
y(8)=i61*1000;
x(9)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))+600;
y(9)=i61*1000;
x(10)=(s3(n1)*1000-50)*cos(theta3(n1));
y(10)=(s3(n1)*1000-50)*sin(theta3(n1));
x(11)=x(10)+25*cos(pi/2-theta3(n1));
y(11)=y(10)-25*sin(pi/2-theta3(n1));
x(12)=x(11)+100*cos(theta3(n1));
y(12)=y(11)+100*sin(theta3(n1));
x(13)=x(12)-50*cos(pi/2-theta3(n1));
y(13)=y(12)+50*sin(pi/2-theta3(n1));
x(14)=x(10)-25*cos(pi/2-theta3(n1));
y(14)=y(10)+25*sin(pi/2-theta3(n1));
x(15)=x(10);
y(15)=y(10);
x(16)=0;
y(16)=0;
x(17)=0;
y(17)=i6*1000;
k=1:2;
plot(x(k),y(k));
hold on;
k=3:4;
plot(x(k),y(k));
hold on;
k=5:9;
plot(x(k),y(k));
hold on;
k=10:15;
plot(x(k),y(k));
hold on;
k=16:17;
plot(x(k),y(k));
hold on;
grid on;
axis ([-500 600 0 650]);
title(‘牛头刨床’);
grid on;
xlabel(‘mm’);
ylabel(‘mm’);
plot(x(1),y(1),’o’);
plot(x(3),y(3),’o’);
plot(x(4),y(4),’o’);
plot(x(6),y(6),’o’);
plot(x(7),y(7),’o’);
hold on;
grid on;
xlabel(‘mm’);
ylabel(‘mm’);
axis([-400 600 0 650]);
function [theta, omega, alpha]=six_bar(theta1,omega1,ii)
i1=ii(1);
i3=ii(2);
i4=ii(3);
i6=ii(4);
i61=ii(5);
%1.计算角位移和线位移
s3 =sqrt((i1*cos(theta1))*(i1*cos(theta1))+(i6+i1*sin(theta1))*(i6+i1*sin(theta1)));
%s3表示滑块2相对于CD杆的位移
theta3 =acos((i1*cos(theta1 ))/s3 );
%theta3表示杆3转过角度
theta4 =pi-asin((i61-i3*sin(theta3 ))/i4);
%theta4表示杆4转过角度
sE =i3*cos(theta3 )+i4*cos(theta4 );
%sE表示杆5的位移
theta(1)=s3;
theta(2)=theta3;
theta(3)=theta4;
theta(4)=sE;
%2.计算角速度和线速度
A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
-cos(theta3 ),s3 *sin(theta3 ),0,0;
0,i3*sin(theta3 ),i4*sin(theta4 ),1;
0,i3*cos(theta3 ),i4*cos(theta4 ),0];
B=[i1*cos(theta1 );i1*sin(theta1 );0;0];
omega=A(omega1*B);
v2 =omega(1);
%滑块2的速度
omega3 =omega(2);
%构件3的角速度
omega4 =omega(3);
%构件4的角速度
vE =omega(4);
%构件5的速度
%3.计算角加速度和加速度
A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
%从动件位置参数矩阵
cos(theta3 ),-s3 *sin(theta3 ),0,0;
0,i3*sin(theta3 ),i4*sin(theta4 ),1;
0,i3*cos(theta3 ),i4*cos(theta4 ),0];
At=[omega3 *cos(theta3 ),(v2 *cos(theta3 )-s3*omega3 *sin(theta3 )),0,0;
-omega3 *sin(theta3 ),(-v2 *sin(theta3 )-s3 *omega3 *cos(theta3 )),0,0;
0,i3*omega3 *cos(theta3 ),i4*omega4 *cos(theta4 ),0;
0,-i3*omega3 *sin(theta3 ),-i4*omega4 *sin(theta4 ),0];
Bt=[-i1*omega1*sin(theta1 );-i1*omega1*cos(theta1 );0;0];
alpha=A(-At*omega+omega1*Bt);
%机构从动件的加速度列阵
a2 =alpha(1);
%a2表示滑块2的加速度
alpha3 =alpha(2);
%alpha3表示杆3的角加速度
alpha4 =alpha(3);
%alpha4表示杆4的角加速度
aE =alpha(4);
%构件5的加速度
end clear;
i1=0.1;
i3=0.59;
i4=0.159;
i6=0.36;
i61=0.578;
omega1=0.87;
alpha=0;
hd=pi/180;
du=180/pi;
%调用子函数,six_bar返回牛头刨床的位移,角速度,角加速度
for n1 = 1:459
theta1(n1)=-2*pi+5.8119+(n1-1)*hd;
ii=[i1,i3,i4,i6,i61];
[theta,omega,alpha]=six_bar(theta1(n1),omega1,ii);
s3(n1)=theta(1);
theta3(n1)=theta(2);
theta4(n1)=theta(3);
sE(n1)=theta(4);
v2(n1)=omega(1);
omega3(n1)=omega(2);
omega4(n1)=omega(3);
vE(n1)=omega(4);
a2(n1)=alpha(1);
alpha3(n1)=alpha(2);
alpha4(n1)=alpha(3);
aE(n1)=alpha(4);
end
%3.位移,角速度，角加速度和牛头刨床图形输出
figure(1);
n1=1:459;
t=(n1-1)*2*pi/360;
subplot(2,2,1);%绘角位移及位移线图
plot(t,theta3*du,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角位移及位移线图’)
axis auto;
[haxes,~,~]=plotyy(t,theta4*du,t,sE); %
grid on;
hold on;
xlabel(‘时间/s’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角位移/circ’);
axes(haxes(2));
ylabel(‘位移/m’);
hold on;
grid on;
text(1.15,-0.65,’theta_3′);
text(3.4,0.27,’theta_4′);
text(2.25,-0.15,’s_E’);
subplot(2,2,2); %绘角速度及速度线图
plot(t,omega3,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角速度及速度线图’)
axis auto;
[haxes,~,~]=plotyy(t,omega4,t,vE); %
grid on;
hold on;
xlabel(‘时间/s’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角速度/radcdots∧{-1}’);
axes(haxes(2));
ylabel(‘速度/mcdots∧{-1}’);
hold on;
grid on;
text(3.1,0.35,’omega_3′);
text(2.1,0.1,’omega_4′);
text(5.5,0.45,’v_E’);
subplot(2,2,3); %绘角加速度及加速度线图
plot(t,alpha3,’r-.’); %’r-.’表示线条是红色的，并且线型是点画线
grid on;
hold on;
title(‘角加速度及加速度线图’)
axis auto;
[haxes,~,~]=plotyy(t,alpha4,t,aE); %
grid on;
hold on;
xlabel(‘时间’);
axes(haxes(1));%haxes(1)表示左边的坐标轴的句柄,这行代码表示下面的代码是在左边的坐标轴上作图,修改
ylabel(‘角加速度/radcdots∧{-2}’);
axes(haxes(2));
ylabel(‘加速度/mcdots∧{-2}’);
hold on;
grid on;
text(1.5,0.3,’alpha_3′);
text(3.5,0.51,’alpha_4′);
text(1.5,-0.11,’a_E’);
subplot(2,2,4); %牛头刨床机构
n1=20;
x(1)=0;y(1)=0;
x(2)=(s3(n1)*1000-50)*cos(theta3(n1));
y(2)=(s3(n1)*1000-50)*sin(theta3(n1));
x(3)=0;y(3)=i6*1000;
x(4)=i1*1000*cos(theta1(n1));
y(4)=s3(n1)*1000*sin(theta3(n1));
x(5)=(s3(n1)*1000+50)*cos(theta3(n1));
y(5)=(s3(n1)*1000+50)*sin(theta3(n1));
x(6)=i3*1000*cos(theta3(n1));
y(6)=i3*1000*sin(theta3(n1));
x(7)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1));
y(7)=i3*1000*sin(theta3(n1))+i4*1000*sin(theta4(n1));
x(8)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))-900;
y(8)=i61*1000;
x(9)=i3*1000*cos(theta3(n1))+i4*1000*cos(theta4(n1))+600;
y(9)=i61*1000;
x(10)=(s3(n1)*1000-50)*cos(theta3(n1));
y(10)=(s3(n1)*1000-50)*sin(theta3(n1));
x(11)=x(10)+25*cos(pi/2-theta3(n1));
y(11)=y(10)-25*sin(pi/2-theta3(n1));
x(12)=x(11)+100*cos(theta3(n1));
y(12)=y(11)+100*sin(theta3(n1));
x(13)=x(12)-50*cos(pi/2-theta3(n1));
y(13)=y(12)+50*sin(pi/2-theta3(n1));
x(14)=x(10)-25*cos(pi/2-theta3(n1));
y(14)=y(10)+25*sin(pi/2-theta3(n1));
x(15)=x(10);
y(15)=y(10);
x(16)=0;
y(16)=0;
x(17)=0;
y(17)=i6*1000;
k=1:2;
plot(x(k),y(k));
hold on;
k=3:4;
plot(x(k),y(k));
hold on;
k=5:9;
plot(x(k),y(k));
hold on;
k=10:15;
plot(x(k),y(k));
hold on;
k=16:17;
plot(x(k),y(k));
hold on;
grid on;
axis ([-500 600 0 650]);
title(‘牛头刨床’);
grid on;
xlabel(‘mm’);
ylabel(‘mm’);
plot(x(1),y(1),’o’);
plot(x(3),y(3),’o’);
plot(x(4),y(4),’o’);
plot(x(6),y(6),’o’);
plot(x(7),y(7),’o’);
hold on;
grid on;
xlabel(‘mm’);
ylabel(‘mm’);
axis([-400 600 0 650]);
function [theta, omega, alpha]=six_bar(theta1,omega1,ii)
i1=ii(1);
i3=ii(2);
i4=ii(3);
i6=ii(4);
i61=ii(5);
%1.计算角位移和线位移
s3 =sqrt((i1*cos(theta1))*(i1*cos(theta1))+(i6+i1*sin(theta1))*(i6+i1*sin(theta1)));
%s3表示滑块2相对于CD杆的位移
theta3 =acos((i1*cos(theta1 ))/s3 );
%theta3表示杆3转过角度
theta4 =pi-asin((i61-i3*sin(theta3 ))/i4);
%theta4表示杆4转过角度
sE =i3*cos(theta3 )+i4*cos(theta4 );
%sE表示杆5的位移
theta(1)=s3;
theta(2)=theta3;
theta(3)=theta4;
theta(4)=sE;
%2.计算角速度和线速度
A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
-cos(theta3 ),s3 *sin(theta3 ),0,0;
0,i3*sin(theta3 ),i4*sin(theta4 ),1;
0,i3*cos(theta3 ),i4*cos(theta4 ),0];
B=[i1*cos(theta1 );i1*sin(theta1 );0;0];
omega=A(omega1*B);
v2 =omega(1);
%滑块2的速度
omega3 =omega(2);
%构件3的角速度
omega4 =omega(3);
%构件4的角速度
vE =omega(4);
%构件5的速度
%3.计算角加速度和加速度
A=[sin(theta3 ),s3 *cos(theta3 ),0,0;
%从动件位置参数矩阵
cos(theta3 ),-s3 *sin(theta3 ),0,0;
0,i3*sin(theta3 ),i4*sin(theta4 ),1;
0,i3*cos(theta3 ),i4*cos(theta4 ),0];
At=[omega3 *cos(theta3 ),(v2 *cos(theta3 )-s3*omega3 *sin(theta3 )),0,0;
-omega3 *sin(theta3 ),(-v2 *sin(theta3 )-s3 *omega3 *cos(theta3 )),0,0;
0,i3*omega3 *cos(theta3 ),i4*omega4 *cos(theta4 ),0;
0,-i3*omega3 *sin(theta3 ),-i4*omega4 *sin(theta4 ),0];
Bt=[-i1*omega1*sin(theta1 );-i1*omega1*cos(theta1 );0;0];
alpha=A(-At*omega+omega1*Bt);
%机构从动件的加速度列阵
a2 =alpha(1);
%a2表示滑块2的加速度
alpha3 =alpha(2);
%alpha3表示杆3的角加速度
alpha4 =alpha(3);
%alpha4表示杆4的角加速度
aE =alpha(4);
%构件5的加速度
end six_bar MATLAB Answers — New Questions
Does the Digilent Toolbox Add-On work with Matlab Library Compiler?
Does the Digilent Toolbox work with the Matlab Library compiler? I’m attempting to deploy a matlab .NET application, but when the deployed application runs the command
AO = daq("digilent")
it returns the exception "Specify an operational vendor" which tells me it doesn’t have the Digilent Toolbox info included in the compiled application.
I also don’t see Digilent Toolbox in the Suggested Support Packages in the Library Compiler application window.
Here is a link to the toolbox I’m referring to: https://www.mathworks.com/matlabcentral/fileexchange/122817-digilent-toolboxDoes the Digilent Toolbox work with the Matlab Library compiler? I’m attempting to deploy a matlab .NET application, but when the deployed application runs the command
AO = daq("digilent")
it returns the exception "Specify an operational vendor" which tells me it doesn’t have the Digilent Toolbox info included in the compiled application.
I also don’t see Digilent Toolbox in the Suggested Support Packages in the Library Compiler application window.
Here is a link to the toolbox I’m referring to: https://www.mathworks.com/matlabcentral/fileexchange/122817-digilent-toolbox Does the Digilent Toolbox work with the Matlab Library compiler? I’m attempting to deploy a matlab .NET application, but when the deployed application runs the command
AO = daq("digilent")
it returns the exception "Specify an operational vendor" which tells me it doesn’t have the Digilent Toolbox info included in the compiled application.
I also don’t see Digilent Toolbox in the Suggested Support Packages in the Library Compiler application window.
Here is a link to the toolbox I’m referring to: https://www.mathworks.com/matlabcentral/fileexchange/122817-digilent-toolbox digilent, compiler, toolbox MATLAB Answers — New Questions
How to use fminsearch for vectors
I searched for prior discussions on fminsearch, including one on "How to use fminsearch for a vector" and did not find what I am looking for.
I have two somewhat independent vectors of measured values, A and B. Is the optimisation a linear domain or a plane? I have had a quick look at optimvar, fminsearch, fmincon, and fminbnd and have not secured a result. Error on fminsearch is " Maximum number of function evaluations has been exceeded". Is there a more appropriate function than fminsearch? Thank you.
rng(‘shuffle’);
A = rand(1,10); % proxy for measured values
B = rand(1,10); % proxy for independent measured values
fun = @(ratio) min(sum(A.*ratio – B));
ratioRange = [1.0,1.2];
ratio = fminsearch(fun,ratioRange);I searched for prior discussions on fminsearch, including one on "How to use fminsearch for a vector" and did not find what I am looking for.
I have two somewhat independent vectors of measured values, A and B. Is the optimisation a linear domain or a plane? I have had a quick look at optimvar, fminsearch, fmincon, and fminbnd and have not secured a result. Error on fminsearch is " Maximum number of function evaluations has been exceeded". Is there a more appropriate function than fminsearch? Thank you.
rng(‘shuffle’);
A = rand(1,10); % proxy for measured values
B = rand(1,10); % proxy for independent measured values
fun = @(ratio) min(sum(A.*ratio – B));
ratioRange = [1.0,1.2];
ratio = fminsearch(fun,ratioRange); I searched for prior discussions on fminsearch, including one on "How to use fminsearch for a vector" and did not find what I am looking for.
I have two somewhat independent vectors of measured values, A and B. Is the optimisation a linear domain or a plane? I have had a quick look at optimvar, fminsearch, fmincon, and fminbnd and have not secured a result. Error on fminsearch is " Maximum number of function evaluations has been exceeded". Is there a more appropriate function than fminsearch? Thank you.
rng(‘shuffle’);
A = rand(1,10); % proxy for measured values
B = rand(1,10); % proxy for independent measured values
fun = @(ratio) min(sum(A.*ratio – B));
ratioRange = [1.0,1.2];
ratio = fminsearch(fun,ratioRange); fminsearch, vectors MATLAB Answers — New Questions
Robotic System Toolbox workspaceGoalRegion prevents me to make a right Pick and Place with a robot manipualtor, surrounded by collision object in the scenario.
I’m using the Robotic System Toolbox in order to plan the path of a Robotic Manipulator in a scenario full of collision objects. I’m using the function manipulatorRRT, then in the function ‘plan’ it asks me to put the start and goal configurations of the robot. I’m using the function ‘workspaceGoalRegion’ instead of giving the Joint Configuration matrix, because using ‘GoalRegion.ReferencePose=trvec2tform[ x y z ]’ I can calculate the Joint configuration of the robot starting from a point I want to reach in the space.
The problem is that it works until a collision object is in front of the robot. In that situation matlab starts a neverending loop and the programme never stops.
Here is the script I am talking about:
scenario = robotScenario(UpdateRate=10);
addMesh(scenario,"Box",Position=[0.35 0 0.2],Size=[0.5 0.9 0.05],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0 -0.6 0.2],Size=[1.3 0.4 0.235],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0.3 -0.25 0.4],Size=[0.8 0.03 0.35],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0.52 0 0.278],Size=[0.06 0.06 0.1],Color=[1 0.5 0.25],Collision="mesh")
addMesh(scenario,"Box",Position=[0.4 -0.1 0.2758],Size=[0.06 0.06 0.1],Color=[1 0.5 0.25],Collision="mesh")
visuals = "on";
collisions = "off";
show3D(scenario,Visuals=visuals,Collisions=collisions);
title("Scene with collision object-based static meshes")
view(81,19)
light
franka = loadrobot("frankaEmikaPanda");
robot = robotPlatform("Manipulator",scenario,RigidBodyTree=franka,ReferenceFrame="ENU");
initialConfig = homeConfiguration(robot.RigidBodyTree);
initialConfig(7) = 0.65;
initialConfig(6) = 0.75;
boxToMove = robotPlatform("Box",scenario,Collision="mesh",InitialBasePosition=[0.5 0.15 0.278]);
updateMesh(boxToMove,"Cuboid",Collision="mesh",Size=[0.06 0.06 0.1],Color=[0.9 0.0 0.0]);
setup(scenario)
move(robot,"joint",initialConfig)
ax = show3D(scenario,Visuals=visuals,Collisions=collisions);
title("Scene with static meshes and robot platform")
view(81,19)
light
goalRegion = workspaceGoalRegion(franka.BodyNames{end});
goalRegion.ReferencePose=trvec2tform([0.5 0.15 0.278+0.05]);
goalRegion.Bounds(1, 🙂 = [0 0]; % X Bounds
goalRegion.Bounds(2, 🙂 = [0 0]; % Y Bounds
goalRegion.Bounds(4, 🙂 = [-pi/2 pi/2];
goalRegion.ReferencePose=goalRegion.ReferencePose*eul2tform([0 pi 0],"ZYX");
pickUpConfig=goalRegion;
%%
planner = manipulatorRRT(robot.RigidBodyTree,scenario.CollisionMeshes);
planner.MaxConnectionDistance = 0.6;
planner.ValidationDistance = 0.01;
planner.SkippedSelfCollisions = "parent";
planner.IgnoreSelfCollision = true;
rng("default")
path = plan(planner,initialConfig,pickUpConfig);
% Number of interpolations between each adjacent configuration.
numState = 10;
interpStates = interpolate(planner,path,numState);
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Approaching towards pick-up location")
for idx = 1:size(interpStates,1)
% Get current joint configuration.
jointConfig = interpStates(idx,:);
% Move manipulator with current joint configuration.
move(robot,"joint",jointConfig)
% Use fast update to move platform visualization frames.
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
%Draw trajectory.
drawnow
% Advance scenario simulation time.
advance(scenario);
end
attach(robot,"Box","panda_hand",ChildToParentTransform=trvec2tform([0.0 0.0 0.1]))
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Picked-up target box")
startConfig=path(end,:);
goalRegion.ReferencePose=trvec2tform([0.4 -0.6 0.4]);
goalRegion.ReferencePose=goalRegion.ReferencePose*eul2tform([0 pi 0],"ZYX");
dropConfig=goalRegion;
planner = manipulatorRRT(robot.RigidBodyTree,scenario.CollisionMeshes);
planner.MaxConnectionDistance = 0.6;
planner.ValidationDistance = 0.01;
planner.SkippedSelfCollisions = "parent";
planner.IgnoreSelfCollision = true;
path = plan(planner,startConfig,dropConfig);
interpStates = interpolate(planner,path,numState);
% Visualize scene
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Approaching towards drop location")
for idx = 1:size(interpStates,1)
% Get current joint configuration.
jointConfig = interpStates(idx,:);
% Move manipulator with current joint configuration.
move(robot,"joint",jointConfig)
% Use fast update to move platform visualization frames.
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
%Draw trajectory.
drawnow
% Advance scenario simulation time.
advance(scenario);
end
detach(robot)
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Placed target box at drop location")
read(robot,"joint")I’m using the Robotic System Toolbox in order to plan the path of a Robotic Manipulator in a scenario full of collision objects. I’m using the function manipulatorRRT, then in the function ‘plan’ it asks me to put the start and goal configurations of the robot. I’m using the function ‘workspaceGoalRegion’ instead of giving the Joint Configuration matrix, because using ‘GoalRegion.ReferencePose=trvec2tform[ x y z ]’ I can calculate the Joint configuration of the robot starting from a point I want to reach in the space.
The problem is that it works until a collision object is in front of the robot. In that situation matlab starts a neverending loop and the programme never stops.
Here is the script I am talking about:
scenario = robotScenario(UpdateRate=10);
addMesh(scenario,"Box",Position=[0.35 0 0.2],Size=[0.5 0.9 0.05],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0 -0.6 0.2],Size=[1.3 0.4 0.235],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0.3 -0.25 0.4],Size=[0.8 0.03 0.35],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0.52 0 0.278],Size=[0.06 0.06 0.1],Color=[1 0.5 0.25],Collision="mesh")
addMesh(scenario,"Box",Position=[0.4 -0.1 0.2758],Size=[0.06 0.06 0.1],Color=[1 0.5 0.25],Collision="mesh")
visuals = "on";
collisions = "off";
show3D(scenario,Visuals=visuals,Collisions=collisions);
title("Scene with collision object-based static meshes")
view(81,19)
light
franka = loadrobot("frankaEmikaPanda");
robot = robotPlatform("Manipulator",scenario,RigidBodyTree=franka,ReferenceFrame="ENU");
initialConfig = homeConfiguration(robot.RigidBodyTree);
initialConfig(7) = 0.65;
initialConfig(6) = 0.75;
boxToMove = robotPlatform("Box",scenario,Collision="mesh",InitialBasePosition=[0.5 0.15 0.278]);
updateMesh(boxToMove,"Cuboid",Collision="mesh",Size=[0.06 0.06 0.1],Color=[0.9 0.0 0.0]);
setup(scenario)
move(robot,"joint",initialConfig)
ax = show3D(scenario,Visuals=visuals,Collisions=collisions);
title("Scene with static meshes and robot platform")
view(81,19)
light
goalRegion = workspaceGoalRegion(franka.BodyNames{end});
goalRegion.ReferencePose=trvec2tform([0.5 0.15 0.278+0.05]);
goalRegion.Bounds(1, 🙂 = [0 0]; % X Bounds
goalRegion.Bounds(2, 🙂 = [0 0]; % Y Bounds
goalRegion.Bounds(4, 🙂 = [-pi/2 pi/2];
goalRegion.ReferencePose=goalRegion.ReferencePose*eul2tform([0 pi 0],"ZYX");
pickUpConfig=goalRegion;
%%
planner = manipulatorRRT(robot.RigidBodyTree,scenario.CollisionMeshes);
planner.MaxConnectionDistance = 0.6;
planner.ValidationDistance = 0.01;
planner.SkippedSelfCollisions = "parent";
planner.IgnoreSelfCollision = true;
rng("default")
path = plan(planner,initialConfig,pickUpConfig);
% Number of interpolations between each adjacent configuration.
numState = 10;
interpStates = interpolate(planner,path,numState);
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Approaching towards pick-up location")
for idx = 1:size(interpStates,1)
% Get current joint configuration.
jointConfig = interpStates(idx,:);
% Move manipulator with current joint configuration.
move(robot,"joint",jointConfig)
% Use fast update to move platform visualization frames.
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
%Draw trajectory.
drawnow
% Advance scenario simulation time.
advance(scenario);
end
attach(robot,"Box","panda_hand",ChildToParentTransform=trvec2tform([0.0 0.0 0.1]))
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Picked-up target box")
startConfig=path(end,:);
goalRegion.ReferencePose=trvec2tform([0.4 -0.6 0.4]);
goalRegion.ReferencePose=goalRegion.ReferencePose*eul2tform([0 pi 0],"ZYX");
dropConfig=goalRegion;
planner = manipulatorRRT(robot.RigidBodyTree,scenario.CollisionMeshes);
planner.MaxConnectionDistance = 0.6;
planner.ValidationDistance = 0.01;
planner.SkippedSelfCollisions = "parent";
planner.IgnoreSelfCollision = true;
path = plan(planner,startConfig,dropConfig);
interpStates = interpolate(planner,path,numState);
% Visualize scene
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Approaching towards drop location")
for idx = 1:size(interpStates,1)
% Get current joint configuration.
jointConfig = interpStates(idx,:);
% Move manipulator with current joint configuration.
move(robot,"joint",jointConfig)
% Use fast update to move platform visualization frames.
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
%Draw trajectory.
drawnow
% Advance scenario simulation time.
advance(scenario);
end
detach(robot)
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Placed target box at drop location")
read(robot,"joint") I’m using the Robotic System Toolbox in order to plan the path of a Robotic Manipulator in a scenario full of collision objects. I’m using the function manipulatorRRT, then in the function ‘plan’ it asks me to put the start and goal configurations of the robot. I’m using the function ‘workspaceGoalRegion’ instead of giving the Joint Configuration matrix, because using ‘GoalRegion.ReferencePose=trvec2tform[ x y z ]’ I can calculate the Joint configuration of the robot starting from a point I want to reach in the space.
The problem is that it works until a collision object is in front of the robot. In that situation matlab starts a neverending loop and the programme never stops.
Here is the script I am talking about:
scenario = robotScenario(UpdateRate=10);
addMesh(scenario,"Box",Position=[0.35 0 0.2],Size=[0.5 0.9 0.05],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0 -0.6 0.2],Size=[1.3 0.4 0.235],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0.3 -0.25 0.4],Size=[0.8 0.03 0.35],Color=[0.7 0.7 0.7],Collision="mesh")
addMesh(scenario,"Box",Position=[0.52 0 0.278],Size=[0.06 0.06 0.1],Color=[1 0.5 0.25],Collision="mesh")
addMesh(scenario,"Box",Position=[0.4 -0.1 0.2758],Size=[0.06 0.06 0.1],Color=[1 0.5 0.25],Collision="mesh")
visuals = "on";
collisions = "off";
show3D(scenario,Visuals=visuals,Collisions=collisions);
title("Scene with collision object-based static meshes")
view(81,19)
light
franka = loadrobot("frankaEmikaPanda");
robot = robotPlatform("Manipulator",scenario,RigidBodyTree=franka,ReferenceFrame="ENU");
initialConfig = homeConfiguration(robot.RigidBodyTree);
initialConfig(7) = 0.65;
initialConfig(6) = 0.75;
boxToMove = robotPlatform("Box",scenario,Collision="mesh",InitialBasePosition=[0.5 0.15 0.278]);
updateMesh(boxToMove,"Cuboid",Collision="mesh",Size=[0.06 0.06 0.1],Color=[0.9 0.0 0.0]);
setup(scenario)
move(robot,"joint",initialConfig)
ax = show3D(scenario,Visuals=visuals,Collisions=collisions);
title("Scene with static meshes and robot platform")
view(81,19)
light
goalRegion = workspaceGoalRegion(franka.BodyNames{end});
goalRegion.ReferencePose=trvec2tform([0.5 0.15 0.278+0.05]);
goalRegion.Bounds(1, 🙂 = [0 0]; % X Bounds
goalRegion.Bounds(2, 🙂 = [0 0]; % Y Bounds
goalRegion.Bounds(4, 🙂 = [-pi/2 pi/2];
goalRegion.ReferencePose=goalRegion.ReferencePose*eul2tform([0 pi 0],"ZYX");
pickUpConfig=goalRegion;
%%
planner = manipulatorRRT(robot.RigidBodyTree,scenario.CollisionMeshes);
planner.MaxConnectionDistance = 0.6;
planner.ValidationDistance = 0.01;
planner.SkippedSelfCollisions = "parent";
planner.IgnoreSelfCollision = true;
rng("default")
path = plan(planner,initialConfig,pickUpConfig);
% Number of interpolations between each adjacent configuration.
numState = 10;
interpStates = interpolate(planner,path,numState);
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Approaching towards pick-up location")
for idx = 1:size(interpStates,1)
% Get current joint configuration.
jointConfig = interpStates(idx,:);
% Move manipulator with current joint configuration.
move(robot,"joint",jointConfig)
% Use fast update to move platform visualization frames.
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
%Draw trajectory.
drawnow
% Advance scenario simulation time.
advance(scenario);
end
attach(robot,"Box","panda_hand",ChildToParentTransform=trvec2tform([0.0 0.0 0.1]))
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Picked-up target box")
startConfig=path(end,:);
goalRegion.ReferencePose=trvec2tform([0.4 -0.6 0.4]);
goalRegion.ReferencePose=goalRegion.ReferencePose*eul2tform([0 pi 0],"ZYX");
dropConfig=goalRegion;
planner = manipulatorRRT(robot.RigidBodyTree,scenario.CollisionMeshes);
planner.MaxConnectionDistance = 0.6;
planner.ValidationDistance = 0.01;
planner.SkippedSelfCollisions = "parent";
planner.IgnoreSelfCollision = true;
path = plan(planner,startConfig,dropConfig);
interpStates = interpolate(planner,path,numState);
% Visualize scene
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Approaching towards drop location")
for idx = 1:size(interpStates,1)
% Get current joint configuration.
jointConfig = interpStates(idx,:);
% Move manipulator with current joint configuration.
move(robot,"joint",jointConfig)
% Use fast update to move platform visualization frames.
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
%Draw trajectory.
drawnow
% Advance scenario simulation time.
advance(scenario);
end
detach(robot)
show3D(scenario,FastUpdate=true,Parent=ax,Visuals=visuals,Collisions=collisions);
title("Placed target box at drop location")
read(robot,"joint") workspacegoalregion, manipulatorrrt, plan MATLAB Answers — New Questions
Is it possible to define a distance from obstacles to RRT*?
Hello everyone,
I wanted to know if it is possible to define a minimum distance from an obstacle for RRT*? I have a code with RRT* and I noticed that it provides me with routes that are very close to obstacles, this creates problems for me to continue my project
Is it possible to define a variable that prevents the algorithm from sticking to the obstacle?
I can’t solve this problem.
I would really appreciate help with this
Thank you very much everyone!! 🙂
I am attaching a figure to illustrate the trajectory given to me by RRT* (I made it discrete)Hello everyone,
I wanted to know if it is possible to define a minimum distance from an obstacle for RRT*? I have a code with RRT* and I noticed that it provides me with routes that are very close to obstacles, this creates problems for me to continue my project
Is it possible to define a variable that prevents the algorithm from sticking to the obstacle?
I can’t solve this problem.
I would really appreciate help with this
Thank you very much everyone!! 🙂
I am attaching a figure to illustrate the trajectory given to me by RRT* (I made it discrete) Hello everyone,
I wanted to know if it is possible to define a minimum distance from an obstacle for RRT*? I have a code with RRT* and I noticed that it provides me with routes that are very close to obstacles, this creates problems for me to continue my project
Is it possible to define a variable that prevents the algorithm from sticking to the obstacle?
I can’t solve this problem.
I would really appreciate help with this
Thank you very much everyone!! 🙂
I am attaching a figure to illustrate the trajectory given to me by RRT* (I made it discrete) rrt*, trajectory, maze, matrix, distance, path planner MATLAB Answers — New Questions
Modelling shell behavior with femodel
I apologize if this question is banal but is there a way to model bending plates with femodel?
Playing around with the object definition I managed to make a 2D rectangular model with uniform surface pressure and fixed edges, I assigned material properties, generated a mesh and solved the model. When I inspect the results, there are only null fields that correspond to membrane behavior and no bending.I apologize if this question is banal but is there a way to model bending plates with femodel?
Playing around with the object definition I managed to make a 2D rectangular model with uniform surface pressure and fixed edges, I assigned material properties, generated a mesh and solved the model. When I inspect the results, there are only null fields that correspond to membrane behavior and no bending. I apologize if this question is banal but is there a way to model bending plates with femodel?
Playing around with the object definition I managed to make a 2D rectangular model with uniform surface pressure and fixed edges, I assigned material properties, generated a mesh and solved the model. When I inspect the results, there are only null fields that correspond to membrane behavior and no bending. femodel, structural static, bending plate MATLAB Answers — New Questions
Scroll through 1 dimension of plotted Array in one figure
Dear all i want to plot the Data in an array and overlay a vektor plot with quiver the Array for the main plot has the deimension 167x176x2000. Two additional Arrays of the same size contain the information for the quiver overlay. I do display the data as XY plots at a givven z value. Obviously i do not want to store 2000 Images. Instead i want a plot, were i can scroll through the 3d dimension of my Array either with the mouse wheel or up and down arrows. The values for the quiver need to be updated as well.
Here is my code wich currently only displays the last image.
quiver_factor=1000000000; %factor to increase the arrows on the quiver images WAS 100 BEFORE
colorscale=0.3; %max abs lateral force in nN. % 0.12 fo const Fz % ORIGINALLY 0.12
for k=1:size(Main_array,3)
clear imagex imagey
[X,Y]=meshgrid(1:size(2nd_array,1),1:size(2nd_array,2));
imagex(:,:)=2nd_array(:,:,k);
imagex=imagex’*quiver_factor;
imagey(:,:)=3nd_array(:,:,k);
imagey=imagey’*quiver_factor;
image_F=Main_array(:,:,k)’;
contourf(X,Y,image_F, 512, ‘LineColor’, ‘none’);
colormap(hsv(1024));
hold on
quiver(X,Y,imagex,imagey);
hold off
end
I found the follwing link, but could make the solution work with my code:
<https://de.mathworks.com/matlabcentral/answers/304335-how-can-i-use-up-and-down-arrow-keys-to-scroll-through-my-plots-in-matlab-similar-to-a-query-on-21 Solution form this forum>
Thanks n advance!Dear all i want to plot the Data in an array and overlay a vektor plot with quiver the Array for the main plot has the deimension 167x176x2000. Two additional Arrays of the same size contain the information for the quiver overlay. I do display the data as XY plots at a givven z value. Obviously i do not want to store 2000 Images. Instead i want a plot, were i can scroll through the 3d dimension of my Array either with the mouse wheel or up and down arrows. The values for the quiver need to be updated as well.
Here is my code wich currently only displays the last image.
quiver_factor=1000000000; %factor to increase the arrows on the quiver images WAS 100 BEFORE
colorscale=0.3; %max abs lateral force in nN. % 0.12 fo const Fz % ORIGINALLY 0.12
for k=1:size(Main_array,3)
clear imagex imagey
[X,Y]=meshgrid(1:size(2nd_array,1),1:size(2nd_array,2));
imagex(:,:)=2nd_array(:,:,k);
imagex=imagex’*quiver_factor;
imagey(:,:)=3nd_array(:,:,k);
imagey=imagey’*quiver_factor;
image_F=Main_array(:,:,k)’;
contourf(X,Y,image_F, 512, ‘LineColor’, ‘none’);
colormap(hsv(1024));
hold on
quiver(X,Y,imagex,imagey);
hold off
end
I found the follwing link, but could make the solution work with my code:
<https://de.mathworks.com/matlabcentral/answers/304335-how-can-i-use-up-and-down-arrow-keys-to-scroll-through-my-plots-in-matlab-similar-to-a-query-on-21 Solution form this forum>
Thanks n advance! Dear all i want to plot the Data in an array and overlay a vektor plot with quiver the Array for the main plot has the deimension 167x176x2000. Two additional Arrays of the same size contain the information for the quiver overlay. I do display the data as XY plots at a givven z value. Obviously i do not want to store 2000 Images. Instead i want a plot, were i can scroll through the 3d dimension of my Array either with the mouse wheel or up and down arrows. The values for the quiver need to be updated as well.
Here is my code wich currently only displays the last image.
quiver_factor=1000000000; %factor to increase the arrows on the quiver images WAS 100 BEFORE
colorscale=0.3; %max abs lateral force in nN. % 0.12 fo const Fz % ORIGINALLY 0.12
for k=1:size(Main_array,3)
clear imagex imagey
[X,Y]=meshgrid(1:size(2nd_array,1),1:size(2nd_array,2));
imagex(:,:)=2nd_array(:,:,k);
imagex=imagex’*quiver_factor;
imagey(:,:)=3nd_array(:,:,k);
imagey=imagey’*quiver_factor;
image_F=Main_array(:,:,k)’;
contourf(X,Y,image_F, 512, ‘LineColor’, ‘none’);
colormap(hsv(1024));
hold on
quiver(X,Y,imagex,imagey);
hold off
end
I found the follwing link, but could make the solution work with my code:
<https://de.mathworks.com/matlabcentral/answers/304335-how-can-i-use-up-and-down-arrow-keys-to-scroll-through-my-plots-in-matlab-similar-to-a-query-on-21 Solution form this forum>
Thanks n advance! matlab plots scrolling bacwards and forwards MATLAB Answers — New Questions
I am not able to open internally saved test harness, blockdiagram.xml is missing. How do I resolved this?
Unable to load Test Harness ‘ThRnwyDet_TemperaturePreconditioning_Harness’.
Caused by:
File ‘D:SIBE_TE_BMSASW_DevspecswCoreAPSWThRnwyrunnable_10msThRnwyThRnwyDet_Updated.slx’ does not contain the requested part: /simulink/ThRnwyDet_TemperaturePreconditioning_Harness_226/blockdiagram.xml
My test harness was internally saved
It is generated / released on r2023b version, but now I am not able to open itUnable to load Test Harness ‘ThRnwyDet_TemperaturePreconditioning_Harness’.
Caused by:
File ‘D:SIBE_TE_BMSASW_DevspecswCoreAPSWThRnwyrunnable_10msThRnwyThRnwyDet_Updated.slx’ does not contain the requested part: /simulink/ThRnwyDet_TemperaturePreconditioning_Harness_226/blockdiagram.xml
My test harness was internally saved
It is generated / released on r2023b version, but now I am not able to open it Unable to load Test Harness ‘ThRnwyDet_TemperaturePreconditioning_Harness’.
Caused by:
File ‘D:SIBE_TE_BMSASW_DevspecswCoreAPSWThRnwyrunnable_10msThRnwyThRnwyDet_Updated.slx’ does not contain the requested part: /simulink/ThRnwyDet_TemperaturePreconditioning_Harness_226/blockdiagram.xml
My test harness was internally saved
It is generated / released on r2023b version, but now I am not able to open it unknown blockdiagram.xml missing, test harness, simulink MATLAB Answers — New Questions
Multivariable Optimization problem with MATLAB
I wish to find a maximum for the following function, the function is m:
x = normrnd(0,1,[1,5]);
global x
xbar=mean(x);
global xbar
n=length(x);
global n
an = @(a0,b0,k0,m0) a0+n/2;
bn = @(a0,b0,k0,m0) (b0 + ((n-1) * sum((x-xbar).^2))/2 + ((k0 * n * ((xbar-m0).^2) )/(2*(k0 + n))));
kn = @(a0,b0,k0,m0) (k0 + n);
mn = @(a0,b0,k0,m0) ((k0 * m0)+(n.*xbar))/(k0 + n);
m=@(a0,b0,k0,m0)(gamma(an(a0,b0,k0,m0))/gamma(a0)) * ((b0^a0)/(bn(a0,b0,k0,m0).^an(a0,b0,k0,m0))) * (sqrt(k0/kn(a0,b0,k0,m0))) * ((2*pi)^(-n/2))
This is not linear. I wish to find a Global maximum. Most of the function I found in matlab handles the polynomials. But I do not find any particular function that gives me a maximum for this. I tried to use ‘fminunc’. But that is showing me an error message
"Failure in initial objective function evaluation. FMINUNC cannot continue."
Looking forward to some help. Thanks in advance.I wish to find a maximum for the following function, the function is m:
x = normrnd(0,1,[1,5]);
global x
xbar=mean(x);
global xbar
n=length(x);
global n
an = @(a0,b0,k0,m0) a0+n/2;
bn = @(a0,b0,k0,m0) (b0 + ((n-1) * sum((x-xbar).^2))/2 + ((k0 * n * ((xbar-m0).^2) )/(2*(k0 + n))));
kn = @(a0,b0,k0,m0) (k0 + n);
mn = @(a0,b0,k0,m0) ((k0 * m0)+(n.*xbar))/(k0 + n);
m=@(a0,b0,k0,m0)(gamma(an(a0,b0,k0,m0))/gamma(a0)) * ((b0^a0)/(bn(a0,b0,k0,m0).^an(a0,b0,k0,m0))) * (sqrt(k0/kn(a0,b0,k0,m0))) * ((2*pi)^(-n/2))
This is not linear. I wish to find a Global maximum. Most of the function I found in matlab handles the polynomials. But I do not find any particular function that gives me a maximum for this. I tried to use ‘fminunc’. But that is showing me an error message
"Failure in initial objective function evaluation. FMINUNC cannot continue."
Looking forward to some help. Thanks in advance. I wish to find a maximum for the following function, the function is m:
x = normrnd(0,1,[1,5]);
global x
xbar=mean(x);
global xbar
n=length(x);
global n
an = @(a0,b0,k0,m0) a0+n/2;
bn = @(a0,b0,k0,m0) (b0 + ((n-1) * sum((x-xbar).^2))/2 + ((k0 * n * ((xbar-m0).^2) )/(2*(k0 + n))));
kn = @(a0,b0,k0,m0) (k0 + n);
mn = @(a0,b0,k0,m0) ((k0 * m0)+(n.*xbar))/(k0 + n);
m=@(a0,b0,k0,m0)(gamma(an(a0,b0,k0,m0))/gamma(a0)) * ((b0^a0)/(bn(a0,b0,k0,m0).^an(a0,b0,k0,m0))) * (sqrt(k0/kn(a0,b0,k0,m0))) * ((2*pi)^(-n/2))
This is not linear. I wish to find a Global maximum. Most of the function I found in matlab handles the polynomials. But I do not find any particular function that gives me a maximum for this. I tried to use ‘fminunc’. But that is showing me an error message
"Failure in initial objective function evaluation. FMINUNC cannot continue."
Looking forward to some help. Thanks in advance. optimization MATLAB Answers — New Questions
Plot bar graphs on map
Hi, I am wondering how to put 2D bar graph on a map. For example, I have some bar graphs of monthly temperature for several different weather station. And I want to put these graphs on a map corresponding to the location of each weather station. Is there anyway to do that with matlab?Hi, I am wondering how to put 2D bar graph on a map. For example, I have some bar graphs of monthly temperature for several different weather station. And I want to put these graphs on a map corresponding to the location of each weather station. Is there anyway to do that with matlab? Hi, I am wondering how to put 2D bar graph on a map. For example, I have some bar graphs of monthly temperature for several different weather station. And I want to put these graphs on a map corresponding to the location of each weather station. Is there anyway to do that with matlab? plot MATLAB Answers — New Questions
Time format is changed to decimal when calling the time data from excel.
In matlab, when i read a time data (hh:mm:ss AM/PM) from excel, the time data becomes decimal number so unable to plot the time in x-axis in hh:mm:ss format. Kindly guideIn matlab, when i read a time data (hh:mm:ss AM/PM) from excel, the time data becomes decimal number so unable to plot the time in x-axis in hh:mm:ss format. Kindly guide In matlab, when i read a time data (hh:mm:ss AM/PM) from excel, the time data becomes decimal number so unable to plot the time in x-axis in hh:mm:ss format. Kindly guide time format, hh:mm:ss MATLAB Answers — New Questions
importing data from .csv files of different lengths and same columns in matlab and getting RMS.
I have a displacement data in 73 .csv files which are comma separated with 7 (fixed) columns and rows (varying) from 1500 to 9000 rows depending on the file. Each file contains 7 Columns as x y z(coordinates of voxels) and u v w(displacements in u v w directions) and sigma in 7th column. I want to create a code to read the .csv files and get the rms error in u v w directions respectively. RMS formula has to applied and i have to get 1 rms value for all the files combined and that too it should consider all the values of u v w in all the files.Can someone help me with the idea.I have attached the code which i used which i not working. I attached some files of the format .csv in the attachments.
data0001 = csvread(‘C://Users//User//Desktop//.csv//data_0001.csv’);
.
.
.
.
.
data0071 = csvread(‘C://Users//User//Desktop//.csv//data_0071.csv’);
data0072 = csvread(‘C://Users//User//Desktop//.csv//data_0072.csv’);
data0073 = csvread(‘C://Users//User//Desktop//.csv//data_0073.csv’);
%Root Mean Squared Error
for ii= 01:01:73
erroru{ii} = (sum((data00ii(:,4)).^2)); %error in u
errorv{ii} = (sum((data00ii(:,5)).^2)); %error in v
errorw{ii} = (sum((data00ii(:,6)).^2)); %error in w
U=erroru{ii}/(length(data00ii));
V=errorv{ii}/(length(data00ii));
W=errorw{ii}/(length(data00ii));
S= sum(U+V+W);
rmsError{ii} = sqrt(S);
show(rmsError{ii})
endI have a displacement data in 73 .csv files which are comma separated with 7 (fixed) columns and rows (varying) from 1500 to 9000 rows depending on the file. Each file contains 7 Columns as x y z(coordinates of voxels) and u v w(displacements in u v w directions) and sigma in 7th column. I want to create a code to read the .csv files and get the rms error in u v w directions respectively. RMS formula has to applied and i have to get 1 rms value for all the files combined and that too it should consider all the values of u v w in all the files.Can someone help me with the idea.I have attached the code which i used which i not working. I attached some files of the format .csv in the attachments.
data0001 = csvread(‘C://Users//User//Desktop//.csv//data_0001.csv’);
.
.
.
.
.
data0071 = csvread(‘C://Users//User//Desktop//.csv//data_0071.csv’);
data0072 = csvread(‘C://Users//User//Desktop//.csv//data_0072.csv’);
data0073 = csvread(‘C://Users//User//Desktop//.csv//data_0073.csv’);
%Root Mean Squared Error
for ii= 01:01:73
erroru{ii} = (sum((data00ii(:,4)).^2)); %error in u
errorv{ii} = (sum((data00ii(:,5)).^2)); %error in v
errorw{ii} = (sum((data00ii(:,6)).^2)); %error in w
U=erroru{ii}/(length(data00ii));
V=errorv{ii}/(length(data00ii));
W=errorw{ii}/(length(data00ii));
S= sum(U+V+W);
rmsError{ii} = sqrt(S);
show(rmsError{ii})
end I have a displacement data in 73 .csv files which are comma separated with 7 (fixed) columns and rows (varying) from 1500 to 9000 rows depending on the file. Each file contains 7 Columns as x y z(coordinates of voxels) and u v w(displacements in u v w directions) and sigma in 7th column. I want to create a code to read the .csv files and get the rms error in u v w directions respectively. RMS formula has to applied and i have to get 1 rms value for all the files combined and that too it should consider all the values of u v w in all the files.Can someone help me with the idea.I have attached the code which i used which i not working. I attached some files of the format .csv in the attachments.
data0001 = csvread(‘C://Users//User//Desktop//.csv//data_0001.csv’);
.
.
.
.
.
data0071 = csvread(‘C://Users//User//Desktop//.csv//data_0071.csv’);
data0072 = csvread(‘C://Users//User//Desktop//.csv//data_0072.csv’);
data0073 = csvread(‘C://Users//User//Desktop//.csv//data_0073.csv’);
%Root Mean Squared Error
for ii= 01:01:73
erroru{ii} = (sum((data00ii(:,4)).^2)); %error in u
errorv{ii} = (sum((data00ii(:,5)).^2)); %error in v
errorw{ii} = (sum((data00ii(:,6)).^2)); %error in w
U=erroru{ii}/(length(data00ii));
V=errorv{ii}/(length(data00ii));
W=errorw{ii}/(length(data00ii));
S= sum(U+V+W);
rmsError{ii} = sqrt(S);
show(rmsError{ii})
end rms, csv MATLAB Answers — New Questions
How to check if my linear constraints and bounds are infeasible & reduce execution time?
Hello,
I’m making a code where it creates an optimized engine usage schedule based on engine’s reliability using genetic algorithm solver. The following is the code that I’ve made:
tic
clear
% Create arrays to store total Hour & daily result
totalHourWeekday = zeros(5,13);
result = zeros(5,13);
% Calculate the maximum engine by adding all of the engine’s status
engineStatus = ones(1,13);
maxEngine = sum(engineStatus);
% Initialize optimization problem
rel = zeros(1,13);
decision = optimproblem("ObjectiveSense","max");
engine = optimvar(‘engine’,13,’Type’,’integer’,’LowerBound’,0,’UpperBound’,1);
% Define days
for day = 1:5
% calculate maximum duration
maxDuration = 24*max(day);
% Initialize demand
for demand = 1:5
disp([‘Day: ‘,num2str(day),’ Demand: ‘,num2str(demand)])
switch demand
case 1
hourDemand = 3;
engineDemand = 9;
% Calculate reliability based on the sum between the hour
% demand & previous hour
if day == 1
timeCalc = hourDemand + totalHourWeekday(demand,:);
else
timeCalc = hourDemand + result(day-1,:);
end
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
% Set reliability as objective function
decision.Objective = relSys;
% Initialize constraint #1: engines should be equal to the
% demand
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
% Initialize constraint #2: engine’s hour should be less
% than max duration
if day == 1
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand,:) <= maxDuration;
else
decision.Constraints.durationConstraint = engine*result(day-1,:) <= maxDuration;
end
% Solve optimization problem
sol = solve(decision);
% Add the hour demand to the chosen engines
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 2
hourDemand = 2;
engineDemand = 6;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 3
hourDemand = 5;
engineDemand = 10;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 4
hourDemand = 8;
engineDemand = maxEngine-1;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 5
% Special case
for i = 1:13
if engineStatus(1,i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + 6;
end
end
end
% Add the total hour from all hour demands within the day
if day == 1
result(day,:) = sum(totalHourWeekday);
else
result(day,:) = result(day-1,:) + sum(totalHourWeekday);
end
end
end
The code ran well until it reached 4th day. When calculating 4th day’s 1st demand, it gave me the following error:
I guessed that, since it gave the infeasible linear constraints and bounds, an error was detected on line 57. That’s why I am wondering how to solve this problem or how to check for linear constraints and bounds infeasibility?
Furthermore, is there any way to reduce the code’s execution time?
Thank you in advanceHello,
I’m making a code where it creates an optimized engine usage schedule based on engine’s reliability using genetic algorithm solver. The following is the code that I’ve made:
tic
clear
% Create arrays to store total Hour & daily result
totalHourWeekday = zeros(5,13);
result = zeros(5,13);
% Calculate the maximum engine by adding all of the engine’s status
engineStatus = ones(1,13);
maxEngine = sum(engineStatus);
% Initialize optimization problem
rel = zeros(1,13);
decision = optimproblem("ObjectiveSense","max");
engine = optimvar(‘engine’,13,’Type’,’integer’,’LowerBound’,0,’UpperBound’,1);
% Define days
for day = 1:5
% calculate maximum duration
maxDuration = 24*max(day);
% Initialize demand
for demand = 1:5
disp([‘Day: ‘,num2str(day),’ Demand: ‘,num2str(demand)])
switch demand
case 1
hourDemand = 3;
engineDemand = 9;
% Calculate reliability based on the sum between the hour
% demand & previous hour
if day == 1
timeCalc = hourDemand + totalHourWeekday(demand,:);
else
timeCalc = hourDemand + result(day-1,:);
end
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
% Set reliability as objective function
decision.Objective = relSys;
% Initialize constraint #1: engines should be equal to the
% demand
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
% Initialize constraint #2: engine’s hour should be less
% than max duration
if day == 1
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand,:) <= maxDuration;
else
decision.Constraints.durationConstraint = engine*result(day-1,:) <= maxDuration;
end
% Solve optimization problem
sol = solve(decision);
% Add the hour demand to the chosen engines
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 2
hourDemand = 2;
engineDemand = 6;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 3
hourDemand = 5;
engineDemand = 10;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 4
hourDemand = 8;
engineDemand = maxEngine-1;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 5
% Special case
for i = 1:13
if engineStatus(1,i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + 6;
end
end
end
% Add the total hour from all hour demands within the day
if day == 1
result(day,:) = sum(totalHourWeekday);
else
result(day,:) = result(day-1,:) + sum(totalHourWeekday);
end
end
end
The code ran well until it reached 4th day. When calculating 4th day’s 1st demand, it gave me the following error:
I guessed that, since it gave the infeasible linear constraints and bounds, an error was detected on line 57. That’s why I am wondering how to solve this problem or how to check for linear constraints and bounds infeasibility?
Furthermore, is there any way to reduce the code’s execution time?
Thank you in advance Hello,
I’m making a code where it creates an optimized engine usage schedule based on engine’s reliability using genetic algorithm solver. The following is the code that I’ve made:
tic
clear
% Create arrays to store total Hour & daily result
totalHourWeekday = zeros(5,13);
result = zeros(5,13);
% Calculate the maximum engine by adding all of the engine’s status
engineStatus = ones(1,13);
maxEngine = sum(engineStatus);
% Initialize optimization problem
rel = zeros(1,13);
decision = optimproblem("ObjectiveSense","max");
engine = optimvar(‘engine’,13,’Type’,’integer’,’LowerBound’,0,’UpperBound’,1);
% Define days
for day = 1:5
% calculate maximum duration
maxDuration = 24*max(day);
% Initialize demand
for demand = 1:5
disp([‘Day: ‘,num2str(day),’ Demand: ‘,num2str(demand)])
switch demand
case 1
hourDemand = 3;
engineDemand = 9;
% Calculate reliability based on the sum between the hour
% demand & previous hour
if day == 1
timeCalc = hourDemand + totalHourWeekday(demand,:);
else
timeCalc = hourDemand + result(day-1,:);
end
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
% Set reliability as objective function
decision.Objective = relSys;
% Initialize constraint #1: engines should be equal to the
% demand
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
% Initialize constraint #2: engine’s hour should be less
% than max duration
if day == 1
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand,:) <= maxDuration;
else
decision.Constraints.durationConstraint = engine*result(day-1,:) <= maxDuration;
end
% Solve optimization problem
sol = solve(decision);
% Add the hour demand to the chosen engines
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 2
hourDemand = 2;
engineDemand = 6;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 3
hourDemand = 5;
engineDemand = 10;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 4
hourDemand = 8;
engineDemand = maxEngine-1;
timeCalc = hourDemand + totalHourWeekday(demand-1,:);
rel = exp(-(4.44*10^-5*timeCalc));
relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
decision.Objective = relSys;
engineMaintOrNot = sum(engine*engineStatus);
decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
sol = solve(decision);
for i = 1:13
if sol.engine(i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
end
end
case 5
% Special case
for i = 1:13
if engineStatus(1,i) == 1
totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + 6;
end
end
end
% Add the total hour from all hour demands within the day
if day == 1
result(day,:) = sum(totalHourWeekday);
else
result(day,:) = result(day-1,:) + sum(totalHourWeekday);
end
end
end
The code ran well until it reached 4th day. When calculating 4th day’s 1st demand, it gave me the following error:
I guessed that, since it gave the infeasible linear constraints and bounds, an error was detected on line 57. That’s why I am wondering how to solve this problem or how to check for linear constraints and bounds infeasibility?
Furthermore, is there any way to reduce the code’s execution time?
Thank you in advance genetic algorithm, optimization, error MATLAB Answers — New Questions
Bit stream (Delta Sigma ADC) to Decimation filter concerns
Hi, I have a second order delta sigma modulator that I want to process the output bitstream to get the analog output
Right now I have the output bitstream (varying from 0~3.3) exported from Cadence Virtuoso
And I passed the bitstream directly to a decimation filter made using filterDesigner setting the decimation factor the same as my OSR (32)
(input analog sine wave, output bitstream after delta sigma modulator, decimated sine wave on Matlab, FFT of the previous decimated plot)
However, the sine wave looks pretty distorted, I am guessing that’s due to the export sampling frequency (the simulation is sampled at much higher rate than my clock frequency) that Cadence Virtuoso does, but I am not certain.
The question is, is there any pre-filter process I need to do in order to decimate a digital data? I couldn’t find any example that process the bitstream directly to the decimation filter.
subplot(4,1,3);
xoutfilter = OUT{:,1};
youtfilter = Hm(OUT{:,2}); %Hm being the system object of FIRDecimator
plot(youtfilter)
xlabel(‘time’);
ylabel(‘voltage’);
Thanks a lot!Hi, I have a second order delta sigma modulator that I want to process the output bitstream to get the analog output
Right now I have the output bitstream (varying from 0~3.3) exported from Cadence Virtuoso
And I passed the bitstream directly to a decimation filter made using filterDesigner setting the decimation factor the same as my OSR (32)
(input analog sine wave, output bitstream after delta sigma modulator, decimated sine wave on Matlab, FFT of the previous decimated plot)
However, the sine wave looks pretty distorted, I am guessing that’s due to the export sampling frequency (the simulation is sampled at much higher rate than my clock frequency) that Cadence Virtuoso does, but I am not certain.
The question is, is there any pre-filter process I need to do in order to decimate a digital data? I couldn’t find any example that process the bitstream directly to the decimation filter.
subplot(4,1,3);
xoutfilter = OUT{:,1};
youtfilter = Hm(OUT{:,2}); %Hm being the system object of FIRDecimator
plot(youtfilter)
xlabel(‘time’);
ylabel(‘voltage’);
Thanks a lot! Hi, I have a second order delta sigma modulator that I want to process the output bitstream to get the analog output
Right now I have the output bitstream (varying from 0~3.3) exported from Cadence Virtuoso
And I passed the bitstream directly to a decimation filter made using filterDesigner setting the decimation factor the same as my OSR (32)
(input analog sine wave, output bitstream after delta sigma modulator, decimated sine wave on Matlab, FFT of the previous decimated plot)
However, the sine wave looks pretty distorted, I am guessing that’s due to the export sampling frequency (the simulation is sampled at much higher rate than my clock frequency) that Cadence Virtuoso does, but I am not certain.
The question is, is there any pre-filter process I need to do in order to decimate a digital data? I couldn’t find any example that process the bitstream directly to the decimation filter.
subplot(4,1,3);
xoutfilter = OUT{:,1};
youtfilter = Hm(OUT{:,2}); %Hm being the system object of FIRDecimator
plot(youtfilter)
xlabel(‘time’);
ylabel(‘voltage’);
Thanks a lot! adc, decimation, bitstream, fir, filter MATLAB Answers — New Questions
Different Sample time for RL environment
Hello Everyone,
I am trying to build a RL agent using DQN currently. My environment model is composed of Nonlinear equations which have faster dynamics around 1ms. I was wondering if it possible to have my RL agent to work at 10ms and let the environment run at 1ms. This will obviously mean that for 10 timesteps the environment will be feeding at a constant control action from the RL agent.
I know I can code this up if I create the agent and environment myself, but currently I’m taking advantage of MATLAB’s Reinforcement Learning Toolbox in MATLAB editor and I would save a lot of time if this is possible in the toolbox itself.
Thanks.Hello Everyone,
I am trying to build a RL agent using DQN currently. My environment model is composed of Nonlinear equations which have faster dynamics around 1ms. I was wondering if it possible to have my RL agent to work at 10ms and let the environment run at 1ms. This will obviously mean that for 10 timesteps the environment will be feeding at a constant control action from the RL agent.
I know I can code this up if I create the agent and environment myself, but currently I’m taking advantage of MATLAB’s Reinforcement Learning Toolbox in MATLAB editor and I would save a lot of time if this is possible in the toolbox itself.
Thanks. Hello Everyone,
I am trying to build a RL agent using DQN currently. My environment model is composed of Nonlinear equations which have faster dynamics around 1ms. I was wondering if it possible to have my RL agent to work at 10ms and let the environment run at 1ms. This will obviously mean that for 10 timesteps the environment will be feeding at a constant control action from the RL agent.
I know I can code this up if I create the agent and environment myself, but currently I’m taking advantage of MATLAB’s Reinforcement Learning Toolbox in MATLAB editor and I would save a lot of time if this is possible in the toolbox itself.
Thanks. rl environment, sample time, rl toolbox MATLAB Answers — New Questions
How to modify CSI-RS in wireless network simulator 5G?
Hi everyone,
I’m currently working on MATLAB’s wireless network simulator toolbox. I have two nodes: one gNB and one UE. I’m generating video traffic and using my own TDL channel. My problem is that when I add the Doppler effect to my channel, the automatic MCS selection is not optimal. I don’t understand why the channel estimation degrades with the Doppler effect. Therefore, I’m looking to adjust the CSI-RS parameters to see if this can improve the channel estimation and consequently the MCS selection.I don’t find any function who can permit me to adjust parameters and I can’t go depper in the code of wireless network simulator.Hi everyone,
I’m currently working on MATLAB’s wireless network simulator toolbox. I have two nodes: one gNB and one UE. I’m generating video traffic and using my own TDL channel. My problem is that when I add the Doppler effect to my channel, the automatic MCS selection is not optimal. I don’t understand why the channel estimation degrades with the Doppler effect. Therefore, I’m looking to adjust the CSI-RS parameters to see if this can improve the channel estimation and consequently the MCS selection.I don’t find any function who can permit me to adjust parameters and I can’t go depper in the code of wireless network simulator. Hi everyone,
I’m currently working on MATLAB’s wireless network simulator toolbox. I have two nodes: one gNB and one UE. I’m generating video traffic and using my own TDL channel. My problem is that when I add the Doppler effect to my channel, the automatic MCS selection is not optimal. I don’t understand why the channel estimation degrades with the Doppler effect. Therefore, I’m looking to adjust the CSI-RS parameters to see if this can improve the channel estimation and consequently the MCS selection.I don’t find any function who can permit me to adjust parameters and I can’t go depper in the code of wireless network simulator. csi-rs, simulation MATLAB Answers — New Questions
Separate graph arrows for two y-axes on the matlab
In the figure below, the circle arrows are used for two y-axes to indicate which axis the graph follows, how do you do this?In the figure below, the circle arrows are used for two y-axes to indicate which axis the graph follows, how do you do this? In the figure below, the circle arrows are used for two y-axes to indicate which axis the graph follows, how do you do this? y-axis MATLAB Answers — New Questions
How to set timezones for datetimes read from parquet files?
I just started playing around with parquet files for storing large tables of test data. All my data has timestamp columns which when saved has a timezone associated with it. When using parquetread(file) to get the whole table the datetimes are read in and display UTC timestamps for the display format with an empty timezone field. Currently I’m using the following to convert back to my local timezone for display purposes:
data = parquetread(file);
data.Timestamp.TimeZone = ‘UTC’;
data.Timestamp.TimeZone = ‘local’;
Is there a more concise way to do this?I just started playing around with parquet files for storing large tables of test data. All my data has timestamp columns which when saved has a timezone associated with it. When using parquetread(file) to get the whole table the datetimes are read in and display UTC timestamps for the display format with an empty timezone field. Currently I’m using the following to convert back to my local timezone for display purposes:
data = parquetread(file);
data.Timestamp.TimeZone = ‘UTC’;
data.Timestamp.TimeZone = ‘local’;
Is there a more concise way to do this? I just started playing around with parquet files for storing large tables of test data. All my data has timestamp columns which when saved has a timezone associated with it. When using parquetread(file) to get the whole table the datetimes are read in and display UTC timestamps for the display format with an empty timezone field. Currently I’m using the following to convert back to my local timezone for display purposes:
data = parquetread(file);
data.Timestamp.TimeZone = ‘UTC’;
data.Timestamp.TimeZone = ‘local’;
Is there a more concise way to do this? parquet datetime, datetime timezones MATLAB Answers — New Questions
Usage of the commande nnls(A,b,tol) in MATLAB new versions
Hello,
I want to solve the following non negative least squares problem : nnls(A,b,tol). Given a matrix A, a vector b and a tolerance tol, I want to fin x such ||A*x-b||^2 < tol* ||b||^2 with x>=0
Does someone have any idea how to do it in matlab ?
Regards !Hello,
I want to solve the following non negative least squares problem : nnls(A,b,tol). Given a matrix A, a vector b and a tolerance tol, I want to fin x such ||A*x-b||^2 < tol* ||b||^2 with x>=0
Does someone have any idea how to do it in matlab ?
Regards ! Hello,
I want to solve the following non negative least squares problem : nnls(A,b,tol). Given a matrix A, a vector b and a tolerance tol, I want to fin x such ||A*x-b||^2 < tol* ||b||^2 with x>=0
Does someone have any idea how to do it in matlab ?
Regards ! non negative least squares, matlab, lsqlin MATLAB Answers — New Questions
Why i’m getting this error
for i =1:N
fx =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaax1jjx’ num2str(i) ‘.txt’];
fy =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaay1jjy’ num2str(i) ‘.txt’];
fu =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaau11jju1’ num2str(i) ‘.txt’];
fv =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaau21jju2’ num2str(i) ‘.txt’];
x = textfile(fx);
y = textfile(fy);
z1 = textfile(fu);
z2 = textfile(fv);
U1 = griddata(x,y,z1,K1,K2);
U2 = griddata(x,y,z2,K1,K2);
zz1 = U1+zz1;
zz2 = U2+zz2;
i ;
end
Error using fscanf
Invalid file identifier. Use fopen to generate a valid file identifier.
Error in untitled>textfile (line 160)
value= fscanf(fileID,’%f’);
this is the error i’m gettingfor i =1:N
fx =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaax1jjx’ num2str(i) ‘.txt’];
fy =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaay1jjy’ num2str(i) ‘.txt’];
fu =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaau11jju1’ num2str(i) ‘.txt’];
fv =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaau21jju2’ num2str(i) ‘.txt’];
x = textfile(fx);
y = textfile(fy);
z1 = textfile(fu);
z2 = textfile(fv);
U1 = griddata(x,y,z1,K1,K2);
U2 = griddata(x,y,z2,K1,K2);
zz1 = U1+zz1;
zz2 = U2+zz2;
i ;
end
Error using fscanf
Invalid file identifier. Use fopen to generate a valid file identifier.
Error in untitled>textfile (line 160)
value= fscanf(fileID,’%f’);
this is the error i’m getting for i =1:N
fx =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaax1jjx’ num2str(i) ‘.txt’];
fy =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaay1jjy’ num2str(i) ‘.txt’];
fu =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaau11jju1’ num2str(i) ‘.txt’];
fv =[ ‘C:UsersSC23M048-BABARMALIKMDownloadsAbaqusHierarchical compositeaau21jju2’ num2str(i) ‘.txt’];
x = textfile(fx);
y = textfile(fy);
z1 = textfile(fu);
z2 = textfile(fv);
U1 = griddata(x,y,z1,K1,K2);
U2 = griddata(x,y,z2,K1,K2);
zz1 = U1+zz1;
zz2 = U2+zz2;
i ;
end
Error using fscanf
Invalid file identifier. Use fopen to generate a valid file identifier.
Error in untitled>textfile (line 160)
value= fscanf(fileID,’%f’);
this is the error i’m getting help MATLAB Answers — New Questions
