Matlab getting stuck , no output

I was trying out this code for heat exchanger optimization but it isn't giving any output , its like getting stuck

function [f1 x] = ftrialr2 (x)
% Ds = x(1);%inside shell diameter
% do = x(2); %outside tube diameter
% B = x(3); %baffle spacing
 a1 = 8000;
 a2 = 259.2;
 a3 = 0.93;

 n = 1; %number of passes
 eff = 0.7; %Pump effficiency
 Nt = 60; %no. of tubes

 % tube side parameters(crude oil)
 mt = 18.80; %tube side mass flow rate
 Tci = 37.8; %cold fluid inlet temperature
 Tco = 76.7; %cold fluid outlet temperature
 rhot = 995.0; %density 
 mut = 0.00358; %viscosity
 muwt = 0.00213;
 cpt = 2.05; %heat capacity
 kt = 0.13; %thermal conductivity
 Rft = 0.00061; %fouling resistance

 st = 1.25*x(2); %tube pitch
 vt = (mt/(rhot*0.8*x(2).^2*pi/4))*n/Nt; %velocity of fluid on tube side
 di = 0.8*x(2); %tube inner diameter
 
 Ret = rhot*vt*0.8*x(2)/mut; %Reynold's number tube side
 ft = 0.079/(Ret^0.25); %Darcy friction factor
 Prt = mut*cpt/kt; %Prandtl number tube side
 
 % shell side parameters(kersoene)
 ms = 5.52;
 Thi = 199.0;
 Tho = 93.3;
 rhos = 850.0;
 cps = 2.47;
 mews = 0.0004;
 mewws = 0.00036;
 ks = 0.13;
 Rfs = 0.00061;
 ce = 0.12; %energy cost
 H = 7000;  %annual operating time in hours

 de = 4*(0.43*st.^2 - (0.5*pi*x(2).^2/4))/(0.5*pi*x(2)); %equivalent dia
 As = x(1)*x(3)*(1-x(2))/(1.25*x(2)); %shell side crooss section area
 vs = ms/(rhos*As); %velocity of fluid on shell side
 Res = ms*de/(As*mews); %shell side Reynold's number
 Prs = mews*cps/ks; %shell side Prandtl's number
 
  %Thermal Calculations
 hs = 0.36*(kt/de)*(Res^0.55)*(Prs^(1/3))*(mews/mewws)^0.14; %shell side heat transfer coefficient
 R = (Thi - Tho)/(Tco - Tci); %correction coefficient
 P = (Tco - Tci)/(Thi - Tci); %efficiency
 F = sqrt((R.^2+1)/(R.^2-1)); %correction factor
 LMTD = ((Thi - Tco) - (Tho - Tci))/(log((Thi-Tco)/(Tho-Tci)));
 Q = ms*cps*(Thi - Tho);
 L = 4.5;
 
 syms ht U A L 
 [ht, U, A, L] = solve((ht==kt/di*(3.657+(0.0677*(Ret*Prt*((di/L).^1.33)).^1/3))/(1+0.1*Prt*(Ret*(di/L)).^0.3)),U ==1/((1/hs)+Rfs+(x(2)/di)*(Rft+(1/ht))), A == Q/(U*F*LMTD), L == A/(pi*x(2)*Nt));

 if Ret < 2300
   ht = kt/di*(3.657 +(0.0677*(Ret*Prt*((di/L).^1.33).^(1/3)))/(1+0.1*Prt*(Ret*(di/L)).^0.3));
  else if Ret > 2300 && Ret < 10000
    ht = kt/di*((1+di/L).^0.67*(ft/8)*(Ret - 1000)*Prt/(1+12.7*((ft/8).^0.5)*((Prt.^(2/3))-1)));
   else if Ret > 10000
     ht = 0.027*kt/do*(Ret.^0.8)*(Prt.^(1/3)).*(mewt/mewwt).^0.14;
   end
  end
 end
 
 U = 1/((1/hs)+Rfs+(x(2)/di)*(Rft+(1/ht)));
 A = Q/(U*F*LMTD);
 L = A/(pi*x(2)*Nt);

 %Pressure drop
 %tube side pressure drop
 p = 4; %p = 2.5; %constant
 pt = 0.5*rhot*vt.^2*(L*ft/(0.8*x(2))+p)*n;
 %shell side pressure drop
 b0 = 0.72;
 fs = 2*b0*Res;
 ps = fs*(rhos*vs.^2/2)*(L/x(2))*(x(1)/de);
 f1 = (a1 + a2*A^a3)+ (ce*H*(mt*pt/rhot + ms*ps/rhos)/eff);
end

function [hist,searchdir] = runfmincon
 global hist
 hist.x = [];
 hist.fval = [];
 searchdir = [];
 x0 = [0.3 0.02 0.2]
 lb = [0.2 0.015 0.05 3];
 ub = [2 0.051 0.5 6];
 options = optimset('fmincon');
 options = optimset(options,'Display','Iter-detailed');
 options = optimset('PlotFcns',@optimplotfvalcust)
 figure
 options = optimset('PlotFcns',@optimplotxcust)
 options = optimset(options,'outputfcn',@outfun);
 nonlcon = [];
 [x,f1val] = fmincon(@ftrialr2,x0,[],[],[],[],lb,ub,nonlcon,options)
end

[Mohamed Elmalah]

Hello Good Day,

Could you find a solution for that ? because it is the same case for me

Answer 1

[Walter Roberson]

function [hist,searchdir] = runfmincon
 hist.x = [];
 hist.fval = [];
 searchdir = [];
 x0 = [0.3 0.02 0.2];
 lb = [0.2 0.015 0.05 3];
 ub = [2 0.051 0.5 6];
 options = optimset('fmincon');
 options = optimset(options,'Display','Iter-detailed');
