fplot: Attempted to access x(2); index out of bounds because numel(x)=1.

Consider plot function.

B = arrayfun(@(x)Lobo(x(1), x(2), x(3), x(4), x(5), x(6)), x(1), 'un',0);
fplot(B, [0.01 0.2]);

Error:

Attempted to access x(2); index out of bounds because numel(x)=1.

[Birdman]

Lobo function?

[Devdatt Thengdi]

function [qout] = Lobo(d0, L, Tg, nb, uo, Ef)
clc
%d0 = 0.01;
%L = 10;
%Tg = 1950;
%nb = 1;
%uo = 10;

%input data
Ti = 525.05; %Fluid in degK
T0 = 659.25; %Fluid/gas out degK
cpL = 0.32; %Cp fluid KJ/kgdegK
flowL = 255.5; %flowrate fluid kg/s
THD = flowL*cpL*(T0 - Ti);%Total Heat Duty  
%Ef = 0.9320; %Fuel efficiency
qfuel = THD/Ef;%Heat released by fuel Kw
LHV = 41.05; %Fuel value KJ/kg
mfuel = qfuel/LHV;
%pi = 22/7;

%absorptivity
alpha = 0.9; 

%Exchange factor
x = pi; %ctc
B = (x^2) - 4;
C = sqrt(B);
F = (C - x + 2*sinh(2/x))/(2*pi);
F
disp(F);

%Heat transfer coefficient
k = 0.25*(10^-3); %Thermal conductivity
%rof = 1320; %density
mu = 130; %Viscosity Pa-s
A0 = (pi*(d0^2))/4;
h0 = ((d0/k)*(0.023*((flowL*d0/mu*A0)^0.8)*((cpL*mu/k))^0.3));
d0
disp(d0);
h0
disp(h0);

%Tube wall temperature
Tref = 41 + 273.15; %Ambient
Tw = Tref + (T0 - Ti)./L; 
Tw
disp(Tw)

%Radiant heat flux
stef = 5.67*10^-8; %Stefan-Boltzmann constant
Tg4 = Tg.*Tg.*Tg.*Tg;
Tw4 = Tw.*Tw.*Tw.*Tw;
dT4 = [Tg4(:,1) - Tw4(:,1)];
dT = Tg(:,1) - Tw(:,1);
qr = stef*alpha*F*(dT4) + h0*(dT);
qr
disp(qr);

%set conditions
%if qr >= (45000) && qr <= (32000)
%qr = qr;
%else
%qr = randi ([32000, 45000], 1);
%qr = qr;
%end

%Air
ea = 0.15; %excess air
cpair = 0.171; %Kcal/kgdegK
Tairin = 114;
atfr = 1.344; %in notebook
mair = mfuel*atfr*(1+ea);

%Flame calculation
cpsteam = 4.1855;
cpco2 = 0.918;
cpn2 = 1.044;
xco2 = 0.054;
xsteam = 0.034; %Total should be 8.8
xn2 = 0.66;
cpgas = xsteam*cpsteam + xco2*cpco2 + xn2*cpn2;
kgas = 0.25*10^-3; %Thermal conductivity
mgas = mfuel + mair;
Lf = 0.0042*(qfuel^0.478); %flame length
Lf
disp(Lf);
tf = sqrt(pi*uo*Lf*cpgas/mgas*kgas); %flame size
tf
disp(tf);
db = 1.5*tf; %Burner tile diameter
dbc = db(:)*(1./(sin(180./nb) + (pi./nb))); %Burner circle diameter
dbc
disp(dbc);

%Cold plane area
cbt = (qfuel*1.055/(4*10^6)) + 1.5;
cbt
disp(cbt);
dtc = cbt + dbc; %tube circle diameter
dtc
disp(dtc);
ntubes = pi*dtc/x;
ntubes
disp(ntubes);
n1 = (ntubes(:)*L);
Ar = x.*n1;
Ar
disp(Ar);

%Heat leaving
Qr = qr * Ar; %total radiant heat absorbed
qair = mair*cpair*(Tairin - Tref);
qloss = 0.05*qfuel;
qout = (qfuel + qair - (Qr + qloss));

if qout >= 0
   qout = qout; 
else
qout = 0;
end
qout
disp(qout);
end

[KSSV]

What is your x value? You have not shown that. I suspect your X is not a vector.....it is a scalar, so the error popping out.

[Devdatt Thengdi]

I am trying to call the 5 variables given here:

