Fsolve stops fsolve stopped because the sum of squared function values, r, has gradient with relative norm 0; this is less than options.Op​timalityTo​lerance

Fsolve not converging:

FUNCTION:

%%Defines all inputs and then goes into a for loop with anonymous functions for the unknown parameters
%%Equation has 7 unknown variables, No. of equations = 8

SCRIPT TO RUN THE FUNCTION:

x = @reactCoil_tuning_rev1
x0 = [ 0.023,0.8,0.4,0.64,0.57,-0.0363,0.42,-.104]
x = fsolve(fun,x0,optimset( 'Algorithm','trust-region-reflective','Display', 'final-detailed'))

ERROR:

fsolve stopped because the sum of squared function values, r, has gradient with
relative norm 0.000000e+00; this is less than options.OptimalityTolerance = 1.000000e-06.
However, r = 1.596198e+03, exceeds sqrt(options.FunctionTolerance) = 1.000000e-03.

Optimization Metric                                          Options
norm(grad r) =   0.00e+00                        OptimalityTolerance =   1e-06 (default)
r =   1.60e+03                               sqrt(FunctionTolerance) =  1.0e-03 (default)

x =

      0.4836    1.4093    3.1116    5.1671    7.3562   -1.4461   -7.1655   -4.6443

fval =

     -2.2000  -19.2098   -6.4000  -14.7215   -6.9000  -14.6362   -3.0000  -19.1888   -4.0000  -17.5961

exitflag =

      -2

Answer 1

[Walter Roberson]

The function output looks constant in the area that was examined, but the function value was not 0

You need to double-check that fun is not producing a constant value.

But first you might want to check to see what fun is, since

x = @reactCoil_tuning_rev1

is putting the function handle into x rather than into fun.

Answer 2

[Roshni Khetan]

sorry that was a typo on this forum - fun = @reactCoil_tuning_rev1

[Walter Roberson]

In any case, the output looked pretty constant when the minimizer looked at it.

Are the values expected to be naturally small? If so then you might need to change the tolerances.

[Roshni Khetan]

Are you referring to the 'error' when you say output or are you referring to the parameter 'x'?

The value of 'x' should be pretty small - of almost the same order as the guess values: [ 0.023,0.8,0.4,0.64,0.57,-0.0363,0.42,-.104]

Also, I tried changing the optimal tolerance by using optimset as below:

   [x,fval] = fsolve(fun,x0,optimset( 'Algorithm','trust-region-reflective','Display', 'final-detailed','MaxIter',1500,'TolFun',1e-50)),

but it gave me the same error and did not converge

[Walter Roberson]

There is no solution near that x0.

My probes so far suggest that some of the outputs have no zeros (but some of them do have singularities, and some of them can go complex valued.)

Answer 3

[Roshni Khetan]

Here is the function which I am trying to tune

%Inputs for tuning

function [o] = react_Coil_tuning_rev1(x)
T_w_in = [363,363,363,363]
T_w_out = [310,315,314,312];
m_w = [1.89,1.47,1.29,1.55];
T_a_in = [294,311,311,294];
T_a_out = [360,356,356,360];
m_a = [3.88,3.88,3.53,3.17];

P_w_in = 400*10^3;
P_a_in = 300*10^3;

%==Unit conversion==%
inch2meter = 0.0254; 

%== User provided geometrical inputs ==%
Tube_dia_out = (3/8)* inch2meter; % Tube outer diameter [m]
Tube_tks = 0.025* inch2meter; % Thickness of the tube [m]
HEX_height = 60*inch2meter ; % Height of the heat exchanger [m]
N_r = 12; % No. of rows deep
Tube_vert = 1*inch2meter; %  Vertical tube spacing [m]
Tube_hor = 0.866*inch2meter; %  Horizontal tube spacing [m]
HEX_width = N_r* (Tube_hor); % Width of the heat exchanger [m]
Length_coil = 108*inch2meter; % length of the HEX [m]
Fin_spacing = 14/inch2meter; % Fins per meter
Fin_pitch = 1/Fin_spacing; %distance betwee centerline of two fins
Fin_tks = 0.0095*inch2meter; % Fin thickness [m]

