Error: Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Hi, I am trying to solve set of nonlinear equations res - I try to calculate flow rate distribution in system of pipes lets say using Kirchoffs law method. Idea is, that based on initial flow rates, it calculates speeds, then lambas from outer file, then pressure drops and then it solves all files using fsolve.

my objective is to give matlab initial q (volume rate) to get calculated volumes within pipes. However, it gives me this Error: Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

function res = Pressure_distribution(x, pars)

% Rozbaleni parametru
QTOT=pars(3);   % [m3/s]      - Prutok negalytu a posilytu svazkem
% Pocet jednoclanku a jejich vlastnosti
Ncell=pars(4);  % [-]         - Pocet clanku
wE=pars(8);     % [m]         - Sirka elektrodoveho prostoru clanku
hE=pars(9);     % [m]         - Vyska elektrodoveho prostoru clanku
dE=pars(10);    % [m]         - Tloustka elektrodoveho prostoru clanku
% Geometricke parametry svazku
lC=pars(11);    % [m]         - Delka kanalku
AC=pars(12);    % [m]         - Prurez kanalku
OC=pars(13);    % [m]         - Smoceny obvod kanalku
lC_R=pars(14);  % [m]         - Delka rovne casti kanalku mezi otockama
kMO_O=pars(15); % [-]         - Soucinitel mistniho odporu pro 180x otocku kanalku
lM=pars(16);    % [m]         - Delka manifoldu mezi stredy kanalku
AM=pars(17);    % [m2]        - Prurez manifoldu
OM=pars(18);    % [m]         - Smoceny obvod manifoldu
ea=pars(30);    % [m]         - absolutni drsnost povrchu
% Vlastnosti elektrolytu
rhoN=pars(21);  % [kg/m3]     - Hustota negalytu 
etaN=pars(22);  % [Pa.s]      - Dynamicka viskozita negalytu
% Vlastnosti plsti
dF=pars(27);    % [m]         - Prumer vlakna plste
etaF=pars(28);  % [-]         - Porozita plste
kKC=pars(29);   % [-]         - Carman-Kozeneho konstanta pro tok v plsti


qC1=x(1);
qC2=x(2);
qC3=x(3);
qC4=x(4);
qC5=x(5);
qIM1=x(6);
qIM2=x(7);
qIM3=x(8);
qIM4=x(9);
qOM1=x(10);
qOM2=x(11);
qOM3=x(12);
qOM4=x(13);


% Vypocet Reynoldse pro poloclanky
% Qcell=x(1:5);
vcell=x(1,1:5)/(wE*dE);
Reycell_N=vcell*dF*rhoN/(1-etaF)/etaN;
% Kontrola Reynoldsu
if Reycell_N>2300
    Povr='error - turbulentnx proudxnx v zxpornxm poloxlxnku, nutno pouxxt jinx korelace!';
    Ploss='error - turbulentnx proudxnx v zxpornxm poloxlxnku, nutno pouxxt jinx korelace!';
    return
end

%definovani konstant pro vypocet lambda
konst_flambda_CN=[AC,OC,rhoN,etaN,ea];
konst_flambda_MN=[AM,OM,rhoN,etaN,ea];

% Tlakova ztrata v kanalcich %note - nebylo by lepsi vyjadrit DPC vic
% obecne a zrusit tyto vzorecky?
% Qcell=x(1:5);
vC=x(1,1:5)/AC;
deC=4*AC/OC;

N_MO=fix(lC/lC_R);       % Urceni poctu 180x otocek kanalku
lambda_CN=Moody_diagram(Qcell,konst_flambda_CN);
DpC_N=(lambda_CN*lC/deC+N_MO*kMO_O)*rhoN*vC^2/2;
% Tlakova ztrata v manifoldu - sel by uvazovat jeste natok/vytok z/do kanalku
% QM=x(6:13);
vM=x(1,6:13)/AM;
deM=4*AM/OM;

lambda_MN=Moody_diagram(QM,konst_flambda_MN);
DpM_N=lambda_MN*(lM*Ncell)/deM*rhoN*vM^2/2;
% Tlakova ztrata ve clanku - sel by uvazovat jeste natok/vytok z/do clanku
kappaE=dF^2/16/kKC*etaF^3/(1-etaF)^2;
DpCell_N=etaN/kappaE*hE*vcell;
%Tlakova ztrata v kanalcich + clanek
DpCtotal_N=2*DpC_N+DpCell_N;
% Vypocet celkove tlakove ztraty a ztratoveho vykonu zpusobeneho tokem elektrolytu
DpTOT_N=2*DpC_N+2*DpM_N+DpCell_N;


