Why do I get an error message "Function value at starting guess must be finite and real." when my starting guess is 8?

3 vues (au cours des 30 derniers jours)
Dursman Mchabe
Dursman Mchabe le 4 Fév 2019
Commenté : Dursman Mchabe le 4 Fév 2019
Hi everyone,
When I run the code given below, I get, I get the error message:
Error using fzero (line 328)
Function value at starting guess must be finite and real.
Error in SlurryCase04Feb2019>SO2_OdeDriver (line 189)
pH = fzero(@(pH)HionpH(pH,b),pH_trial);
Error in SlurryCase04Feb2019>kinetics (line 155)
y = SO2_OdeDriver(y0,b);
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in SlurryCase04Feb2019 (line 77)
[b]=lsqcurvefit(@kinetics,b0,tdata,ydata);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
The full code is given below:
function SlurryCase04Feb2019
tdata = [0
6000
12000
18000
24000
30000
36000
42000
48000
54000
60000
66000
72000
78000
84000
90000
96000
102000
108000
114000
120000
126000
132000
138000
144000
150000
156000
162000
168000
174000
180000];
ydata = [0.076298787 0 0 0 0 47.61362376 0
0.019461171 0.029472654 6.668799696 16.23121199 9.959999844 35.99050802 6.300684797
0.020167931 0.042887117 11.84949115 24.04825592 14.75631477 29.46171885 9.33482827
0.019107791 0.051014128 17.25229685 29.01821217 17.80549194 25.47239217 11.26373435
0.019107791 0.051014128 17.25229685 29.01821217 17.80549194 25.47239217 11.26373435
0.019343503 0.063743181 22.80780219 36.81976241 22.59189684 18.84582265 14.29160875
0.019343503 0.078528466 28.5615147 46.11097327 28.29200863 10.86857406 17.89749311
0.019696506 0.081955518 34.59038619 49.09849644 30.12402642 9.128367823 19.05642553
0.020850931 0.088613792 40.44036333 53.76834115 32.98847955 5.388343709 20.86847538
0.024902168 0.092530424 46.03725425 56.87434597 34.89336965 3.160034015 22.07350674
0.024902168 0.092530424 46.03725425 56.87434597 34.89336965 3.160034015 22.07350674
0.031049924 0.097915792 51.29729887 61.13770121 37.50806004 0 23.72755696
0.035572586 0.097915792 56.06755651 61.90018578 37.97509734 0 24.02300424
0.038601396 0.097915792 60.70353205 62.73848974 38.48857555 0 24.34782998
0.042662438 0.097915792 65.10035551 63.65803386 39.05181501 0 24.70413464
0.047985767 0.097915792 68.80560315 64.36554967 39.48518283 0 24.97828263
0.047985767 0.097915792 68.80560315 64.36554967 39.48518283 0 24.97828263
0.051943095 0.097915792 72.44446145 65.33369776 40.07819321 0 25.35342034
0.056347712 0.097915792 76.00361827 66.58555956 40.84498401 0 25.83849134
0.060163613 0.097915792 78.89271653 67.55634934 41.43961248 0 26.21465265
0.063414181 0.097915792 81.66324156 68.6913213 42.13480587 0 26.65443122
0.06666475 0.097915792 84.24007982 69.93531357 42.89677643 0 27.13645294
0.06666475 0.097915792 84.24007982 69.93531357 42.89677643 0 27.13645294
0.069208558 0.097915792 86.52721476 71.15222604 43.64216005 0 27.60798179
0.071540791 0.097915792 88.61843725 72.41217703 44.41390563 0 28.09618718
0.073377839 0.097915792 90.46782755 73.63940812 45.16560962 0 28.57171429
0.074673315 0.097915792 92.12354968 74.82827386 45.89381405 0 29.0323756
0.07512098 0.097915792 92.26733175 74.82827386 45.89381405 0 29.0323756
0.07512098 0.097915792 92.26733175 74.82827386 45.89381405 0 29.0323756
0.075709695 0.097915792 92.33853494 74.82827386 45.89381405 0 29.0323756
0.076063075 0.097915792 92.36679575 74.82827386 45.89381405 0 29.0323756];
% Initial guesses
% b0=[8.314;149;5.15e+3; 1.39e-9;2.89e-9;3.53e-9;1.7e-3;6.55e-8;6.24;5.68e-5;5.3E-8;3.1e-4;10;8.825e-3;12.54;100.0869;8.4e-3;2703;...
% 258.30;2540;2e-08;1.01E+05;323.15;4.e-4;1.66667e-5;2e-3;3.257E-02;0.0; 0.0;4.74e-5;6e-4;9.598e-5];
%b0 = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b0 = [1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23;1e-23];
[b]=lsqcurvefit(@kinetics,b0,tdata,ydata);
fprintf(1,'\tParameters:\n')
for k1 = 1:length(b)
fprintf(1, '\t\tP(%d) = %8.5f\n', k1, b(k1))
end
end
function Results = kinetics(y,b)
global S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 Trial CH ...
nloop Delt t CH_trial CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter CCO3_inter pHtrial pi pHInter
format shorte
%% Initial conditions and Relationships
pHtrial = 7;
Ca_total = 0; %4.879E-02;
C_total = 0; %4.879E-02;
S_total = 0;
CCaCO3 = 47.61362376; %48.624- Ca_total;
CCaSO3 = 0;
fr = 1.;
SO2_g = fr * b(27);
CO2_g = b(29);
%CO2_g = 0.2;
CCO2_inter = CO2_g ;
pg = 200;
CSO2_inter = pg/b(2);
CH_trial = 1.e-7;
pH1 = 7.0;
pH = fzero(@(pH)HionpH(pH,b),pH1);
CH = 10^(3-pH);
CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + b(9)*CH + b(9)*b(10));
CHSO3 = (S_total*b(9)*CH)/(CH^2 + b(9)*CH + b(9)*b(10));
CSO3 = (S_total*b(9)*b(10))/(CH^2 + b(9)*CH + b(9)*b(10));
CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + b(7)*CH + b(7)*b(8));
CHCO3 = (C_total*b(7)*CH)/(CH^2 + b(7)*CH + b(7)*b(8));
CCO3 = (C_total*b(7)*b(8))/(CH^2 + b(7)*CH + b(7)*b(8));
pH2 = 6.0;
pHInter = fzero (@(pH)HionStar(pH,b),pH2);
CH = 10^(3-pH) ;
y0 = [ SO2_g CO2_g S_total C_total Ca_total CCaCO3 CCaSO3 ];
y0 = y0';
Trial(1) = (SO2_g * b(1) *b(23) ) / b(2);
Trial(2) = 1;
Trial(3) = CH ;
%% Time
t0 = 0.0;
Delt = 10;
t = t0;
tmax = 5806;
nloop = tmax/Delt ;
nstep = 31;
%% Results
Results = [ t y0' pH pHInter];
%% Loop1
for k = 1:nstep
y = SO2_OdeDriver(y0,b);
pH = 3 - log10 (CH) ;
tf = t;
Results = [Results; tf y' pH pHInter];
y0 = y;
end
Results;
pH = 3 - log10 (Results(:,9));
pHInt = 3 - log10 (Results(:,10));
end
function ph = HionpH (pH,b)
global S_total C_total Ca_total
CH = 10^(3-pH);
ph = CH + 2 * Ca_total - ((S_total*b(9)*CH)/(CH^2 + b(9)*CH + b(9)*b(10)))...
- 2*((S_total*b(9)*b(10))/(CH^2 + b(9)*CH + b(9)*b(10)))...
-((C_total*b(7)*CH)/(CH^2 + b(7)*CH + b(7)*b(8)))...
- 2*((C_total*b(7)*b(8))/(CH^2 + b(7)*CH + b(7)*b(8)))- b(11) /CH ;
end
function y = SO2_OdeDriver (y0,b)
global S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 CH nloop Delt t CH_trial Tag pHInter
for k = 1:nloop
pH_trial = 8.0;
pH = fzero(@(pH)HionpH(pH,b),pH_trial);
CH = 10^(3-pH);
pH = 3 -log10 (CH) ;
CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + b(9)*CH + b(9)*b(10));
CHSO3 = (S_total*b(9)*CH)/(CH^2 + b(9)*CH + b(9)*b(10));
CSO3 = (S_total*b(9)*b(10))/(CH^2 + b(9)*CH + b(9)*b(10));
CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + b(7)*CH + b(7)*b(8));
CHCO3 = (C_total*b(7)*CH)/(CH^2 + b(7)*CH + b(7)*b(8));
CCO3 = (C_total*b(7)*b(8))/(CH^2 + b(7)*CH + b(7)*b(8));
dydt = odes (t, y0,b);
y = y0 + dydt' * Delt;
SO2_g = y(1);
CO2_g = y(2);
S_total = y(3) ;
C_total = y(4) ;
Ca_total = y(5);
CCaCO3 = y(6);
CCaSO3 = y(7);
y0 = y ;
t = t +Delt;
CH_trial = CH;
end
end
function dcdt = odes (t,c,b)
global kCaSO3 CSO2_bulkslurry CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 CH S_total C_total Ca_total Trial CSO2_inter ...
CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter CCO3_inter pHtrial pi pHInter pg
S_total = c(3) ;
C_total = c(4) ;
Ca_total = c(5) ;
pH1 = 7;
pH = fzero (@(pH)HionpH(pH,b),pH1);
CH = 10^(3-pH) ;
CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + b(9)*CH + b(9)*b(10));
CHSO3 = (S_total*b(9)*CH)/(CH^2 + b(9)*CH + b(9)*b(10));
CSO3 = (S_total*b(9)*b(10))/(CH^2 + b(9)*CH + b(9)*b(10));
CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + b(7)*CH + b(7)*b(8));
CHCO3 = (C_total*b(7)*CH)/(CH^2 + b(7)*CH + b(7)*b(8));
CCO3 = (C_total*b(7)*b(8))/(CH^2 + b(7)*CH + b(7)*b(8));
CH;
pH;
p_exit = c(1)*b(1)*b(23) ;
pg = p_exit;
Trial ;
Tnew = fsolve (@(Tr)NSub(Tr,b),Trial);
pHInter = 3-log10(Tnew(3) );
pstar = Tnew(1) * b(2);
NSO2 = b(30) *(pg-pstar) ;
Trial = Tnew;
dcdt(1) = (1/b(24)) * (b(25) * b(27) - b(25) * c(1) ) - NSO2 ;
E_CO2 = 1;
NCO2 = 9.598e-5*E_CO2*( c(2)*b(1)*b(23)/b(3) - CCO2_bulkslurry);
dcdt(2) = (1/b(24)) * (b(25) * b(29) - b(25) * c(2)) - NCO2 ;
n = 3;
sat = c(5)*CSO3/ b(12) ;
if ( sat < 1.0 )
RCaSO3 = 0.0;
else
RCaSO3 = b(33) * ( sat -1)^n ;
end
RCaCO3 = b(14)*b(15)*b(16)*CCaCO3*CH *(1 - (b(17)*CH)/(1+b(17)*CH));
dcdt(3) = NSO2 - RCaSO3;
dcdt(4) = NCO2 + RCaCO3;
dcdt(5) = RCaCO3 - RCaSO3;
if ( c(6) < 1.0e-3 )
c(6) = 0.0;
dcdt(6) = 0.0;
Time_c = t ;
Tag = 1;
else
dcdt(6) = - RCaCO3 ;
end
dcdt(7) = RCaSO3;
end
function discrep = NSub (Tr,b)
global S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter CCO3_inter ...
CCO2_bulkslurry CHCO3 CCO3 pg pHInter
Tr;
CSO2_inter = Tr(1) ;
pstar = CSO2_inter * b(2) ;
S1Total = Tr(2) ;
CH = Tr(3) ;
S_total = CHSO3 + CSO3 + CSO2_bulkslurry ;
CHSO3_inter = b(9)*CSO2_inter / CH ;
CSO3_inter = b(10) *CHSO3_inter /CH;
Ctotal = CHCO3 + CCO3 ;
CHCO3_inter = Ctotal * ( 1+ b(8) / CH) ;
CCO3_inter = b(8) *CHCO3_inter /CH ;
S_total_inter = CSO2_inter + CHSO3_inter + CSO3_inter;
df = S_total_inter- S_total ;
pstar;
pg;
Rate_gas = b(30) * (pg-pstar) ;
Rate_liq = b(31) * df;
discrep(1) = S1Total - (CHSO3_inter + CSO3_inter) ;
discrep(2) = Rate_gas - Rate_liq ;
mCH_cal = 2*Ca_total -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter + b(11) /CH ) ;
discrep(3) = CH + mCH_cal ;
end
function ph1 = HionStar (pH,b)
global S_total C_total Ca_total CSO2_bulkslurry CHSO3 CSO3 CSO2_inter CHSO3_inter CSO3_inter CCO2_inter CHCO3_inter ...
CCO3_inter CCO2_bulkslurry CHCO3 CCO3 pg
CH = 10^(3-pH);
CSO2_inter ;
CHSO3_inter = b(9)*CSO2_inter / CH ;
CSO3_inter = b(10) *CHSO3_inter /CH;
Stotal = CHSO3 + CSO3;
Ctotal = CHCO3 + CCO3 ;
CHSO3_inter = Stotal * ( 1+ b(10) / CH) ;
CSO3_inter = b(10) * CHSO3_inter /CH;
CHCO3_inter = Ctotal * ( 1+ b(8) / CH);
CCO3_inter = b(8) *CHCO3_inter /CH ;
Ca = Ca_total;
ph1 = CH + 2*Ca -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter + b(11) /CH );
end
I have also attached the code.
  3 commentaires