Tube_no =(HEX_height/Tube_vert)*N_r; % No. of tubes in each row or tubes high* no. of rows
serpentine = 2/3; %Circuitry/serpentine is the ratio of the number of tubes on any coil being fed from the header to the total number of tubes in a row
No_circuit = (HEX_height/Tube_vert)*serpentine; % no of circuits
Tube_passes = Tube_no/No_circuit; % No. of tube side fluid PASSES or no of tubes in each circuit
Tube_conduct = 400; % Thermal conductivity of the material of the tube- copper [W/mK]
Fin_conduct = 205; % Thermal conductivity of the material of the fin-aluminium [W/mK]

%==Tube geometry calculations==%
TubeEachRow = Tube_no/N_r; % No.  of tubes high
%Number of passes of water through tubes is essentially no. of tubes per circuit
Tube_dia_in=Tube_dia_out-Tube_tks*2; %Tube internal diameter Since this is dia, tks is (*2) [m]
Tube_cs_area=pi*(Tube_dia_in^2)/4; %Tube cross sectional area [m^2]
a_i = pi*Tube_dia_in*Length_coil; % inner tube surface area [m^2]
a_o = pi* Tube_dia_out*Length_coil; %outer tube surface area [m^2]

rho_w=1000;
nu_w= 3.26*10^-7
Pr_w=1.96
k_w=0.6729;
cp_w = 4200;

Tube_freeflowarea = (pi*(Tube_dia_in^2)/4)*(Tube_no/Tube_passes); % Flow area for the tube side fluid [m^2]
Q_w = m_w./rho_w; %Water flow rate [m^3/sec]
Vw=Q_w./Tube_freeflowarea; %Water velocity [m/sec]
G_w = m_w./Tube_freeflowarea; % Mass velocity [kg/m^2 sec]
Re_w=Vw.*Tube_dia_in/nu_w; %Reynolds number of water flow

rho_air=4.8;
Pr_air=0.71;
k_air=0.026;
cp_air=1010;
nu_air=3.87*10^-6;

A_fin = 2*((HEX_width*HEX_height) - ((pi*Tube_dia_out^2)/4)*Tube_no)*Fin_spacing*Length_coil + ...
    2*HEX_height*Fin_tks*Fin_spacing*Length_coil; % secondary area of heat transfer [m^2]
A_primary = pi*Tube_dia_out*(Length_coil - Fin_tks*Fin_spacing*Length_coil)*Tube_no...
    + 2*(HEX_width*HEX_height - ((pi*Tube_dia_out^2)/4)*Tube_no); % primary area of heat transfer or surface area of tubes [m^2]
A_tot = A_fin + A_primary; % total heat transfer surface area [m^2]
A_min = ((Tube_vert - Tube_dia_out)*Length_coil - (Tube_vert - Tube_dia_out)*Fin_tks*Fin_spacing*Length_coil)...
    *(HEX_height/Tube_vert); % Minimum free flow area for air [m^2]

V_air  = m_a./(rho_air*A_min); % velocity of air [m/sec]
g_max = m_a./A_min; % mass velocity [kg/m^2 sec]
V = HEX_height * HEX_width* Length_coil; % [ vol. of the heat exchanger] [m^3]
alpha = A_tot/V; % surface area density [1/m]
A_frontal = HEX_height * Length_coil; % Frontal/face area for air side [m^2]
sigma = A_min/A_frontal; % ratio of free flow to frontal area
D_h = 4*sigma/alpha; % Hydraulic diameter [m]
Re_air = g_max.*D_h./nu_air; % Reynolds number  

r = Tube_dia_in/2;
x_l =   ((Tube_vert/2)^2 + (Tube_hor^2)/2)^0.5 ;
x_m = Tube_vert/2;
rrat = 1.27*(x_m/r)*(((x_l/x_m) - 0.3)^0.5);
phi = (rrat - 1)*(1+0.35*log(rrat));

