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

1 次查看（过去 30 天）
Martin Pecha
Martin Pecha 2018-11-19
评论： Torsten 2018-11-20
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

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2018-11-19
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.
  8 个评论
显示 6更早的评论
Martin Pecha
Martin Pecha 2018-11-19
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
Torsten 2018-11-20
What do you mean ? The values for the x-vector in "Pressure_distribution" don't change ?
Did you try different initial values for x0 ?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Oil, Gas & Petrochemical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by