res(1)=QTOT-qOM1-qC1;
res(2)=qC1-qIM1;
res(3)=qOM1-qC2-qOM2;
res(4)=qC2+qIM1-qIM2;
res(5)=qOM2-qC3-qOM3;
res(6)=qC3+qIM2-qIM3;
res(7)=qOM3-qC4-qOM4;
res(8)=qC4+qIM3-qIM4;
res(9)=qOM4-qC5;
res(10)=DpM_N(5)+DpCtotal_N2(2)-DpM_N(1)-DpCtotal_N2(1);
res(11)=DpM_N(6)+DpCtotal_N2(3)-DpM_N(2)-DpCtotal_N2(2);
res(12)=DpM_N(7)+DpCtotal_N2(4)-DpM_N(3)-DpCtotal_N2(3);
res(13)=DpM_N(8)+DpCtotal_N2(5)-DpM_N(4)-DpCtotal_N2(4);
res=res';
end

This is main file

x0=zeros(1,13);
x0(:)=QTOT;

x = fsolve(@Pressure_distribution,pars,x0);

This is calculation of lambda from moody diagram

function lambda = Moody_diagram(Q,konst)
% Rozbaleni parametru
A=konst(1);    % [m]         - Prurez kanalku
O=konst(2);    % [m]         - Smoceny obvod kanalku
rho=konst(3);  % [kg/m3]     - Hustota negalytu 
eta=konst(4);  % [Pa.s]      - Dynamicka viskozita negalytu
ea=konst(5);    % [m]         - absolutni drsnost povrchu
%%Vypocet Reynoldse

% Vypocet Reynoldse
v=Q/A
de=4*A/O;
Rey=de*v*rho/eta;
% Vypocet lambdy
if Rey<2300
    lambda=64/Rey;
elseif Rey>=2300
    lambda=0.25/(log10((6.81/Rey)^0.9+((ea/de)/3.7)))^2;
end
end

Answer 1

[Torsten]

Print the "res"-vector, and you'll probably find the problem. I guess some elements of the res-vector will be NaN or Inf values.

[Martin Pecha]

what does it mean? - how I can print res vector? How can I fixed that so there will be no NaN or Inf values - I am not skilled in that

[Torsten]

Remove the semicolon behind the line

res=res';

and run the code.

If there are NaN or Inf values, you have to adjust your initial guess vector x0 such that the res vector gives physical values.

[Martin Pecha]

I tried, but nothing has changed - still same error and no clue how to try to eliminate potential errors

update: res(10)=DpM_N(5)+DpCtotal_N2(2)-DpM_N(1)-DpCtotal_N2(1);....

should be DpCtotal_N(2) instead

however it does not change the error

[Torsten]

What is the output for "res" after you removed the ";" in the line "res=res';" and ran the code ?

In the command window, there should appear something like

res =

(13 numbers)

[Martin Pecha]

unfortunately it is not running I guess

it gives only this in command

Warning: Struct field assignment overwrites a value with class "double". See MATLAB R14SP2 Release Notes,
Assigning Nonstructure Variables As Structures Displays Warning, for details. 
> In createOptionFeedback at 33
  In prepareOptionsForSolver at 31
  In fsolve at 141
  In main_Design at 125 
Error using Pressure_distribution (line 4)
Not enough input arguments.

Error in fsolve (line 219)
            fuser = feval(funfcn{3},x,varargin{:});

Error in main_Design (line 125)
x = fsolve(@Pressure_distribution,pars,x0);


Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

[Torsten]

Substitute

x = fsolve(@Pressure_distribution,pars,x0);

by

x = fsolve(@(x)Pressure_distribution(x,pars),x0);

[Martin Pecha]

so somehow it works, I changed some operations to scalar, however, I think it does not calculate well - like it takes same values all the time and not working with the new one

[Torsten]

What do you mean ? The values for the x-vector in "Pressure_distribution" don't change ?

Did you try different initial values for x0 ?