Cair=cp_air.*m_a; % air side capacity rate (W/K)
Cw=cp_w.*m_w; % water side capacity rate (W/K)
Cmin=min(Cair,Cw); Cmax=max(Cair,Cw); % Min. and Max. capacity rates (W/K)
C_r = Cmin./Cmax; % Specific heat ratio

for i = 1:5
Nu_w{i} =@(x) x(1)*(Re_w(i)^x(2))*(Pr_w^x(3)); %Turbulent HT correlation for pipe flow
h_w{i}= @(x) k_w*Nu_w{i}(x)/Tube_dia_in; %Water side convection heat transfer coefficient (W/m2.K)

J4{i} = @(x) x(4)*(Re_air(i)^((x(5))))*((Fin_pitch/D_h)^((x(6))))*((Fin_pitch/Tube_dia_out)^(x(7)))*((Fin_pitch/Tube_vert)^(x(8))) ;
h_o{i} = @(x) J4{i}(x)*g_max(i)*cp_air/(Pr_air^(2/3)); % Heat transfer coefficient

m{i} =  @(x) (2*h_o{i}(x)/(Fin_conduct*Fin_tks))^0.5;
eff_fin{i} =  @(x) tanh(m{i}(x)*r*phi)/(m{i}(x)*r*phi); % efficiency of a fin
n_ov{i} =  @(x) 1- (A_fin/A_tot)*(eff_fin{i}(x)); % Extended/overall surface efficiency

UA{i} =  @(x) 1/((1/(n_ov{i}(x)*h_o{i}(x)*A_tot)) + (1/(h_w{i}(x)*a_i*Tube_no))); % overall heat transfer coefficient [W/m^2K] 
NTU{i} =  @(x) (UA{i}(x)/Cmin(i)); % Number of  transfer units
effectiveness{i} =  @(x) 1-exp((NTU{i}(x)^0.22)*((exp(-C_r(i)*(NTU{i}(x)^0.78))-1)/C_r(i))); % Effectiveness equation for a cross flow heat exchanger 
error_air(i) = (T_a_out(i) - (T_a_in(i) -(effectiveness{i}(x)*Cmin(i)*(T_a_in(i) - T_w_in(i)))/(m_a(i)*cp_air)));  % Leaving air temperature [K]
error_water(i) = (T_w_out(i) - (T_w_in(i) + (effectiveness{i}(x)*Cmin(i)*(T_a_in(i) - T_w_in(i)))/(m_w(i)*cp_w) )); % Leaving water temperature [K]
o(2*i-1) = error_air(i);
o(2*i) = error_water(i);

end
end

[Walter Roberson]

You have

for i = 1:5

but i cannot exceed 4 when you have 4 items in those vectors you initialize at the beginning.

I am having a look at the function; it might take a while.

[Roshni Khetan]

Yeah that should be 4. I realized that I should not use anonymous functions here since I am passing 'x' as an argument - ran the code after making the changes but no change in fsolve result. i will wait for your comments to see if there is anything else which can be done

Answer 4

[DHABALESWAR MOHAPATRA]

How to fix the following error while using fsolve,

fsolve stopped because the sum of squared function values, r, is changing by less

than options.FunctionTolerance = 1.000000e-05 relative to its initial value.

However, r = 1.085832e+01, exceeds sqrt(options.FunctionTolerance) = 3.162278e-03.

[Walter Roberson]

@DHABALESWAR MOHAPATRA

Your search has found a location where the function is somewhat close to "flat" but which is not close enough to zero. For example if the function where @(x) x.^2 + 1 then if you are searching over real values the function cannot be less than 1.0 .

Sometimes you can get around that by using better initial starting points --- that is, sometimes the problem is that you have hit a local minima.

But sometimes it can be very difficult to find a starting point that works, if the function is badly scaled or is very "bumpy". Sometimes rescaling can help but sometimes there just isn't any practical numeric solution that you would be able to find using these methods.

And sometimes... there just isn't any solution.

What can help sometimes is, instead of trying to fsolve(FUNCTION) to instead construct FUNCTION2 = @(x) FUNCTION(x).^2 and then ask fmincon or ga to minimize FUNCTION2 .