%options = optimset(options, 'PlotFcns', {@optimplotfvalcust, @optimplotxcust});
 figure
 %options = optimset(options,'outputfcn',@outfun);
 nonlcon = [];
 [x,f1val] = fmincon(@ftrialr2,x0,[],[],[],[],lb,ub,nonlcon,options);
 hist.x = x;
 hist.fval = f1val;
end

Please re-check lb and ub: you have four entries in them but x0 is only 3 elements long and the function does not appear to expect 4 values.

I had to comment out the plotting because I do not have your customized plotting functions.

function [f1, x] = ftrialr2 (x)
    % Ds = x(1);%inside shell diameter
    % do = x(2); %outside tube diameter
    % B = x(3); %baffle spacing
    a1 = 8000;
    a2 = 259.2;
    a3 = 0.93;

    n = 1; %number of passes
    eff = 0.7; %Pump effficiency
    Nt = 60; %no. of tubes

    % tube side parameters(crude oil)
    mt = 18.80; %tube side mass flow rate
    Tci = 37.8; %cold fluid inlet temperature
    Tco = 76.7; %cold fluid outlet temperature
    rhot = 995.0; %density
    mut = 0.00358; %viscosity
    %muwt = 0.00213;
    cpt = 2.05; %heat capacity
    kt = 0.13; %thermal conductivity
    Rft = 0.00061; %fouling resistance

    st = 1.25*x(2); %tube pitch
    vt = (mt/(rhot*0.8*x(2).^2*pi/4))*n/Nt; %velocity of fluid on tube side
    di = 0.8*x(2); %tube inner diameter

    Ret = rhot*vt*0.8*x(2)/mut; %Reynold's number tube side
    ft = 0.079/(Ret^0.25); %Darcy friction factor
    Prt = mut*cpt/kt; %Prandtl number tube side

    % shell side parameters(kersoene)
    ms = 5.52;
    Thi = 199.0;
    Tho = 93.3;
    rhos = 850.0;
    cps = 2.47;
    mews = 0.0004;
    mewws = 0.00036;
    ks = 0.13;
    Rfs = 0.00061;
    ce = 0.12; %energy cost
    H = 7000;  %annual operating time in hours

    de = 4*(0.43*st.^2 - (0.5*pi*x(2).^2/4))/(0.5*pi*x(2)); %equivalent dia
    As = x(1)*x(3)*(1-x(2))/(1.25*x(2)); %shell side crooss section area
    vs = ms/(rhos*As); %velocity of fluid on shell side
    Res = ms*de/(As*mews); %shell side Reynold's number
    Prs = mews*cps/ks; %shell side Prandtl's number

    %Thermal Calculations
    hs = 0.36*(kt/de)*(Res^0.55)*(Prs^(1/3))*(mews/mewws)^0.14; %shell side heat transfer coefficient
    R = (Thi - Tho)/(Tco - Tci); %correction coefficient
    %P = (Tco - Tci)/(Thi - Tci); %efficiency
    F = sqrt((R.^2+1)/(R.^2-1)); %correction factor
    LMTD = ((Thi - Tco) - (Tho - Tci))/(log((Thi-Tco)/(Tho-Tci)));
    Q = ms*cps*(Thi - Tho);
    %L = 4.5;

    syms ht U A L
    ht = sym('ht');
    U = sym('U');
    A = sym('A');
    L = sym('L');
    %[ht, U, A, L] = ...
    sol = vpasolve((ht==kt/di*(3.657+(0.0677*(Ret*Prt*((di/L).^1.33)).^1/3))/(1+0.1*Prt*(Ret*(di/L)).^0.3)),U ==1/((1/hs)+Rfs+(x(2)/di)*(Rft+(1/ht))), A == Q/(U*F*LMTD), L == A/(pi*x(2)*Nt));
    L = double(sol.L);

    if Ret < 2300
        ht = kt/di*(3.657 +(0.0677*(Ret*Prt*((di/L).^1.33).^(1/3)))/(1+0.1*Prt*(Ret*(di/L)).^0.3));
    elseif Ret > 2300 && Ret < 10000
            ht = kt/di*((1+di/L).^0.67*(ft/8)*(Ret - 1000)*Prt/(1+12.7*((ft/8).^0.5)*((Prt.^(2/3))-1)));
    elseif Ret > 10000
            ht = 0.027*kt/do*(Ret.^0.8)*(Prt.^(1/3)).*(mewt/mewwt).^0.14;
    else
        ht = inf;
    end
 

    U = 1/((1/hs)+Rfs+(x(2)/di)*(Rft+(1/ht)));
    A = Q/(U*F*LMTD);
    L = A/(pi*x(2)*Nt);

    %Pressure drop
    %tube side pressure drop
    p = 4; %p = 2.5; %constant
    pt = 0.5*rhot*vt.^2*(L*ft/(0.8*x(2))+p)*n;
    %shell side pressure drop
    b0 = 0.72;
    fs = 2*b0*Res;
    ps = fs*(rhos*vs.^2/2)*(L/x(2))*(x(1)/de);
    f1 = (a1 + a2*A^a3)+ (ce*H*(mt*pt/rhot + ms*ps/rhos)/eff);
end

Notice in particular the change from solve() to vpasolve()

[Mohamed Elmalah]

[Mohamed Elmalah]

please I have run the code but there is one remaining error as shown. Could you please help

[Walter Roberson]

You need to use the runfmincon function.

If you need to test the calculation function by itself then you need to run it from the command line passing in a vector of three values.