Two variable Optimization: Variables is defined in a tertiary equation. Design of a Furnace.

function qr = Sol(~)
%Assume Q and d0
Q = 7000000; %Kcal/m^2degC
d0 = 0.033; %Tube diameter m
THD = 3773.09; %Total heat duty Gcal/h
%Target fluid
Ti = 251.9; %Fluid in degC
T0 = 386.1; %Fluid/gas out degC
cpl = 0.674; %Cp fluid Kcal/kgdegC
LHV = 10.515; %Fuel value Kcal/kg
flowL = 12.9; %flowrate fluid
%Air
ea = 0.15; %excess air
cpair = 0.171; %Kcal/kgdegK
cpstm = 0.5125; %Kcal/kgdegK
Tairin = 114;
Tstmin = 141;
Tref = 41; %Ambient
%Fuel
Hea = 111.5; %Heat released by fuel Gcal/h 
Ef = 0.9320; %Fuel efficiency
atfr = 10.344; %in notebook
stfr = 0.5; %literature
%Dimensions
Tg = 869; %degC
Tw = 100+0.5*(Ti+T0); %Tube wall temperature
ctc = 0.3048; %ctc in m
Lh = 19.236; %Exposed tube length- horizontal m
Lv = 18.431; %Exposed tube length- vertical m 
pi = 22/7;
dstk = 4.435;
lstk = 42.405;
h = 2; %approximation
sig = 0.173*10^-08; 
%Radiant section Calcuations
qfuel = THD/Ef;
mfuel = qfuel/LHV;
mreqair = mfuel*atfr;
mair = mreqair*(1+ea);
qair = mair*cpair*(Tairin - Tref);
mstm = mfuel*stfr;
qstm = mstm*cpstm*(Tstmin - Tref);
A = pi*d0*(Lh+Lv);
Tavg = (Ti + T0)/2; %Assumption
qr = h*A*(Tw - Tavg);
qout = 0.1*Q;
qlibr = qfuel+qair+qstm-qr-qout;
nt = qlibr/A*Q;
Acp = (Lh + Lv)*nt*ctc;
Ar = Acp;
r = ctc/d0;
al = 1 - (0.0277+0.0927*(r - 1))*r - 1;
Vrad = 4*A; %4cuft/sqft of radiant area (neglecting convection section)
%Vstk = pi*dstk*lstk;
%VF = Vrad + Vstk;
Di = 10.000; %Furnace Diameter
Lm = Di; %Assuming L/D >=2
P = 0.288 - 0.229*ea + 0.09*ea^2;
pl = P*Lm;
z = (Tg+460)/1000;
a = 0.47916 - 0.19847*z + 0.022569*z*z;
b = 0.047029 + 0.0699*z - 0.01528*z*z;
c = 0.000803 - 0.00726*z + 0.00159*z*z;
e = a + b*(pl) + c*pl*pl;
Aw = 2*pi*Di*Di - Acp;
%Exchange factor
y = Aw/(al*Ar);
u = 0.00064 + 0.0591*y + 0.00101*y*y;
v = 1.0256 + 0.4908*y - 0.058*y*y;
w = -0.144 - 0.552*y + 0.040*y*y;
F = u + e*v + w*e*e;
Qr = ((Tg)^4 - (Tw)^4) + h*A*(Tg - Tw)/(sig*al*Ar*F);
qr = Qr/(sig*al*Ar*F);  
disp(qr);
disp(d0);
disp(Q);
if abs(qr - Q) > 1000    
return
end
%disp(F);
%disp(e);
%disp(pl);
%disp(a);disp(b);disp(c);
%disp(y);
%disp(u);
%disp(v);
%disp(w);
%disp(Tw);
%disp(pl);
%disp(y);
%disp(e);
  end

Here, we are assuming a value of Q to calculate it again later using obtained sub-variables and compare if it is valid enough. If its not we can either assume a new value of Q or d0. Thus, I Need to optimize qr over Q and d0. Both of them are not defined directly in the objective function but, in secondary or tertiary equations. Do I need to include them in the objective function? If yes, how? if no, please suggest me a suitable optimzation solvers for the same. Are we getting 2 objective functions here?

[Walter Roberson]

I am not clear what your optimization is?

You construct qr in terms of Q and d0. You test whether the difference between qr and Q is at least 1000 and return if so. Are you trying to maximize the difference between qr and Q? If so then my exploration suggests that there are an infinite number of locations (a curve) at which there is a denominator of 0, leading to infinite differences. This is especially notable for Q < 1/10.

[Devdatt Thengdi]

@walter Sorry. That's a typo. I need to minimize the difference. The process is like this:

