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 SlurryCase04Feb2019tdata = [06000120001800024000300003600042000480005400060000660007200078000840009000096000102000108000114000120000126000132000138000144000150000156000162000168000174000180000];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))endendfunction 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));endfunction ph = HionpH (pH,b)global S_total C_total Ca_totalCH = 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 ;endfunction 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; endendfunction 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;endfunction 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 ;endfunction 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.
Best Answer