Afficher 1 commentaire plus ancien
Dursman Mchabe
Dursman Mchabe le 4 Fév 2019
Thanks a lot for the comment.
HionPh is supposed to be evaluated by lines 171 through 180. i.e:
function ph = HionpH (pH,b)
global S_total C_total Ca_total
CH = 10^(3-pH);
ph = CH + 2 * Ca_total - ((S_total*b(9)*CH)/(CH^2 + b(9)*CH + b(9)*b(10)))...
- 2*((S_total*b(9)*b(10))/(CH^2 + b(9)*CH + b(9)*b(10)))...
-((C_total*b(7)*CH)/(CH^2 + b(7)*CH + b(7)*b(8)))...
- 2*((C_total*b(7)*b(8))/(CH^2 + b(7)*CH + b(7)*b(8)))- b(11) /CH ;
end
The initial value of HionPH is given by line 110. i.e.
pH = fzero(@(pH)HionpH(pH,b),pH1);
I am sceptical to use , pH = 7 directly, because it could override the function in lines 171 through 180.
I am too affraid to temper with the global variables, because I am relatively new in matlab, and the code was developed by my teacher/mentor. If they are a problem I will be forced to temper with them.
Dursman Mchabe
Dursman Mchabe le 4 Fév 2019
I have just tried:
%pH1 = 7.0;
%pH = fzero(@(pH)HionpH(pH,b),pH1);
pH = 7;
CH = 10^(3-pH);
For the initial value of HionpH, but the error still remain.