1. We are assuming a value of qr itself, renaming it to Q.
2. We carry out the whole procedure based on the the assumed value.
3. Once done, we calculate new value of qr based on the values of terms derived from assumed qr (Q) using the model equation given.
4. Then we compare the assumed value to the calculated value and aim for minimum difference between them.
5. Now, if the condition is satisfied it's fine and dandy but otherwise, I would have to either, a - assume a new value of qr (Q) and repeat or, b - change the value of d0.
6. This gives us (as I can see) 2 objective functions i.e. model equation and the difference between Q and qr and two variables Q and d0.
7. It is easy to construct the model equation in the form of d0 (simple substitution of values) but, I am held up while solving for Q.

Get it?

Answer 1

[Walter Roberson]

If you have the symbolic toolbox, use

syms d0 nonnegative
syms Q real

and skip the initialization of those by numeric values and run the rest of your code up to the point

if abs(qr - Q) > 1000

This gives you

qr = (166227300197011511522099200000000*d0^2*((1967832025192261675700914749440000000*d0^3)/(944520154538697205399*Q*(353187/(12500000*d0) - 13/200)*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736)*((301779270272994367636807587565564375*d0^2*((2667*Q*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736))/(27500*d0) + 4400/7))/(53670583394087919043943355748712448*Q*(353187/(12500000*d0) - 13/200)*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736)) - (422947675392560621728619988744831396484375*d0^4*((2667*Q*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736))/(27500*d0) + 4400/7)^2)/(20451048334681640592886886342576975118336*Q^2*(353187/(12500000*d0) - 13/200)^2*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736)^2) + 953597232042240918401963381071191/1980704062856608439838598758400000)) + 539446471200))/(8500681390848274848591*Q*(353187/(12500000*d0) - 13/200)*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736)*((301779270272994367636807587565564375*d0^2*((2667*Q*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736))/(27500*d0) + 4400/7))/(53670583394087919043943355748712448*Q*(353187/(12500000*d0) - 13/200)*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736)) - (422947675392560621728619988744831396484375*d0^4*((2667*Q*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736))/(27500*d0) + 4400/7)^2)/(20451048334681640592886886342576975118336*Q^2*(353187/(12500000*d0) - 13/200)^2*(Q/10 + (118382*d0)/5 - 4884955544470871/68719476736)^2) + 953597232042240918401963381071191/1980704062856608439838598758400000))

For any given numeric d0 such as 0.033

temp = subs(qr - Q, d0, 0.033);
soltemp = solve(temp, 'returnconditions', true);
sol = solve(soltemp.conditions);

This will give 6 roots out of an equation of degree 8 (the other two are imaginary and so were ruled out by our declaring that Q had to be real).

Then for a target initial Q such as 700000

[~, minidx] = min(abs(sol - 700000));
bestQ  = double(sol(minidx));

These are exact solutions to within numeric round-off.

You can also

soltemp = solve(qr-Q,Q, 'returnconditions', true);
sol = solve(soltemp.conditions,soltemp.parameters, 'returnconditions', true);

and sol.x will then be the solutions (with each one only valid if the corresponding sol.conditions entry is true)

What this tells us is that it is not useful to optimize d0 and Q simultaneously, because no matter what d0 you choose (within reason), you can find an exact Q.

Likewise if you substitute a particular numeric Q into the qr-Q equation, then you can solve for d0. This gives a polynomial of degree 12. For Q = 700000 then four of the roots of the polynomial are imaginary and 2 are negative, giving 6 positive real roots; the smallest non-negative of which is 0.0458456070692645991358438592521 .

If you have a constraint on d0 (notice this 0.0458 is about 40% larger than your coded d0), then it would probably make the most sense to substitute in the boundary d0 and solve for Q.

Optimizing over 2 variables only starts to make sense if you can come up with a formula for penalty away from a "preferred" d0 compared to penalty away from a "preferred" Q.

[Devdatt Thengdi]

"Function 'lt' is not implemented for MuPAD symbolic objects."

I got this error.

[Devdatt Thengdi]

Sadly I have R2011b. Symbolic toolbox is in R2017. What should I do now?

[Walter Roberson]

Symbolic Toolbox has existed since the days of MATLAB 4 or so, but it is optional. On the other hand you would have received a different error message if you did not have it. I suspect you allowed the statement "if abs(qr - Q) > 1000 " to execute; you should put a breakpoint there.

However in R2011b Symbolic Toolbox, the option 'returnconditions' did not exist. Try putting the breakpoint and skipping those two initialization lines and doing the "syms" I showed, and running to the breakpoint, and doing

temp = subs(qr - Q, d0, 0.033);
solve(temp)

[Devdatt Thengdi]

New error: Too many output arguments. Thanks. at -

Error in Sol (line 92)
sol = solve(soltemp.conditions);

[Devdatt Thengdi]

Undefined function 'min' for input arguments of type 'sym'.

Error in Sol (line 93)
minidx = min(abs(soltemp - 700000));

[Devdatt Thengdi]

The error was resolved thanks. But, there is a new one.

Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit.  Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.

[Walter Roberson]

Please show your current code and indicate the point at which it encountered the recursion error.