Lobo(d0, L, Tg, nb, uo)

But, as you said, they are being called as scalars (x(i)). How do I call them as vectors? I short, How do I plot the value of Lobo wrt to all 5 variables.

[Devdatt Thengdi]

What I did initially was that I extracted the values from the workspace (optimresults.x) and plot them using excel sheet. But, I need to know how do I that here?

[Walter Roberson]

Plotting Lobo with respect to all 5 variables would require a 6 dimensional plot. In my experience, 5 dimensions is the maximum you can encode and still have a semi-understandable ranking.

[Devdatt Thengdi]

Also, before plotting I am optimizing the same function. Now, in this line:

qfuel = THD/Ef;%Heat released by fuel Kw

I've observed that optimized values of x(1) to x(5) vary with Ef. How can I plot the 'optimized values' vs Ef (or vice versa)?

Answer 1

[Walter Roberson]

You wrote

B = arrayfun(@(x)Lobo(x(1), x(2), x(3), x(4), x(5), x(6)), x(1), 'un',0);

Break this down:

   Array_to_operate_over = x(1);
   Function_to_call = @(x)Lobo(x(1), x(2), x(3), x(4), x(5), x(6));
   B = arrayfun(Function_to_call, Array_to_operate_over, 'un', 0);

This requires that a vector or array x existed before the B = statement that you coded. The first element of that vector or array is extracted because of your x(1); and a function is executed for each member of that scalar because of the arrayfun(). The function you are calling tries to extract 6 elements from the scalar and pass the 6 of them into Lobo, but a scalar only has 1 elements so you get an index out of range error.

arrayfun() always proceeds one element at a time, so even if you had passed an array in, instead of the scalar x(1), then what would reach the Function_to_call would always be a scalar.

If somehow it all worked, then because of the 'un', 0 you would be returning a cell array of values. But you then proceed to try to fplot() that cell array of values, which is not possible.

[Devdatt Thengdi]

@Walter Roberson Okay. let's put it this way then. I've used the same function to optimize x(1) to x(5). What if wanted to plot Ef from Lobo:

qfuel = THD/Ef;%Heat released by fuel Kw

I've observed that optimized values of x(1) to x(5) vary with Ef. How can I plot the 'optimized values' vs Ef ?

Since you can not create a function inside a function and I tried to loop the function when I call it in the Solution file but, it's not defined in the function file.

[Walter Roberson]

How could you have used a function with six inputs to "optimize" 5 of the inputs? What does it mean to "optimize" an input?

"Since you can not create a function inside a function"

I have no idea why you say that.

function outer
   t = 0;
   function out = inner(in)
     out = factorial(in);
   end
   for K = 1 : 10
     t = t + inner(K);
   end
end

This defines function "inner" inside function "outer".

https://www.mathworks.com/help/matlab/matlab_prog/nested-functions.html

Nested functions have been available since MATLAB 7.0 (R14) (not R2014*, I mean R14, released in 2004.)

[Devdatt Thengdi]

Oh sorry! I removed one of the variables. Let me put it this way; I minimized the function Lobo (which is now a 5 variable function). I determined the values of x which bring out the minimal values of function Lobo. I need to plot Ef(0.01:0.99) vs x(1 to 5). How do I do that? I'll try nested functions now.

[Walter Roberson]

Plotting Lobo with respect to all 5 variables would require a 6 dimensional plot. In my experience, 5 dimensions is the maximum you can encode and still have a semi-understandable ranking.

With what you have posted so far, Ef and Lobo appear to be independent.

[Devdatt Thengdi]

:

function [qout] = Lobo(d0, L, Tg, nb, uo)
clc
%d0 = 0.01;
%L = 10;
%Tg = 1950;
%nb = 1;
%uo = 10;

%input data
Ti = 525.05; %Fluid in degK
T0 = 659.25; %Fluid/gas out degK
cpL = 0.32; %Cp fluid KJ/kgdegK
flowL = 255.5; %flowrate fluid kg/s
THD = flowL*cpL*(T0 - Ti);%Total Heat Duty  
Ef = 0.9320; %Fuel efficiency
qfuel = THD/Ef;%Heat released by fuel Kw
LHV = 41.05; %Fuel value KJ/kg
mfuel = qfuel/LHV;
%pi = 22/7;

%absorptivity
alpha = 0.9; 

%Exchange factor
x = pi; %ctc
B = (x^2) - 4;
C = sqrt(B);
F = (C - x + 2*sinh(2/x))/(2*pi);
F
disp(F);