Connectez-vous pour commenter.

Connectez-vous pour répondre à cette question.

Réponse acceptée

Guillaume
Guillaume le 4 Fév 2019
Modifié(e) : Guillaume le 4 Fév 2019
Well, as said, global variables makes it difficult to debug. I'm certainly not going to bother trying to figure what are the values of the global variables at the point where you're calling your function.
Clearly though, you're not understanding the problem, pH is not an input, it's the output that is non-finite (or non-real). The problem is that fzero is calling HionPH with an initial of pH_trial and b, whatever their actual value is. When called with these values, HionPH either return NaN or +/-Inf.
So either pH_Trial or b are incorrect and outside the domain or the function, or somethin is wrong with the function itself. Maybe CA_total, or C_Total, or S_Total is wrong. You'll have to track where they're created and modified to check they're what you're expect. Since they're global variables they may have been modified inadvertently by any function.
Possibly the best thing to do is to issue dbstop if error at the command line and run the code. It will break into the debugger when it errors. At this point, look at the values of
ph_Trial
b
CA_Total
C_Total
S_Total
and see the result of
HionPH(pH_Trial, b) %will most likely return Nan or +/- Inf
and work from there to figure out what went wrong

Plus de réponses (0)

Catégories

En savoir plus sur Function Definition dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by