Tag Archives: matlab
How to check the convexity of the objective function?
Hi Everyone, I have the following objective function: Max 1-exp(- SNR_threshold/average SNR) and subject to power constraints,0>p< P_max. How can I check if it convex problem or not?Hi Everyone, I have the following objective function: Max 1-exp(- SNR_threshold/average SNR) and subject to power constraints,0>p< P_max. How can I check if it convex problem or not? Hi Everyone, I have the following objective function: Max 1-exp(- SNR_threshold/average SNR) and subject to power constraints,0>p< P_max. How can I check if it convex problem or not? wireless communication optimization MATLAB Answers — New Questions
Improving optimization results(Fmincon)
Hi guys. I’m trying to solve an optimization problem but the results I’m getting from fmincon() don’t have the accuracy that I’m looking for. I have tried to change the alghorithm and step tolerance but they didn’t affect the results. Is there any way to improve fmincon() results?? Here is my optimization problem code:
clear variables
clc
Objective=@MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
kL=500;
kH=500;
p0=[kL,kH];
options = optimoptions(‘fmincon’,’Display’,’iter’,’Algorithm’,’sqp-legacy’,’StepTolerance’,1e-11,"MaxFunctionEvaluations",2e3);
nlcon =[];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub, nlcon, options);
disp(k)
function MTE=MassTransferErrors_Closed_loop(p)
kL=p(1);
kH=p(2);
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
%Initial Condition
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T1=[230 230 230 180 200 230 230]; %T for different cases;
Kc_Total_gases=1;
tauI_Total_gases=1;
MTE_j=zeros(1,7);
Experiments = {[ 0.2985 0.6498 0.6147 0.43917 0.40398],[0.68662 1.6373 1.4260 1.4437 1.53169],[2.90493 5.68662 5.75704 2.65845 1.00352],[3.50352 11.3908 6.77817 3.46831 2.2007],[4.73592 10.8979 4.48944 3.01056 2.76408],[4.80634 9.45423 6.60211 4.03169 2.83451],[4.41901 10.4754 7.09507 4.13732 2.27113]};
for i=1:7
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization(t,y,Q0(i),T1(i),kL,kH,Kc_Total_gases,tauI_Total_gases),[0 18000],y0);
a = zeros(numel(t),1);
for q=1:numel(t)
a(q) = sum(y(q,1:7));
end
% Create plots for the gas mols in the reactor
% Total Gas mol in the reactor %[=mol]
%subplot(2,1,1)
%figure(i)
%plot(t,a,’r’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s]’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
% for t=100s
%subplot(2,1,2)
%t_new=1:50; %time vector for interpolation for plotting for the first 300 seconds
%aint =interp1(t,a,t_new,’pchip’); %interpolated state matrix for plotting
%plot(t_new,aint,’b’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s] (first 50s)’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
%hold off
% moles of products in gas and liquid at end
molGend=y(end,1:7);
molLend=y(end,8:14);
% product masses [g] in gas and liquid at end
massGend=molGend’.*moleWt;
massLend=molLend’.*moleWt;
%Total mass
TotalProduct = zeros(1,7);
for j=1:7
TotalProduct(j) = massGend(j) + massLend(j); %Sum of the liquid and gas phase products(g)
end
Experiment_i = cell2mat(Experiments(i)); %Converting Experiments set to matrix
MTE_i = ((TotalProduct(2)-Experiment_i(1))/Experiment_i(1))^2+((TotalProduct(3)-Experiment_i(2))/Experiment_i(2))^2+((TotalProduct(4)-Experiment_i(3))/Experiment_i(3))^2+((TotalProduct(5)-Experiment_i(4))/Experiment_i(4))^2+((TotalProduct(6)-Experiment_i(5))/Experiment_i(5))^2; %Defining an Mass Transfer Error relation
MTE_j(i) = MTE_i; %Defines a Mass Transfer Error Vector(1*7) that contains the error for each case
end
MTE = sum(MTE_j(1:7)); %Objective function(Sum of the all arrays in MTE_j Vector) what I need to minimize is each array that is on the MTE_j Vector but since I can’t return a Vector as an objective function I sum all the arrays as my objective function.
end
Also here is my code for the function that I have used in ode23s:
function S = Sec_model_fun_for_optimization(t,x,Q0,T1,kL,kH,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;…
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
Also the results are changing dramatically when I change initial values for kL and kH. I know it’s normal since fmincon() doesn’t compute global maximum/minimum but it’s just so weird to get different resluts whenever I change kL and kH values!Hi guys. I’m trying to solve an optimization problem but the results I’m getting from fmincon() don’t have the accuracy that I’m looking for. I have tried to change the alghorithm and step tolerance but they didn’t affect the results. Is there any way to improve fmincon() results?? Here is my optimization problem code:
clear variables
clc
Objective=@MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
kL=500;
kH=500;
p0=[kL,kH];
options = optimoptions(‘fmincon’,’Display’,’iter’,’Algorithm’,’sqp-legacy’,’StepTolerance’,1e-11,"MaxFunctionEvaluations",2e3);
nlcon =[];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub, nlcon, options);
disp(k)
function MTE=MassTransferErrors_Closed_loop(p)
kL=p(1);
kH=p(2);
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
%Initial Condition
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T1=[230 230 230 180 200 230 230]; %T for different cases;
Kc_Total_gases=1;
tauI_Total_gases=1;
MTE_j=zeros(1,7);
Experiments = {[ 0.2985 0.6498 0.6147 0.43917 0.40398],[0.68662 1.6373 1.4260 1.4437 1.53169],[2.90493 5.68662 5.75704 2.65845 1.00352],[3.50352 11.3908 6.77817 3.46831 2.2007],[4.73592 10.8979 4.48944 3.01056 2.76408],[4.80634 9.45423 6.60211 4.03169 2.83451],[4.41901 10.4754 7.09507 4.13732 2.27113]};
for i=1:7
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization(t,y,Q0(i),T1(i),kL,kH,Kc_Total_gases,tauI_Total_gases),[0 18000],y0);
a = zeros(numel(t),1);
for q=1:numel(t)
a(q) = sum(y(q,1:7));
end
% Create plots for the gas mols in the reactor
% Total Gas mol in the reactor %[=mol]
%subplot(2,1,1)
%figure(i)
%plot(t,a,’r’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s]’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
% for t=100s
%subplot(2,1,2)
%t_new=1:50; %time vector for interpolation for plotting for the first 300 seconds
%aint =interp1(t,a,t_new,’pchip’); %interpolated state matrix for plotting
%plot(t_new,aint,’b’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s] (first 50s)’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
%hold off
% moles of products in gas and liquid at end
molGend=y(end,1:7);
molLend=y(end,8:14);
% product masses [g] in gas and liquid at end
massGend=molGend’.*moleWt;
massLend=molLend’.*moleWt;
%Total mass
TotalProduct = zeros(1,7);
for j=1:7
TotalProduct(j) = massGend(j) + massLend(j); %Sum of the liquid and gas phase products(g)
end
Experiment_i = cell2mat(Experiments(i)); %Converting Experiments set to matrix
MTE_i = ((TotalProduct(2)-Experiment_i(1))/Experiment_i(1))^2+((TotalProduct(3)-Experiment_i(2))/Experiment_i(2))^2+((TotalProduct(4)-Experiment_i(3))/Experiment_i(3))^2+((TotalProduct(5)-Experiment_i(4))/Experiment_i(4))^2+((TotalProduct(6)-Experiment_i(5))/Experiment_i(5))^2; %Defining an Mass Transfer Error relation
MTE_j(i) = MTE_i; %Defines a Mass Transfer Error Vector(1*7) that contains the error for each case
end
MTE = sum(MTE_j(1:7)); %Objective function(Sum of the all arrays in MTE_j Vector) what I need to minimize is each array that is on the MTE_j Vector but since I can’t return a Vector as an objective function I sum all the arrays as my objective function.
end
Also here is my code for the function that I have used in ode23s:
function S = Sec_model_fun_for_optimization(t,x,Q0,T1,kL,kH,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;…
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
Also the results are changing dramatically when I change initial values for kL and kH. I know it’s normal since fmincon() doesn’t compute global maximum/minimum but it’s just so weird to get different resluts whenever I change kL and kH values! Hi guys. I’m trying to solve an optimization problem but the results I’m getting from fmincon() don’t have the accuracy that I’m looking for. I have tried to change the alghorithm and step tolerance but they didn’t affect the results. Is there any way to improve fmincon() results?? Here is my optimization problem code:
clear variables
clc
Objective=@MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
kL=500;
kH=500;
p0=[kL,kH];
options = optimoptions(‘fmincon’,’Display’,’iter’,’Algorithm’,’sqp-legacy’,’StepTolerance’,1e-11,"MaxFunctionEvaluations",2e3);
nlcon =[];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub, nlcon, options);
disp(k)
function MTE=MassTransferErrors_Closed_loop(p)
kL=p(1);
kH=p(2);
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
%Initial Condition
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T1=[230 230 230 180 200 230 230]; %T for different cases;
Kc_Total_gases=1;
tauI_Total_gases=1;
MTE_j=zeros(1,7);
Experiments = {[ 0.2985 0.6498 0.6147 0.43917 0.40398],[0.68662 1.6373 1.4260 1.4437 1.53169],[2.90493 5.68662 5.75704 2.65845 1.00352],[3.50352 11.3908 6.77817 3.46831 2.2007],[4.73592 10.8979 4.48944 3.01056 2.76408],[4.80634 9.45423 6.60211 4.03169 2.83451],[4.41901 10.4754 7.09507 4.13732 2.27113]};
for i=1:7
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization(t,y,Q0(i),T1(i),kL,kH,Kc_Total_gases,tauI_Total_gases),[0 18000],y0);
a = zeros(numel(t),1);
for q=1:numel(t)
a(q) = sum(y(q,1:7));
end
% Create plots for the gas mols in the reactor
% Total Gas mol in the reactor %[=mol]
%subplot(2,1,1)
%figure(i)
%plot(t,a,’r’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s]’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
% for t=100s
%subplot(2,1,2)
%t_new=1:50; %time vector for interpolation for plotting for the first 300 seconds
%aint =interp1(t,a,t_new,’pchip’); %interpolated state matrix for plotting
%plot(t_new,aint,’b’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s] (first 50s)’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
%hold off
% moles of products in gas and liquid at end
molGend=y(end,1:7);
molLend=y(end,8:14);
% product masses [g] in gas and liquid at end
massGend=molGend’.*moleWt;
massLend=molLend’.*moleWt;
%Total mass
TotalProduct = zeros(1,7);
for j=1:7
TotalProduct(j) = massGend(j) + massLend(j); %Sum of the liquid and gas phase products(g)
end
Experiment_i = cell2mat(Experiments(i)); %Converting Experiments set to matrix
MTE_i = ((TotalProduct(2)-Experiment_i(1))/Experiment_i(1))^2+((TotalProduct(3)-Experiment_i(2))/Experiment_i(2))^2+((TotalProduct(4)-Experiment_i(3))/Experiment_i(3))^2+((TotalProduct(5)-Experiment_i(4))/Experiment_i(4))^2+((TotalProduct(6)-Experiment_i(5))/Experiment_i(5))^2; %Defining an Mass Transfer Error relation
MTE_j(i) = MTE_i; %Defines a Mass Transfer Error Vector(1*7) that contains the error for each case
end
MTE = sum(MTE_j(1:7)); %Objective function(Sum of the all arrays in MTE_j Vector) what I need to minimize is each array that is on the MTE_j Vector but since I can’t return a Vector as an objective function I sum all the arrays as my objective function.
end
Also here is my code for the function that I have used in ode23s:
function S = Sec_model_fun_for_optimization(t,x,Q0,T1,kL,kH,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;…
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
Also the results are changing dramatically when I change initial values for kL and kH. I know it’s normal since fmincon() doesn’t compute global maximum/minimum but it’s just so weird to get different resluts whenever I change kL and kH values! optimization, ode MATLAB Answers — New Questions
in y axis i want 10^(-2000). can anyone plz help
set(gca, ‘YScale’, ‘log’)
ylim ([0 1e-308])
Instead of 1e-308 i want 1e-2000. how can i set the limit plzset(gca, ‘YScale’, ‘log’)
ylim ([0 1e-308])
Instead of 1e-308 i want 1e-2000. how can i set the limit plz set(gca, ‘YScale’, ‘log’)
ylim ([0 1e-308])
Instead of 1e-308 i want 1e-2000. how can i set the limit plz format y axis MATLAB Answers — New Questions
Thorlabs Linear Stage LTS150 Long Travel Stage with Matlab
Hi,
I have some trouble getting this stage to work in Matlab.
It works using the Throlabs App—Hence hardware is ok
I am able to read the serial number of the linear stage—-> The test computer does see the hardware.
But when trying to connect to the device—Matlab thinks the stage is not connected.
Anyone else may have found a solution for this issue?
Thanks,
——————- Code for future reference ——–
% Load .dll files associated with the Thorlabs motor
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.DeviceManagerCLI.dll’);
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.GenericMotorCLI.dll’);
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.IntegratedStepperMotorsCLI.dll’)
%NET.addAssembly(‘C:Program FilesThorlabsKinesisThorLabs.MotionControl.VerticalStageCLI.dll’)
import Thorlabs.MotionControl.DeviceManagerCLI.*
import Thorlabs.MotionControl.GenericMotorCLI.*
import Thorlabs.MotionControl.GenericMotorCLI.ControlParameters.*
import Thorlabs.MotionControl.GenericMotorCLI.AdvancedMotor.*
import Thorlabs.MotionControl.GenericMotorCLI.Settings.*
import Thorlabs.MotionControl.IntegratedStepperMotorsCLI.*
%import ThorLabs.MotionControl.VerticalStageCLI.*
% Generate device list
Thorlabs.MotionControl.DeviceManagerCLI.DeviceManagerCLI.BuildDeviceList() % Build device list
serialNumbersNet = Thorlabs.MotionControl.DeviceManagerCLI.DeviceManagerCLI.GetDeviceList() % Get device list
serialNumbers = cell(ToArray(serialNumbersNet)) % Convert serial numbers to cell array
serial_no = serialNumbers{1} % shows correct serial number (‘45357364’)
device = LongTravelStage.CreateLongTravelStage(serial_no);%The output of this line must be suppressed
device.Connect(serial_no);
—————-Hi,
I have some trouble getting this stage to work in Matlab.
It works using the Throlabs App—Hence hardware is ok
I am able to read the serial number of the linear stage—-> The test computer does see the hardware.
But when trying to connect to the device—Matlab thinks the stage is not connected.
Anyone else may have found a solution for this issue?
Thanks,
——————- Code for future reference ——–
% Load .dll files associated with the Thorlabs motor
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.DeviceManagerCLI.dll’);
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.GenericMotorCLI.dll’);
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.IntegratedStepperMotorsCLI.dll’)
%NET.addAssembly(‘C:Program FilesThorlabsKinesisThorLabs.MotionControl.VerticalStageCLI.dll’)
import Thorlabs.MotionControl.DeviceManagerCLI.*
import Thorlabs.MotionControl.GenericMotorCLI.*
import Thorlabs.MotionControl.GenericMotorCLI.ControlParameters.*
import Thorlabs.MotionControl.GenericMotorCLI.AdvancedMotor.*
import Thorlabs.MotionControl.GenericMotorCLI.Settings.*
import Thorlabs.MotionControl.IntegratedStepperMotorsCLI.*
%import ThorLabs.MotionControl.VerticalStageCLI.*
% Generate device list
Thorlabs.MotionControl.DeviceManagerCLI.DeviceManagerCLI.BuildDeviceList() % Build device list
serialNumbersNet = Thorlabs.MotionControl.DeviceManagerCLI.DeviceManagerCLI.GetDeviceList() % Get device list
serialNumbers = cell(ToArray(serialNumbersNet)) % Convert serial numbers to cell array
serial_no = serialNumbers{1} % shows correct serial number (‘45357364’)
device = LongTravelStage.CreateLongTravelStage(serial_no);%The output of this line must be suppressed
device.Connect(serial_no);
—————- Hi,
I have some trouble getting this stage to work in Matlab.
It works using the Throlabs App—Hence hardware is ok
I am able to read the serial number of the linear stage—-> The test computer does see the hardware.
But when trying to connect to the device—Matlab thinks the stage is not connected.
Anyone else may have found a solution for this issue?
Thanks,
——————- Code for future reference ——–
% Load .dll files associated with the Thorlabs motor
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.DeviceManagerCLI.dll’);
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.GenericMotorCLI.dll’);
NET.addAssembly(‘C:Program FilesThorlabsKinesisThorlabs.MotionControl.IntegratedStepperMotorsCLI.dll’)
%NET.addAssembly(‘C:Program FilesThorlabsKinesisThorLabs.MotionControl.VerticalStageCLI.dll’)
import Thorlabs.MotionControl.DeviceManagerCLI.*
import Thorlabs.MotionControl.GenericMotorCLI.*
import Thorlabs.MotionControl.GenericMotorCLI.ControlParameters.*
import Thorlabs.MotionControl.GenericMotorCLI.AdvancedMotor.*
import Thorlabs.MotionControl.GenericMotorCLI.Settings.*
import Thorlabs.MotionControl.IntegratedStepperMotorsCLI.*
%import ThorLabs.MotionControl.VerticalStageCLI.*
% Generate device list
Thorlabs.MotionControl.DeviceManagerCLI.DeviceManagerCLI.BuildDeviceList() % Build device list
serialNumbersNet = Thorlabs.MotionControl.DeviceManagerCLI.DeviceManagerCLI.GetDeviceList() % Get device list
serialNumbers = cell(ToArray(serialNumbersNet)) % Convert serial numbers to cell array
serial_no = serialNumbers{1} % shows correct serial number (‘45357364’)
device = LongTravelStage.CreateLongTravelStage(serial_no);%The output of this line must be suppressed
device.Connect(serial_no);
—————- thorlabs, linearstage, lts150, .net MATLAB Answers — New Questions
The sum of the 2 branches flow is different than pump flow
Hi there,
I am trying a very simple hydraulic circuit, consisting of a Fixed-Displacement Pump and two hydraulic branches, each one with a pressure relief valve. When I measure the flow after the pressure relief valves I found that the sum of both is different than the pump flow.
My conceptual idea is that the sum of both flows after relife valve should be equal (or closely) to the pump flow. Where is my mistake?Hi there,
I am trying a very simple hydraulic circuit, consisting of a Fixed-Displacement Pump and two hydraulic branches, each one with a pressure relief valve. When I measure the flow after the pressure relief valves I found that the sum of both is different than the pump flow.
My conceptual idea is that the sum of both flows after relife valve should be equal (or closely) to the pump flow. Where is my mistake? Hi there,
I am trying a very simple hydraulic circuit, consisting of a Fixed-Displacement Pump and two hydraulic branches, each one with a pressure relief valve. When I measure the flow after the pressure relief valves I found that the sum of both is different than the pump flow.
My conceptual idea is that the sum of both flows after relife valve should be equal (or closely) to the pump flow. Where is my mistake? simscape, simscape isothermal liquid, flow rate sensor il MATLAB Answers — New Questions
Skeletonisation until edges of a beam
Hello MATLAB community,
I am doing some image processing with MATLAB and some issues with my coding. I just like to warn you that I am very new at coding and MATLAB so I apologise in advance for my low level and I would be very glad to have some help as I have hitted a wall, and can’t find a solution to my problem.
Context: I have a video of beams, that move right to left over time. The base is fixed, only the beam moves. I converted the video to images, and my MATLAB program is going through the image file and treating every image in it. Here are two image examples:
and
I want to measure the following things:
a. The coordinates between the 2 extremities of the beam (length of the beam, without its base), let’s call them A and B.
b. The bending deformation E (L0-Lt/L0 *100), obtained by calculating the distance between A and B, called Lt.
c. the curvature of the beam (1/R), obtained by extracting the radius R of a circle fitting the curvature of the beam.
d. The angle between a vertical line passing through A, and the line AB.
What I have done so far:
My approach has been to transform my image into an rgbimage, then binaryImage, then have the complementary image, apply some modifications/corrections to the image, and then skeletonize it. And from then, I extract the coordinates of A and B, the distance between A&B (Lt), the radius of the beam R, and the angle between A&B (T).
My main issue is the skeletonisation. Because my beam is quite thick, it shortens up too much my beam, and in an inconcistent manner. So then my results are completly wrong. Here is an image of the different images and operation I have done and the result:
So as you can see, the length is shorter. I would like to have a skeleton that meets the edges of the beam to calculate the end points.
I have tried "bwskel(BW, ‘MinBranchLength’, 30)" and "bwmorph(BW, ‘thin’, inf)", and this: https://uk.mathworks.com/matlabcentral/fileexchange/11123-better-skeletonization. But the problem remains the same. I have tried regionpropos, but the major axis they return is too long, I have tried bwferet(), but the maxlength is in diagonal of the beam… I have running out of ideas.
Problem: So I guess my main problem is how can I get a skeletonisation that goes to the edges of the beam?
Here is my code:
for i = TrackingStart:TrackingEnd
FileRGB(:,:,i) = rgb2gray(imread(IMG)); % Convert to grayscale
croppedRGB = FileRGB(y3left:y3right, x3left:x3right, i);
binaryImage = imbinarize(croppedRGB, ‘adaptive’, ‘ForegroundPolarity’,’dark’,’Sensitivity’, 0.50);
out = nnz(~binaryImage);
while out <= 4300 % Change threshold if needed
for j = 1:50
sensitivity = 0.50 + j * 0.01;
binaryImage = imbinarize(croppedRGB, ‘adaptive’, ‘ForegroundPolarity’, ‘dark’, ‘Sensitivity’, sensitivity);
out = nnz(~binaryImage);
if out >= 4325
break; % Exit the loop if the condition is met
end
end
end
% Create a line Model
BW = imcomplement(binaryImage);
BW(y1left:y1right, x1left:x1right) = 1; % there is always sample at the junction area (between beam and base)
BW(y2left:y2right, x2left:x2right) = 0; % Always = 0 if no sample here
BW = bwmorph(BW, ‘close’, Inf);
BW = bwmorph(BW, ‘bridge’);
BW = bwareafilt(BW, 1);
s = regionprops(BW, ‘FilledImage’);
BW = s.FilledImage;
BW = bwskel(BW, ‘MinBranchLength’, 30);
endpoints = bwmorph(BW, ‘endpoints’);
[y_end, x_end] = find(endpoints == 1);
%Degree of bending deformation method
Lt = sqrt(power(x_end(1)-x_end(2),2)+power(y_end(1)-y_end(2),2));
if x_end(2) > x_end(1)
Lt = -Lt;
end
Lstore(i) = Lt;
%Curvature method
[row_dots_cir, col_dots_cir, val] = find(BW == 1);
[xc(i),yc(i),Rstore(i),a] = circfit(col_dots_cir,row_dots_cir);
%Angle method
slope_endpoints = (x_end(1) – x_end(2)) / (y_end(1) – y_end(2));
angle_radians = atan(slope_endpoints);
angle_degrees = rad2deg(angle_radians);
if x_end(2) > x_end(1)
angle_degrees = -angle_degrees;
end
Tstore(i) = angle_degrees;
i
endHello MATLAB community,
I am doing some image processing with MATLAB and some issues with my coding. I just like to warn you that I am very new at coding and MATLAB so I apologise in advance for my low level and I would be very glad to have some help as I have hitted a wall, and can’t find a solution to my problem.
Context: I have a video of beams, that move right to left over time. The base is fixed, only the beam moves. I converted the video to images, and my MATLAB program is going through the image file and treating every image in it. Here are two image examples:
and
I want to measure the following things:
a. The coordinates between the 2 extremities of the beam (length of the beam, without its base), let’s call them A and B.
b. The bending deformation E (L0-Lt/L0 *100), obtained by calculating the distance between A and B, called Lt.
c. the curvature of the beam (1/R), obtained by extracting the radius R of a circle fitting the curvature of the beam.
d. The angle between a vertical line passing through A, and the line AB.
What I have done so far:
My approach has been to transform my image into an rgbimage, then binaryImage, then have the complementary image, apply some modifications/corrections to the image, and then skeletonize it. And from then, I extract the coordinates of A and B, the distance between A&B (Lt), the radius of the beam R, and the angle between A&B (T).
My main issue is the skeletonisation. Because my beam is quite thick, it shortens up too much my beam, and in an inconcistent manner. So then my results are completly wrong. Here is an image of the different images and operation I have done and the result:
So as you can see, the length is shorter. I would like to have a skeleton that meets the edges of the beam to calculate the end points.
I have tried "bwskel(BW, ‘MinBranchLength’, 30)" and "bwmorph(BW, ‘thin’, inf)", and this: https://uk.mathworks.com/matlabcentral/fileexchange/11123-better-skeletonization. But the problem remains the same. I have tried regionpropos, but the major axis they return is too long, I have tried bwferet(), but the maxlength is in diagonal of the beam… I have running out of ideas.
Problem: So I guess my main problem is how can I get a skeletonisation that goes to the edges of the beam?
Here is my code:
for i = TrackingStart:TrackingEnd
FileRGB(:,:,i) = rgb2gray(imread(IMG)); % Convert to grayscale
croppedRGB = FileRGB(y3left:y3right, x3left:x3right, i);
binaryImage = imbinarize(croppedRGB, ‘adaptive’, ‘ForegroundPolarity’,’dark’,’Sensitivity’, 0.50);
out = nnz(~binaryImage);
while out <= 4300 % Change threshold if needed
for j = 1:50
sensitivity = 0.50 + j * 0.01;
binaryImage = imbinarize(croppedRGB, ‘adaptive’, ‘ForegroundPolarity’, ‘dark’, ‘Sensitivity’, sensitivity);
out = nnz(~binaryImage);
if out >= 4325
break; % Exit the loop if the condition is met
end
end
end
% Create a line Model
BW = imcomplement(binaryImage);
BW(y1left:y1right, x1left:x1right) = 1; % there is always sample at the junction area (between beam and base)
BW(y2left:y2right, x2left:x2right) = 0; % Always = 0 if no sample here
BW = bwmorph(BW, ‘close’, Inf);
BW = bwmorph(BW, ‘bridge’);
BW = bwareafilt(BW, 1);
s = regionprops(BW, ‘FilledImage’);
BW = s.FilledImage;
BW = bwskel(BW, ‘MinBranchLength’, 30);
endpoints = bwmorph(BW, ‘endpoints’);
[y_end, x_end] = find(endpoints == 1);
%Degree of bending deformation method
Lt = sqrt(power(x_end(1)-x_end(2),2)+power(y_end(1)-y_end(2),2));
if x_end(2) > x_end(1)
Lt = -Lt;
end
Lstore(i) = Lt;
%Curvature method
[row_dots_cir, col_dots_cir, val] = find(BW == 1);
[xc(i),yc(i),Rstore(i),a] = circfit(col_dots_cir,row_dots_cir);
%Angle method
slope_endpoints = (x_end(1) – x_end(2)) / (y_end(1) – y_end(2));
angle_radians = atan(slope_endpoints);
angle_degrees = rad2deg(angle_radians);
if x_end(2) > x_end(1)
angle_degrees = -angle_degrees;
end
Tstore(i) = angle_degrees;
i
end Hello MATLAB community,
I am doing some image processing with MATLAB and some issues with my coding. I just like to warn you that I am very new at coding and MATLAB so I apologise in advance for my low level and I would be very glad to have some help as I have hitted a wall, and can’t find a solution to my problem.
Context: I have a video of beams, that move right to left over time. The base is fixed, only the beam moves. I converted the video to images, and my MATLAB program is going through the image file and treating every image in it. Here are two image examples:
and
I want to measure the following things:
a. The coordinates between the 2 extremities of the beam (length of the beam, without its base), let’s call them A and B.
b. The bending deformation E (L0-Lt/L0 *100), obtained by calculating the distance between A and B, called Lt.
c. the curvature of the beam (1/R), obtained by extracting the radius R of a circle fitting the curvature of the beam.
d. The angle between a vertical line passing through A, and the line AB.
What I have done so far:
My approach has been to transform my image into an rgbimage, then binaryImage, then have the complementary image, apply some modifications/corrections to the image, and then skeletonize it. And from then, I extract the coordinates of A and B, the distance between A&B (Lt), the radius of the beam R, and the angle between A&B (T).
My main issue is the skeletonisation. Because my beam is quite thick, it shortens up too much my beam, and in an inconcistent manner. So then my results are completly wrong. Here is an image of the different images and operation I have done and the result:
So as you can see, the length is shorter. I would like to have a skeleton that meets the edges of the beam to calculate the end points.
I have tried "bwskel(BW, ‘MinBranchLength’, 30)" and "bwmorph(BW, ‘thin’, inf)", and this: https://uk.mathworks.com/matlabcentral/fileexchange/11123-better-skeletonization. But the problem remains the same. I have tried regionpropos, but the major axis they return is too long, I have tried bwferet(), but the maxlength is in diagonal of the beam… I have running out of ideas.
Problem: So I guess my main problem is how can I get a skeletonisation that goes to the edges of the beam?
Here is my code:
for i = TrackingStart:TrackingEnd
FileRGB(:,:,i) = rgb2gray(imread(IMG)); % Convert to grayscale
croppedRGB = FileRGB(y3left:y3right, x3left:x3right, i);
binaryImage = imbinarize(croppedRGB, ‘adaptive’, ‘ForegroundPolarity’,’dark’,’Sensitivity’, 0.50);
out = nnz(~binaryImage);
while out <= 4300 % Change threshold if needed
for j = 1:50
sensitivity = 0.50 + j * 0.01;
binaryImage = imbinarize(croppedRGB, ‘adaptive’, ‘ForegroundPolarity’, ‘dark’, ‘Sensitivity’, sensitivity);
out = nnz(~binaryImage);
if out >= 4325
break; % Exit the loop if the condition is met
end
end
end
% Create a line Model
BW = imcomplement(binaryImage);
BW(y1left:y1right, x1left:x1right) = 1; % there is always sample at the junction area (between beam and base)
BW(y2left:y2right, x2left:x2right) = 0; % Always = 0 if no sample here
BW = bwmorph(BW, ‘close’, Inf);
BW = bwmorph(BW, ‘bridge’);
BW = bwareafilt(BW, 1);
s = regionprops(BW, ‘FilledImage’);
BW = s.FilledImage;
BW = bwskel(BW, ‘MinBranchLength’, 30);
endpoints = bwmorph(BW, ‘endpoints’);
[y_end, x_end] = find(endpoints == 1);
%Degree of bending deformation method
Lt = sqrt(power(x_end(1)-x_end(2),2)+power(y_end(1)-y_end(2),2));
if x_end(2) > x_end(1)
Lt = -Lt;
end
Lstore(i) = Lt;
%Curvature method
[row_dots_cir, col_dots_cir, val] = find(BW == 1);
[xc(i),yc(i),Rstore(i),a] = circfit(col_dots_cir,row_dots_cir);
%Angle method
slope_endpoints = (x_end(1) – x_end(2)) / (y_end(1) – y_end(2));
angle_radians = atan(slope_endpoints);
angle_degrees = rad2deg(angle_radians);
if x_end(2) > x_end(1)
angle_degrees = -angle_degrees;
end
Tstore(i) = angle_degrees;
i
end transferred, image analysis, bending, deformation MATLAB Answers — New Questions
Adaptive pipelining design cannot insert the required number of registers in a feedback loop with integral modules
In my design I need to use the following design approach:
But after generating HDL code using HDLCoder, the report suggests that:
"Adaptive pipelining is not supported if the block is in a feedback loop. Number of registers required: 2; Number of registers inserted: 0.(from MATLABr21b)” and "Unable to insert required number of pipeline registers because the block, gain, is in a feedback loop and there are not enough latency budget at the output of the block. Number of registers required: 2; number of registers inserted: 0. Consider increasing the latency budget by adding more design delays in the feedback loop or using clock-rate pipelining.(from MATLABr23b)"
At the moment I don’t know how to go about modifying it, I’d like to know if it’s feasible for me to do so, and also how this needs to be implemented for the clock rate pipeline, I don’t have the slightest clue, and I didn’t understand what I should do from reading the official instructions. I hope to get the guidance of professionals, thank you very much.
If you need to provide more instructions, I am very happy to add!
Here is a screenshot of the reportIn my design I need to use the following design approach:
But after generating HDL code using HDLCoder, the report suggests that:
"Adaptive pipelining is not supported if the block is in a feedback loop. Number of registers required: 2; Number of registers inserted: 0.(from MATLABr21b)” and "Unable to insert required number of pipeline registers because the block, gain, is in a feedback loop and there are not enough latency budget at the output of the block. Number of registers required: 2; number of registers inserted: 0. Consider increasing the latency budget by adding more design delays in the feedback loop or using clock-rate pipelining.(from MATLABr23b)"
At the moment I don’t know how to go about modifying it, I’d like to know if it’s feasible for me to do so, and also how this needs to be implemented for the clock rate pipeline, I don’t have the slightest clue, and I didn’t understand what I should do from reading the official instructions. I hope to get the guidance of professionals, thank you very much.
If you need to provide more instructions, I am very happy to add!
Here is a screenshot of the report In my design I need to use the following design approach:
But after generating HDL code using HDLCoder, the report suggests that:
"Adaptive pipelining is not supported if the block is in a feedback loop. Number of registers required: 2; Number of registers inserted: 0.(from MATLABr21b)” and "Unable to insert required number of pipeline registers because the block, gain, is in a feedback loop and there are not enough latency budget at the output of the block. Number of registers required: 2; number of registers inserted: 0. Consider increasing the latency budget by adding more design delays in the feedback loop or using clock-rate pipelining.(from MATLABr23b)"
At the moment I don’t know how to go about modifying it, I’d like to know if it’s feasible for me to do so, and also how this needs to be implemented for the clock rate pipeline, I don’t have the slightest clue, and I didn’t understand what I should do from reading the official instructions. I hope to get the guidance of professionals, thank you very much.
If you need to provide more instructions, I am very happy to add!
Here is a screenshot of the report adaptive pipelining design, registers MATLAB Answers — New Questions
how to set parameters for ble connection in matlab
Hello,
I want to implement a ble connection between a peripheral board and Matlab, which should acts as a central device. I’m using ble function to establish the connection and then write and read function and relative characteristics. Everything is working but I want to set for my connection the parameters such as connection interval, phy, slave latency etc and I don’t know how to do it. Is it present any example which shows how to set configuration parameters?Hello,
I want to implement a ble connection between a peripheral board and Matlab, which should acts as a central device. I’m using ble function to establish the connection and then write and read function and relative characteristics. Everything is working but I want to set for my connection the parameters such as connection interval, phy, slave latency etc and I don’t know how to do it. Is it present any example which shows how to set configuration parameters? Hello,
I want to implement a ble connection between a peripheral board and Matlab, which should acts as a central device. I’m using ble function to establish the connection and then write and read function and relative characteristics. Everything is working but I want to set for my connection the parameters such as connection interval, phy, slave latency etc and I don’t know how to do it. Is it present any example which shows how to set configuration parameters? bluetooth, ble, bluetooth low energy, matlab, toolbox MATLAB Answers — New Questions
How can I get the average color of a ROI (mask) ?
Hello !
I’m doing a project for school and I have a picture that I divided in several ROI.
I used roipoly in matlab to select the ROI. Now, I would like to get the average color inside the area of the mask and replace each pixels but the average color. Here is my code :
averageColors = zeros(sum(mask(:)), 3);
for i = 1:3
channel = image(:,:,i);
averageColors(:,i) = accumarray(find(mask(:)), channel(mask(:)), [], @mean);
end
result = image;
for i = 1:3
result(:,:,i) = accumarray(find(mask(:)), averageColors(:,i), [], @mean);
end
I get an error on ligne 4 :
Unable to perform assignment because the size of the left side is 30683-by-1 and the size of the right side is 421368-by-1.
I don’t understand… find(mask(:)) and channel(mask(:)) have the same size…
Thank you for reading !Hello !
I’m doing a project for school and I have a picture that I divided in several ROI.
I used roipoly in matlab to select the ROI. Now, I would like to get the average color inside the area of the mask and replace each pixels but the average color. Here is my code :
averageColors = zeros(sum(mask(:)), 3);
for i = 1:3
channel = image(:,:,i);
averageColors(:,i) = accumarray(find(mask(:)), channel(mask(:)), [], @mean);
end
result = image;
for i = 1:3
result(:,:,i) = accumarray(find(mask(:)), averageColors(:,i), [], @mean);
end
I get an error on ligne 4 :
Unable to perform assignment because the size of the left side is 30683-by-1 and the size of the right side is 421368-by-1.
I don’t understand… find(mask(:)) and channel(mask(:)) have the same size…
Thank you for reading ! Hello !
I’m doing a project for school and I have a picture that I divided in several ROI.
I used roipoly in matlab to select the ROI. Now, I would like to get the average color inside the area of the mask and replace each pixels but the average color. Here is my code :
averageColors = zeros(sum(mask(:)), 3);
for i = 1:3
channel = image(:,:,i);
averageColors(:,i) = accumarray(find(mask(:)), channel(mask(:)), [], @mean);
end
result = image;
for i = 1:3
result(:,:,i) = accumarray(find(mask(:)), averageColors(:,i), [], @mean);
end
I get an error on ligne 4 :
Unable to perform assignment because the size of the left side is 30683-by-1 and the size of the right side is 421368-by-1.
I don’t understand… find(mask(:)) and channel(mask(:)) have the same size…
Thank you for reading ! roipoly, average, color, pixels, mask, roi, region of interest, image processing MATLAB Answers — New Questions
Help with simulation hanging
the simulation reach 100% and hange what is the solution and I try to change the solver but the same problem?the simulation reach 100% and hange what is the solution and I try to change the solver but the same problem? the simulation reach 100% and hange what is the solution and I try to change the solver but the same problem? simulation MATLAB Answers — New Questions
How do i solve bilevel problem?
I have bilevel optimization, I write upper level with intlinprog and write lower level in linprog and now i want to solve it together because they depend to each other on some variable, so how do i solve them together?I have bilevel optimization, I write upper level with intlinprog and write lower level in linprog and now i want to solve it together because they depend to each other on some variable, so how do i solve them together? I have bilevel optimization, I write upper level with intlinprog and write lower level in linprog and now i want to solve it together because they depend to each other on some variable, so how do i solve them together? optimization MATLAB Answers — New Questions
Oxford Battery degradtion dataset
Hi All,
I am struggling to load all the files in .csv format can some one help me in this. I have loaded the data of each cell (numbered 1-8). But unable to get the data which every cell has that is:
cyc0000 1×1 struct 1×1 struct
C1ch 1×1 struct 1×1 struct
t 3510×1 double 3510×1 double
v 3510×1 double 3510×1 double
q 3510×1 double 3510×1 double
T 3510×1 double 3510×1 double
C1dc 1×1 struct 1×1 struct
OCVch 1×1 struct 1×1 struct
OCVdc 1×1 struct 1×1 struct
cyc0100 1×1 struct 1×1 struct
cyc0200 1×1 struct 1×1 struct
cyc0300 1×1 struct 1×1 struct
cyc0400 1×1 struct 1×1 struct….
and so on.Hi All,
I am struggling to load all the files in .csv format can some one help me in this. I have loaded the data of each cell (numbered 1-8). But unable to get the data which every cell has that is:
cyc0000 1×1 struct 1×1 struct
C1ch 1×1 struct 1×1 struct
t 3510×1 double 3510×1 double
v 3510×1 double 3510×1 double
q 3510×1 double 3510×1 double
T 3510×1 double 3510×1 double
C1dc 1×1 struct 1×1 struct
OCVch 1×1 struct 1×1 struct
OCVdc 1×1 struct 1×1 struct
cyc0100 1×1 struct 1×1 struct
cyc0200 1×1 struct 1×1 struct
cyc0300 1×1 struct 1×1 struct
cyc0400 1×1 struct 1×1 struct….
and so on. Hi All,
I am struggling to load all the files in .csv format can some one help me in this. I have loaded the data of each cell (numbered 1-8). But unable to get the data which every cell has that is:
cyc0000 1×1 struct 1×1 struct
C1ch 1×1 struct 1×1 struct
t 3510×1 double 3510×1 double
v 3510×1 double 3510×1 double
q 3510×1 double 3510×1 double
T 3510×1 double 3510×1 double
C1dc 1×1 struct 1×1 struct
OCVch 1×1 struct 1×1 struct
OCVdc 1×1 struct 1×1 struct
cyc0100 1×1 struct 1×1 struct
cyc0200 1×1 struct 1×1 struct
cyc0300 1×1 struct 1×1 struct
cyc0400 1×1 struct 1×1 struct….
and so on. matlab MATLAB Answers — New Questions
Help with wireless communication
Hello, everyone
I need help on using matlab for wireless communication, for example implementation of alamouti schemeHello, everyone
I need help on using matlab for wireless communication, for example implementation of alamouti scheme Hello, everyone
I need help on using matlab for wireless communication, for example implementation of alamouti scheme wireless MATLAB Answers — New Questions
DC-DC multiphase buck converter
Hi everyone,
may you help in finding the schematic of a 4-phase DC-DC buck converter?
Thank you so much.
Regards,
FrancescoHi everyone,
may you help in finding the schematic of a 4-phase DC-DC buck converter?
Thank you so much.
Regards,
Francesco Hi everyone,
may you help in finding the schematic of a 4-phase DC-DC buck converter?
Thank you so much.
Regards,
Francesco multiphase buck converter MATLAB Answers — New Questions
Rte_Type.h in model reference
I have a autosar classic component that has a model reference is using datatypes defined in the component. Code generation is successfull but, I can not get the metrics report.
Error messages as follows
In component code generation report, section static code metrics: (3 Lines error)
component.h line 65 identifier "XYZ" is undefined, XYZ is a enum class
component.c line 32 identifier "ABC" is undefined, ABC is a enum value
model reference.h Line 34 expected a declaration
In model reference code generation report, section static code metrics: (1 Line error)
model reference.h Line 26 cannot open source file "Rte_Type.h"I have a autosar classic component that has a model reference is using datatypes defined in the component. Code generation is successfull but, I can not get the metrics report.
Error messages as follows
In component code generation report, section static code metrics: (3 Lines error)
component.h line 65 identifier "XYZ" is undefined, XYZ is a enum class
component.c line 32 identifier "ABC" is undefined, ABC is a enum value
model reference.h Line 34 expected a declaration
In model reference code generation report, section static code metrics: (1 Line error)
model reference.h Line 26 cannot open source file "Rte_Type.h" I have a autosar classic component that has a model reference is using datatypes defined in the component. Code generation is successfull but, I can not get the metrics report.
Error messages as follows
In component code generation report, section static code metrics: (3 Lines error)
component.h line 65 identifier "XYZ" is undefined, XYZ is a enum class
component.c line 32 identifier "ABC" is undefined, ABC is a enum value
model reference.h Line 34 expected a declaration
In model reference code generation report, section static code metrics: (1 Line error)
model reference.h Line 26 cannot open source file "Rte_Type.h" autosar simulink MATLAB Answers — New Questions
Optimize loop to reduce run time
My code works well but it takes a lot of time to run since the vectors are larger in size. I know it can be optimized by using vector instead of for loop but I am unable to do that properly. Please help
prob_det_glrt_signoise_fix_pfa_2d = zeros(num_trials,length(noise_sigma));
prob_det_glrt_signoise_fix_pfaa_2d = zeros(1,length(noise_sigma));
inv_sqrt = 1/sqrt(2);
for jj = 1:length(noise_sigma)
for ii = 1:num_trials
complex_awgn = noise_sigma(jj)*(randn(numrangebins,numdoppbins)+1i*randn(numrangebins,numdoppbins))*inv_sqrt;
s_matrix_2D = u_matrix_2D + complex_awgn; % received signal with noise
pow_val_2D = (abs (s_matrix_2D(bin_rng_2d,bin_dopp_2d)) )^2;
prob_det_glrt_signoise_fix_pfa_2d(ii,jj) = pow_val_2D > threshold(jj); % Hypothesis H1
end
prob_det_glrt_signoise_fix_pfaa_2d(jj) = mean(prob_det_glrt_signoise_fix_pfa_2d(:,jj));
endMy code works well but it takes a lot of time to run since the vectors are larger in size. I know it can be optimized by using vector instead of for loop but I am unable to do that properly. Please help
prob_det_glrt_signoise_fix_pfa_2d = zeros(num_trials,length(noise_sigma));
prob_det_glrt_signoise_fix_pfaa_2d = zeros(1,length(noise_sigma));
inv_sqrt = 1/sqrt(2);
for jj = 1:length(noise_sigma)
for ii = 1:num_trials
complex_awgn = noise_sigma(jj)*(randn(numrangebins,numdoppbins)+1i*randn(numrangebins,numdoppbins))*inv_sqrt;
s_matrix_2D = u_matrix_2D + complex_awgn; % received signal with noise
pow_val_2D = (abs (s_matrix_2D(bin_rng_2d,bin_dopp_2d)) )^2;
prob_det_glrt_signoise_fix_pfa_2d(ii,jj) = pow_val_2D > threshold(jj); % Hypothesis H1
end
prob_det_glrt_signoise_fix_pfaa_2d(jj) = mean(prob_det_glrt_signoise_fix_pfa_2d(:,jj));
end My code works well but it takes a lot of time to run since the vectors are larger in size. I know it can be optimized by using vector instead of for loop but I am unable to do that properly. Please help
prob_det_glrt_signoise_fix_pfa_2d = zeros(num_trials,length(noise_sigma));
prob_det_glrt_signoise_fix_pfaa_2d = zeros(1,length(noise_sigma));
inv_sqrt = 1/sqrt(2);
for jj = 1:length(noise_sigma)
for ii = 1:num_trials
complex_awgn = noise_sigma(jj)*(randn(numrangebins,numdoppbins)+1i*randn(numrangebins,numdoppbins))*inv_sqrt;
s_matrix_2D = u_matrix_2D + complex_awgn; % received signal with noise
pow_val_2D = (abs (s_matrix_2D(bin_rng_2d,bin_dopp_2d)) )^2;
prob_det_glrt_signoise_fix_pfa_2d(ii,jj) = pow_val_2D > threshold(jj); % Hypothesis H1
end
prob_det_glrt_signoise_fix_pfaa_2d(jj) = mean(prob_det_glrt_signoise_fix_pfa_2d(:,jj));
end for loop, optimization MATLAB Answers — New Questions
Adjusting the spacing between the bars of different groups within a category using boxchart.
I am using boxchart to plot data as in the example: openExample(‘graphics/UseCombinationOfGroupingVariablesExample’)
How do I reduce the spacing between the groups (blue and red) so that they touch each other for every category? This is required for me as I have six groups. The spacing between the goups makes it difficult to differentiate between boxes of different categories. I want the boxes of each category to be clubbed together. Any suggestions would be highly appreciated.I am using boxchart to plot data as in the example: openExample(‘graphics/UseCombinationOfGroupingVariablesExample’)
How do I reduce the spacing between the groups (blue and red) so that they touch each other for every category? This is required for me as I have six groups. The spacing between the goups makes it difficult to differentiate between boxes of different categories. I want the boxes of each category to be clubbed together. Any suggestions would be highly appreciated. I am using boxchart to plot data as in the example: openExample(‘graphics/UseCombinationOfGroupingVariablesExample’)
How do I reduce the spacing between the groups (blue and red) so that they touch each other for every category? This is required for me as I have six groups. The spacing between the goups makes it difficult to differentiate between boxes of different categories. I want the boxes of each category to be clubbed together. Any suggestions would be highly appreciated. matlab, boxchart, plot, box plot MATLAB Answers — New Questions
Stateflow parameter defined in dspace data dictionary
Is it possible to link a stateflow parameter to a variable defined in the dspace data dictionary? If, so how do I go about it?
Thanks, AlanIs it possible to link a stateflow parameter to a variable defined in the dspace data dictionary? If, so how do I go about it?
Thanks, Alan Is it possible to link a stateflow parameter to a variable defined in the dspace data dictionary? If, so how do I go about it?
Thanks, Alan state flow, data dictionary, parameter MATLAB Answers — New Questions
how to get value of a Variable from dSPACE Data Dictionary using MATLAB Command ?
I have one DD opened and I would like to check each variable’s Value if it is defined or not using only the MATLAB commands. So, Does anyone have any idea about this ?I have one DD opened and I would like to check each variable’s Value if it is defined or not using only the MATLAB commands. So, Does anyone have any idea about this ? I have one DD opened and I would like to check each variable’s Value if it is defined or not using only the MATLAB commands. So, Does anyone have any idea about this ? dd, dspacedd, dsdd, dsddman, matlab, simulink, targetlink MATLAB Answers — New Questions
Graphene conductance and dielectric permittivity
Hello,
I’m trying to reproduce the conductance and the dielectric permittivity of graphene, as described in this link: https://support.lumerical.com/hc/en-us/articles/360042244874-Graphene-surface-conductivity-material-model .
To do this, I wrote the following code in MATLAB:
wls_interval=0.3:0.01:5; % wavelength in microns
c=3*10^8; % speed of light
f=(c./(wls_interval*10^-6)); %Hz
w=2*pi.*f; %rad/s
f=f.*10^-12; % frequency in THz
E=1240./(wls_interval*10^3); % energy expressed as 1240/(wls in nanometers)
htag=6.582*10^-16; % reduced Planck constant in eV*s
hconst=htag*2*pi; % Planck constant in eV*s
%w=2*pi*c./(wls_interval.*10^-6);
epsr=1; % relative permittivity of vacuum
e=-1.602*10^-19; % Electron’s charge in C
kb=8.617*10^-5; % Boltzmann constant in eV/K
T=300; % Temperature in °K
mu_c=0.2; % Chemical potential in eV
eps0=8.854*10^-12; % Vacuum dielectric permittivity (F/m)
n_sheets=1; %number of sheets by which the graphene layer is made of
t_sheet=0.336*10^-9; % thickness of a graphene sheet in nm
tg=n_sheets*t_sheet; % total thickness of the graphene layer
Gamma=0.41*10^-3; % Scattering rate in eV
parfor i=1:length(E)
%—– Interband contribution numerical —–
FD=@(x) 1./(exp((x-mu_c)./(kb*T))+1); %Fermi-Dirac distribution
den=@(x) (E(i)+1i*2*Gamma).^2-4*(x./htag).^2;
fun=@(x) (FD(-x)-FD(x))./den(x);
eta(i)=1i*e^2*(E(i)+1i*2*Gamma)/(pi*htag^2);
int_term(i)=integral(fun,0,E(i));
sigma_inter(i)=eta(i)*int_term(i);
%—————————————————————————-
%—– Intraband contribution analytical —–
alpha(i)=1i*e^2*kb*T/(pi*htag^2*(E(i)+2*1i*Gamma));
beta=mu_c/(kb*T);
sigma_intra(i)=alpha(i)*(beta+2*log(exp(-beta)+1));
sigma_tot(i)=sigma_intra(i)+sigma_inter(i);
eps_graph_par(i)=eps0*(epsr(i)+1i*sigma_tot(i)/(eps0*E(i)*tg));
refind(i)=sqrt(eps_graph_par(i));
end
The results I’m supposed to get are plotted in panel "a" and "b", while the ones I actually get are in panel "d-g". Even though you may notice that the real interband conductivity may resemble to some extent the real total conductivty I’m supposed to get, this is not the case for the imaginary part, that is totally different. Moreover, the order of magnitude is completely out of scale (see panel "c,d"), therefore, the contribution of the interband part is minimal also in the visible range, while it should be significant.
Can you please help me in finding the error?
Thank you so much!Hello,
I’m trying to reproduce the conductance and the dielectric permittivity of graphene, as described in this link: https://support.lumerical.com/hc/en-us/articles/360042244874-Graphene-surface-conductivity-material-model .
To do this, I wrote the following code in MATLAB:
wls_interval=0.3:0.01:5; % wavelength in microns
c=3*10^8; % speed of light
f=(c./(wls_interval*10^-6)); %Hz
w=2*pi.*f; %rad/s
f=f.*10^-12; % frequency in THz
E=1240./(wls_interval*10^3); % energy expressed as 1240/(wls in nanometers)
htag=6.582*10^-16; % reduced Planck constant in eV*s
hconst=htag*2*pi; % Planck constant in eV*s
%w=2*pi*c./(wls_interval.*10^-6);
epsr=1; % relative permittivity of vacuum
e=-1.602*10^-19; % Electron’s charge in C
kb=8.617*10^-5; % Boltzmann constant in eV/K
T=300; % Temperature in °K
mu_c=0.2; % Chemical potential in eV
eps0=8.854*10^-12; % Vacuum dielectric permittivity (F/m)
n_sheets=1; %number of sheets by which the graphene layer is made of
t_sheet=0.336*10^-9; % thickness of a graphene sheet in nm
tg=n_sheets*t_sheet; % total thickness of the graphene layer
Gamma=0.41*10^-3; % Scattering rate in eV
parfor i=1:length(E)
%—– Interband contribution numerical —–
FD=@(x) 1./(exp((x-mu_c)./(kb*T))+1); %Fermi-Dirac distribution
den=@(x) (E(i)+1i*2*Gamma).^2-4*(x./htag).^2;
fun=@(x) (FD(-x)-FD(x))./den(x);
eta(i)=1i*e^2*(E(i)+1i*2*Gamma)/(pi*htag^2);
int_term(i)=integral(fun,0,E(i));
sigma_inter(i)=eta(i)*int_term(i);
%—————————————————————————-
%—– Intraband contribution analytical —–
alpha(i)=1i*e^2*kb*T/(pi*htag^2*(E(i)+2*1i*Gamma));
beta=mu_c/(kb*T);
sigma_intra(i)=alpha(i)*(beta+2*log(exp(-beta)+1));
sigma_tot(i)=sigma_intra(i)+sigma_inter(i);
eps_graph_par(i)=eps0*(epsr(i)+1i*sigma_tot(i)/(eps0*E(i)*tg));
refind(i)=sqrt(eps_graph_par(i));
end
The results I’m supposed to get are plotted in panel "a" and "b", while the ones I actually get are in panel "d-g". Even though you may notice that the real interband conductivity may resemble to some extent the real total conductivty I’m supposed to get, this is not the case for the imaginary part, that is totally different. Moreover, the order of magnitude is completely out of scale (see panel "c,d"), therefore, the contribution of the interband part is minimal also in the visible range, while it should be significant.
Can you please help me in finding the error?
Thank you so much! Hello,
I’m trying to reproduce the conductance and the dielectric permittivity of graphene, as described in this link: https://support.lumerical.com/hc/en-us/articles/360042244874-Graphene-surface-conductivity-material-model .
To do this, I wrote the following code in MATLAB:
wls_interval=0.3:0.01:5; % wavelength in microns
c=3*10^8; % speed of light
f=(c./(wls_interval*10^-6)); %Hz
w=2*pi.*f; %rad/s
f=f.*10^-12; % frequency in THz
E=1240./(wls_interval*10^3); % energy expressed as 1240/(wls in nanometers)
htag=6.582*10^-16; % reduced Planck constant in eV*s
hconst=htag*2*pi; % Planck constant in eV*s
%w=2*pi*c./(wls_interval.*10^-6);
epsr=1; % relative permittivity of vacuum
e=-1.602*10^-19; % Electron’s charge in C
kb=8.617*10^-5; % Boltzmann constant in eV/K
T=300; % Temperature in °K
mu_c=0.2; % Chemical potential in eV
eps0=8.854*10^-12; % Vacuum dielectric permittivity (F/m)
n_sheets=1; %number of sheets by which the graphene layer is made of
t_sheet=0.336*10^-9; % thickness of a graphene sheet in nm
tg=n_sheets*t_sheet; % total thickness of the graphene layer
Gamma=0.41*10^-3; % Scattering rate in eV
parfor i=1:length(E)
%—– Interband contribution numerical —–
FD=@(x) 1./(exp((x-mu_c)./(kb*T))+1); %Fermi-Dirac distribution
den=@(x) (E(i)+1i*2*Gamma).^2-4*(x./htag).^2;
fun=@(x) (FD(-x)-FD(x))./den(x);
eta(i)=1i*e^2*(E(i)+1i*2*Gamma)/(pi*htag^2);
int_term(i)=integral(fun,0,E(i));
sigma_inter(i)=eta(i)*int_term(i);
%—————————————————————————-
%—– Intraband contribution analytical —–
alpha(i)=1i*e^2*kb*T/(pi*htag^2*(E(i)+2*1i*Gamma));
beta=mu_c/(kb*T);
sigma_intra(i)=alpha(i)*(beta+2*log(exp(-beta)+1));
sigma_tot(i)=sigma_intra(i)+sigma_inter(i);
eps_graph_par(i)=eps0*(epsr(i)+1i*sigma_tot(i)/(eps0*E(i)*tg));
refind(i)=sqrt(eps_graph_par(i));
end
The results I’m supposed to get are plotted in panel "a" and "b", while the ones I actually get are in panel "d-g". Even though you may notice that the real interband conductivity may resemble to some extent the real total conductivty I’m supposed to get, this is not the case for the imaginary part, that is totally different. Moreover, the order of magnitude is completely out of scale (see panel "c,d"), therefore, the contribution of the interband part is minimal also in the visible range, while it should be significant.
Can you please help me in finding the error?
Thank you so much! graphene, conductivity, kubo formula MATLAB Answers — New Questions