%Heat transfer coefficient
k = 0.25*(10^-3); %Thermal conductivity
%rof = 1320; %density
mu = 130; %Viscosity Pa-s
A0 = (pi*(d0^2))/4;
h0 = ((d0/k)*(0.023*((flowL*d0/mu*A0)^0.8)*((cpL*mu/k))^0.3));
d0
disp(d0);
h0
disp(h0);

%Tube wall temperature
Tref = 41 + 273.15; %Ambient
Tw = Tref + (T0 - Ti)./L; 
Tw
disp(Tw)

%Radiant heat flux
stef = 5.67*10^-8; %Stefan-Boltzmann constant
Tg4 = Tg.*Tg.*Tg.*Tg;
Tw4 = Tw.*Tw.*Tw.*Tw;
dT4 = [Tg4(:,1) - Tw4(:,1)];
dT = Tg(:,1) - Tw(:,1);
qr = stef*alpha*F*(dT4) + h0*(dT);
qr
disp(qr);

%set conditions
%if qr >= (45000) && qr <= (32000)
%qr = qr;
%else
%qr = randi ([32000, 45000], 1);
%qr = qr;
%end

%Air
ea = 0.15; %excess air
cpair = 0.171; %Kcal/kgdegK
Tairin = 114;
atfr = 1.344; %in notebook
mair = mfuel*atfr*(1+ea);

%Flame calculation
cpsteam = 4.1855;
cpco2 = 0.918;
cpn2 = 1.044;
xco2 = 0.054;
xsteam = 0.034; %Total should be 8.8
xn2 = 0.66;
cpgas = xsteam*cpsteam + xco2*cpco2 + xn2*cpn2;
kgas = 0.25*10^-3; %Thermal conductivity
mgas = mfuel + mair;
Lf = 0.0042*(qfuel^0.478); %flame length
Lf
disp(Lf);
tf = sqrt(pi*uo*Lf*cpgas/mgas*kgas); %flame size
tf
disp(tf);
db = 1.5*tf; %Burner tile diameter
dbc = db(:)*(1./(sin(180./nb) + (pi./nb))); %Burner circle diameter
dbc
disp(dbc);

%Cold plane area
cbt = (qfuel*1.055/(4*10^6)) + 1.5;
cbt
disp(cbt);
dtc = cbt + dbc; %tube circle diameter
dtc
disp(dtc);
ntubes = pi*dtc/x;
ntubes
disp(ntubes);
n1 = (ntubes(:)*L);
ctc = pi*d0/2;
Ar = ctc.*n1;
Ar
disp(Ar);

%Heat leaving
Qr = qr * Ar; %total radiant heat absorbed
qair = mair*cpair*(Tairin - Tref);
qloss = 0.05*qfuel;
qout = (qfuel + qair - (Qr + qloss));

if qout >= 0
   qout = qout; 
else
qout = 0;
end
qout
disp(qout);
end

Solution:

clc
x0 = [0.01, 1, 10, 8, 1];
%options = optimset('PlotFcns',@optimplotfval);
fitnessfcn = @(x)Lobo(x(1), x(2), x(3), x(4), x(5));
nonlcon1 = @(x)nonlcon(x(1), x(2), x(3), x(4), x(5));
lb = [0.01, 0.1, 1, 8, 1];
ub = [1, 18.3, 1950 ,20, 50];
[X, fval] = fmincon(fitnessfcn, x0, [], [], [], [], lb, ub, nonlcon1);
X
disp(X);

This is the line where Ef is used:

   qfuel = THD/Ef;%Heat released by fuel Kw

Now, for each of value of Ef(0.01:0.99) we get a different set of [X] in the solution function (it might oscillate too.) I need to plot Ef vs [X]. Get it?

[Torsten]

I don't know how you want to plot Ef against 5 solution components, but here is the code:

Solution:

Ef = 0.01:0.01:0.99;
Xmat = zeros(numel(Ef),5);
for i=1:numel(Ef)
  Ef_actual = Ef(i);
  x0 = [0.01, 1, 10, 8, 1];
  %options = optimset('PlotFcns',@optimplotfval);
  fitnessfcn = @(x)Lobo(x(1), x(2), x(3), x(4), x(5), Ef_actual);
  nonlcon1 = @(x)nonlcon(x(1), x(2), x(3), x(4), x(5), Ef_actual);
  lb = [0.01, 0.1, 1, 8, 1];
  ub = [1, 18.3, 1950 ,20, 50];
  [X, fval] = fmincon(fitnessfcn, x0, [], [], [], [], lb, ub, nonlcon1);
  Xmat(i,:) = X(:);