[Devdatt Thengdi]

Current code:

        function qr = Sol(~)
        syms d0 nonnegative
        syms Q real
        %syms Qr real
        %Kcal/m^2degC  
        for Q = 0:100000:9000000
        THD = 3773.09; %Total heat duty Gcal/h
        %Target fluid
        Ti = 251.9; %Fluid in degC
        T0 = 386.1; %Fluid/gas out degC
        %cpl = 0.674; %Cp fluid Kcal/kgdegC
        LHV = 10.515; %Fuel value Kcal/kg
        %flowL = 12.9; %flowrate fluid
        %Air
        ea = 0.15; %excess air
        cpair = 0.171; %Kcal/kgdegK
        cpstm = 0.5125; %Kcal/kgdegK
        Tairin = 114;
        Tstmin = 141;
        Tref = 41; %Ambient
        %Fuel
        %Hea = 111.5; %Heat released by fuel Gcal/h 
        Ef = 0.9320; %Fuel efficiency
        atfr = 10.344; %in notebook
        stfr = 0.5; %literature
        %Dimensions
        Tg = 869; %degC
        Tw = 100+0.5*(Ti+T0); %Tube wall temperature
        ctc = 0.3048*20; %ctc in m
        Lh = 19.236; %Exposed tube length- horizontal m
        Lv = 18.431; %Exposed tube length- vertical m 
        d0 = 0.168; %Tube diameter m
        pi = 22/7;
        %dstk = 4.435;
        %lstk = 42.405;
        h = 2; %approximation
        sig = 0.173*10^-08; 

        %Radiant section Calcuations

        qfuel = THD/Ef;
        mfuel = qfuel/LHV;
        mreqair = mfuel*atfr;
        mair = mreqair*(1+ea);
        qair = mair*cpair*(Tairin - Tref);
        mstm = mfuel*stfr;
        qstm = mstm*cpstm*(Tstmin - Tref);
        A = pi*d0*(Lh + Lv);
        Tavg = (Ti + T0)/2; %Assumption
        qr = h*A*(Tw - Tavg);
        qout = 0.1*Q;
        qlibr = qfuel+qair+qstm-qr-qout;

        nt = qlibr/A.*Q;
        Acp = (Lh + Lv)*nt*ctc;
        Ar = Acp;
        r = ctc/d0;
        al = 1 - (0.0277+0.0927*(r - 1))*r - 1;
        %Vrad = 4*A; %4cuft/sqft of radiant area (neglecting convection section)
        %Vstk = pi*dstk*lstk;
        %VF = Vrad + Vstk;
        Di = 10.000; %Furnace Diameter
        Lm = Di; %Assuming L/D >=2

        P = 0.288 - 0.229*ea + 0.09*ea^2;
        pl = P*Lm;

        z = (Tg+460)/1000;
        a = 0.47916 - 0.19847*z + 0.022569*z*z;
        b = 0.047029 + 0.0699*z - 0.01528*z*z;
        c = 0.000803 - 0.00726*z + 0.00159*z*z;
        e = a + b*(pl) + c*pl*pl;

        Aw = 2*pi*Di*Di - Acp;

        %Exchange factor
        y = Aw/(al*Ar);
        u = 0.00064 + 0.0591*y + 0.00101*y*y;
        v = 1.0256 + 0.4908*y - 0.058*y*y;
        w = -0.144 - 0.552*y + 0.040*y*y;
        F = u + e*v + w*e*e;

        Qr = ((Tg+460)/1000.^4 - (Tw+460)/1000.^4) + h*A*pi*d0*Q*((Tg+460)/1000 - (Tw+460)/1000)/(sig*al*qlibr*ctc*F);
        qr = Qr*pi*d0*Q/(sig*al*qlibr*ctc*F);  

        temp = subs(qr - Q, d0, 0.033);
        soltemp =(solve(temp));

        minidx = min(double(abs(Sol - 700000)));
        bestQ = double(soltemp(min));

        %Q
        disp(Q);
        %d0    
        disp(d0);
        %F 
        %disp(F);
        %e 
        %disp(e);
        %pl
        %disp(pl);
        %a
        %disp(a);
        %b
        %disp(b);
        %c
        %disp(c);
        %y
        %disp(y);
        %u
        %disp(u);
        %v
        %disp(v);
        %w
        %disp(w);
        %Tw
        %disp(Tw);
        %qr    
        disp(qr);

        end

        end

Error:

Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit.  Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.

Error in iscellstr

[Walter Roberson]

I suspect you might have your own function named solve() on your path.

Your code cannot work in its current form. You have syms d0 nonnegative syms Q real

but a few lines later

for Q = 0:100000:9000000

so you overwrite Q with a specific numeric value there, and in the dimensions section of your code you have

d0 = 0.168; %Tube diameter m

So you have no remaining symbolic variables, and so cannot use solve() in that situation.