end

Lobo:
function [qout] = Lobo(d0, L, Tg, nb, uo, Ef)
clc
%d0 = 0.01;
%L = 10;
%Tg = 1950;
%nb = 1;
%uo = 10;

%input data
Ti = 525.05; %Fluid in degK
T0 = 659.25; %Fluid/gas out degK
cpL = 0.32; %Cp fluid KJ/kgdegK
flowL = 255.5; %flowrate fluid kg/s
THD = flowL*cpL*(T0 - Ti);%Total Heat Duty  
%Ef = 0.9320; %Fuel efficiency
qfuel = THD/Ef;%Heat released by fuel Kw
LHV = 41.05; %Fuel value KJ/kg
mfuel = qfuel/LHV;
%pi = 22/7;

%absorptivity
alpha = 0.9; 

%Exchange factor
x = pi; %ctc
B = (x^2) - 4;
C = sqrt(B);
F = (C - x + 2*sinh(2/x))/(2*pi);
F
disp(F);

%Heat transfer coefficient
k = 0.25*(10^-3); %Thermal conductivity
%rof = 1320; %density
mu = 130; %Viscosity Pa-s
A0 = (pi*(d0^2))/4;
h0 = ((d0/k)*(0.023*((flowL*d0/mu*A0)^0.8)*((cpL*mu/k))^0.3));
d0
disp(d0);
h0
disp(h0);

%Tube wall temperature
Tref = 41 + 273.15; %Ambient
Tw = Tref + (T0 - Ti)./L; 
Tw
disp(Tw)

%Radiant heat flux
stef = 5.67*10^-8; %Stefan-Boltzmann constant
Tg4 = Tg.*Tg.*Tg.*Tg;
Tw4 = Tw.*Tw.*Tw.*Tw;
dT4 = [Tg4(:,1) - Tw4(:,1)];
dT = Tg(:,1) - Tw(:,1);
qr = stef*alpha*F*(dT4) + h0*(dT);
qr
disp(qr);

%set conditions
%if qr >= (45000) && qr <= (32000)
%qr = qr;
%else
%qr = randi ([32000, 45000], 1);
%qr = qr;
%end

%Air
ea = 0.15; %excess air
cpair = 0.171; %Kcal/kgdegK
Tairin = 114;
atfr = 1.344; %in notebook
mair = mfuel*atfr*(1+ea);

%Flame calculation
cpsteam = 4.1855;
cpco2 = 0.918;
cpn2 = 1.044;
xco2 = 0.054;
xsteam = 0.034; %Total should be 8.8
xn2 = 0.66;
cpgas = xsteam*cpsteam + xco2*cpco2 + xn2*cpn2;
kgas = 0.25*10^-3; %Thermal conductivity
mgas = mfuel + mair;
Lf = 0.0042*(qfuel^0.478); %flame length
Lf
disp(Lf);
tf = sqrt(pi*uo*Lf*cpgas/mgas*kgas); %flame size
tf
disp(tf);
db = 1.5*tf; %Burner tile diameter
dbc = db(:)*(1./(sin(180./nb) + (pi./nb))); %Burner circle diameter
dbc
disp(dbc);

%Cold plane area
cbt = (qfuel*1.055/(4*10^6)) + 1.5;
cbt
disp(cbt);
dtc = cbt + dbc; %tube circle diameter
dtc
disp(dtc);
ntubes = pi*dtc/x;
ntubes
disp(ntubes);
n1 = (ntubes(:)*L);
ctc = pi*d0/2;
Ar = ctc.*n1;
Ar
disp(Ar);

%Heat leaving
Qr = qr * Ar; %total radiant heat absorbed
qair = mair*cpair*(Tairin - Tref);
qloss = 0.05*qfuel;
qout = (qfuel + qair - (Qr + qloss));

if qout >= 0
   qout = qout; 
else
qout = 0;
end
qout
disp(qout);
end

nonlcon:

function [c ceq] = nonlcon(d0, L, Tg, nb, uo, Ef)

...

[Devdatt Thengdi]

So, if want to plot it against the first variable only,

plot(Ef, Xmat(i:1));

right?

Thanks man.

[Torsten]

plot(Ef, Xmat(:,1));