@ Dursman Mchabe: your question is very difficult to read with so much code. Please edit your question, delete the code, upload the code by clicking the paperclip button.
Too many output arguments error when using lsqcurvefit
3 views (last 30 days)
Show older comments
I have a set of 7 ordinary differential equations and 1 algebraic which I use to calculate the concentrations of 8 species with time. I have also conducted experiments to measure the concentrations of such 8 species. The challenge I have now that the concentrations predicted by the equations are not the same as the concentrations measured during experiments. As a result, I am looking at ways the regress the 4 parameters and 8 initial conditions, such that the calculated concentrations can match the experimental concentrations. I have found an example on the internet, from this link:
https://www.mathworks.com/matlabcentral/answers/329901-solving-coupled-differential-equations
in which the approach of using matlab's lsqcurvefit was followed. I Have tried to follow the step by step, however, I get the error message:
Error using EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3/ode15ifun Too many output arguments.
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3 (line 155) B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations)
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
My current code is:
function EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3
global ExperimentalTime ExperimentalConcentrations
function ode15ifun
ExperimentalTime = [600
1200
1800
2400
3000
10200
17400
24600
31800
39000
46200
53400
60600
67800
75000
82200
89400
96600
103800
111000
118200
125400
132600
139800
147000
154200
161400
168600
175800
183000];
ExperimentalConcentrations = [ 0.00E+00 6.97E-08 3.34E-01 4.88E-06 4.88E-06 4.86E+01 0.00E+00 7.41E-06
1.72E-02 8.76E-08 6.79E-01 4.88E-06 4.88E-06 4.92E+01 0.00E+00 9.33E-06
1.71E-02 1.18E-07 1.04E+00 5.11E-06 5.11E-06 4.94E+01 0.00E+00 1.20E-05
1.70E-02 1.03E-07 1.41E+00 5.11E-06 5.11E-06 4.96E+01 0.00E+00 1.05E-05
1.68E-02 1.77E-07 6.00E+00 5.35E-06 5.35E-06 3.78E+01 6.30E+00 1.74E-05
1.74E-02 3.65E-07 1.07E+01 5.35E-06 5.35E-06 3.09E+01 9.33E+00 3.72E-05
1.65E-02 3.50E-07 1.56E+01 5.36E-06 5.36E-06 2.67E+01 1.13E+01 3.55E-05
1.67E-02 8.83E-07 2.07E+01 5.36E-06 5.36E-06 1.98E+01 1.43E+01 1.00E-04
1.67E-02 5.01E-06 2.59E+01 5.07E-06 5.07E-06 1.14E+01 1.79E+01 4.07E-02
1.70E-02 5.06E-06 3.14E+01 5.07E-06 5.07E-06 9.58E+00 1.91E+01 2.45E-01
1.80E-02 4.93E-06 3.67E+01 4.93E-06 4.93E-06 5.65E+00 2.09E+01 6.17E-01
2.15E-02 5.29E-06 4.18E+01 5.30E-06 5.30E-06 3.32E+00 2.21E+01 1.32E+00
2.68E-02 5.30E-06 4.66E+01 5.30E-06 5.30E-06 0.00E+00 2.37E+01 2.29E+00
3.07E-02 5.30E-06 5.09E+01 5.30E-06 5.30E-06 0.00E+00 2.40E+01 2.34E+00
3.33E-02 9.45E-06 5.51E+01 9.46E-06 9.46E-06 0.00E+00 2.43E+01 2.40E+00
3.68E-02 9.45E-06 5.90E+01 9.46E-06 9.46E-06 0.00E+00 2.47E+01 1.82E+00
4.14E-02 1.13E-05 6.23E+01 1.13E-05 1.13E-05 0.00E+00 2.50E+01 1.38E+00
4.48E-02 1.15E-05 6.56E+01 1.15E-05 1.15E-05 0.00E+00 2.54E+01 2.09E+00
4.87E-02 1.10E-05 6.87E+01 1.10E-05 1.10E-05 0.00E+00 2.58E+01 1.82E+00
5.19E-02 1.10E-05 7.12E+01 1.10E-05 1.10E-05 0.00E+00 2.62E+01 1.58E+00
5.47E-02 1.04E-05 7.36E+01 1.04E-05 1.04E-05 0.00E+00 2.67E+01 2.29E+00
5.75E-02 1.04E-05 7.59E+01 1.04E-05 1.04E-05 0.00E+00 2.71E+01 1.62E+00
5.97E-02 1.04E-05 7.78E+01 1.04E-05 1.04E-05 0.00E+00 2.76E+01 1.12E+00
6.17E-02 1.03E-05 7.95E+01 1.03E-05 1.03E-05 0.00E+00 2.81E+01 8.91E-01
6.33E-02 1.04E-05 8.11E+01 1.04E-05 1.04E-05 0.00E+00 2.86E+01 8.51E-01
6.44E-02 1.05E-05 8.24E+01 1.06E-05 1.06E-05 0.00E+00 2.90E+01 7.59E-01
6.48E-02 1.13E-05 8.24E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 8.71E-01
6.53E-02 1.13E-05 8.23E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 1.12E+00
6.57E-02 1.13E-05 8.21E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 1.00E+00
6.59E-02 1.03E-05 8.20E+01 1.03E-05 1.03E-05 0.00E+00 2.90E+01 8.51E-01];
%-------------------------------------------------------------------------------------------------------------------------------------
% b(1:4) = Parameters, b(5:12) = Initial Conditions % MAPPING: b(1) = kga, b(2) = kLa_SO2, b(3) = ktot, b(4) = Kad, % b(5)= CSO2_gas(t0), b(6)= CCO2_gas(t0), b(7)= S_total(t0), b(8)= C_total(t0), b(9)= Ca2_total(t0),b(10)= CCaCO3(t0),b(11)= CCaSO3(t0),b(12)= CH(t0)
syms b V_Headspace F_rate CSO2_in CCO2_in R T HSO2 DCa2 DSO2 DHSO3 DSO32 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 BETCaCO3 MWCaCO3 KCO2 KHCO3 KSO2 KHSO3 Kw CCa2
eqn1 = diff(b(5),t)== F_rate/ V_Headspace * ( CSO2_in - b(5)) - (b(5) * R * T - HSO2 * ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 * b(5) * R * T / b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 * HSO2 * b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))));
eqn2 = diff(b(6) ,t)== F_rate/ V_Headspace * ( CCO2_in - b(6)) - (kLa_CO2 * ((((DCO2*(HCO2 * b(6) * R * T - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))) + DHCO3*(((KCO2 * HCO2 * b(6) * R * T/b(12)) - ((C_total (t) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * b(6) * R * T/b(12)^2) - ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * b(6) * R * T - ((C_total (t) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))))*(b(6) * R * T/HCO2 - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))));
eqn3 = diff(b(7),t)== (b(5) * R * T - HSO2*((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2*(HSO2* b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) + DHSO3*(((KSO2 * HSO2 * b(5) * R * T/b(12)) - ((b(7) * KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))) - 0.162 * exp(-5153/T) * (((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn4 = diff(b(8),t)== (kLa_CO2 * ((((DCO2*(HCO2 * b(6) * R * T - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))) + DHCO3*(((KCO2 * HCO2 * b(6) * R * T/b(12)) - ((C_total (t) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * b(6) * R * T/b(12)^2) - ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * b(6) * R * T - ((C_total (t) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))))))*(b(6) * R * T /HCO2 - ((b(8) * b(12)^2)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)))) + (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) *(1 - (b(4) * b(12))/(1 + b(4) * b(12))));
eqn5 = diff(b(9),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) *(1 - (b(4) * b(12))/(1 + b(4) * b(12)))) - (0.162 * exp(-5153/T) * (((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(b(10),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * b(10) * b(12) * (1 - (b(4) * b(12))/(1 + b(4) * b(12))));
eqn7 = diff(b(11),t)== 0.162 * exp(-5153/T) * ((((CCa2 * ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* b(5) * R * T/b(12)^2) - ((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * b(5) * R * T - ((b(7) * b(12)^2)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(b(12),t)== b(12) + 2 * CCa2 - ((b(7) *KSO2 * b(12))/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3))- 2*((b(7) * KSO2 * KHSO3)/(b(12)^2 + KSO2 * b(12) + KSO2 * KHSO3)) - ((b(8) * KCO2 * b(12))/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3)) - 2 * ((b(8) * KCO2 * KHCO3)/(b(12)^2 + KCO2 * b(12) + KCO2 * KHCO3))- Kw / b(12);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [b(5); b(6); b(7); b(8); b(9); b(10); b(11); b(12)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, b(1), DCa2, DSO2, DHSO3, DSO32, b(2), kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, b(3), BETCaCO3, MWCaCO3, b(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
V_Headspace = 1.5e-6; F_rate = 1.66667e-5; CSO2_in = 6.51332e-2; CCO2_in = 0; CCa2 = 10; R = 8.314; T = 323.15; HSO2 = 149; b(1) = 4.14e-6; DCa2 = 1.39e-9; DSO2 = 2.89e-9; DHSO3 = 2.89e-9; DSO32 = 2.89e-9; b(2) = 8.4e-4; kLa_CO2 = 9.598e-4; HCO2 = 5.15e3; DCO2 = 3.53e-9; DHCO3 = 2.89e-9; DCO32 = 2.89e-9; KSPCaSO3 = 1.07e-7; BETCaSO3 = 10; b(3) = 8.825e-3; BETCaCO3 = 12.54; MWCaCO3 = 100.0869; b(4) = 0.84; KCO2 = 1.7e-3; KHCO3 = 6.55e-8; KSO2 = 6.24; KHSO3 = 5.68e-5; Kw = 5.3e-8;
F = @(t, Y, YP) f(t, Y, YP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, b(1), DCa2, DSO2, DHSO3, DSO32, b(2), kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, b(3), BETCaCO3, MWCaCO3, b(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
vars;
%b0est = [0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
b0est = [1;1;1;1;1;1;1;1;];
bp0est = zeros(8,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[b0, bp0] = decic(F, 0, b(5:12), [], bp0est, [], opt);
[tSol,eqns] = ode15i(F, ExperimentalTime, b0, bp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
end
eqns0 = [4.14e-6; 8.4e-4; 8.825e-3; 0.84; 0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
%B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations)
fprintf(1, '\n\tParameters:\n\t\tkga = %11.3E\n\t\tkLa = %11.3E\n\t\tktot = %11.3E\n\t\tKad = %11.3E\n\tInitial Conditions:\n\t\tCSO2gas = %11.3E\n\t\tCCO2gas = %11.3E\n\t\tStotal = %11.3E\n\t\tCtotal = %11.3E\n\t\tCatotal = %11.3E\n\t\tCCaCO3 = %11.3E\n\t\tCCaSO3 = %11.3E\n\t\tCH = %11.3E\n', B);
F = @(t, B, BP) f(t, B, BP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, DCa2, DSO2, DHSO3, DSO32, kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, BETCaCO3, MWCaCO3, KCO2, KHCO3, KSO2, KHSO3, Kw);
[tSol,eqns] = ode15i(F, ExperimentalTime, B(5:12), bp0, opt);
figure(2)
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
hold on
plot(tSol, ySol, '--')
hold off
grid
legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
end
I can not spot where I am going wrong. I do not know where to fix. What can I try?
When I run the code without including the lsqcurvefit, the code does not give errors, however, the experimental results does not match the calculated values. I think that when I regress 4 input values, the calculated values will fit the experimental values. Below is the code on which I have commented-out the lsqcurvefit:
%function EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT2
function ode15ifun
ExperimentalTime = [600
1200
1800
2400
3000
10200
17400
24600
31800
39000
46200
53400
60600
67800
75000
82200
89400
96600
103800
111000
118200
125400
132600
139800
147000
154200
161400
168600
175800
183000];
ExperimentalConcentrations = [ 0.00E+00 6.97E-08 3.34E-01 4.88E-06 4.88E-06 4.86E+01 0.00E+00 7.41E-06
1.72E-02 8.76E-08 6.79E-01 4.88E-06 4.88E-06 4.92E+01 0.00E+00 9.33E-06
1.71E-02 1.18E-07 1.04E+00 5.11E-06 5.11E-06 4.94E+01 0.00E+00 1.20E-05
1.70E-02 1.03E-07 1.41E+00 5.11E-06 5.11E-06 4.96E+01 0.00E+00 1.05E-05
1.68E-02 1.77E-07 6.00E+00 5.35E-06 5.35E-06 3.78E+01 6.30E+00 1.74E-05
1.74E-02 3.65E-07 1.07E+01 5.35E-06 5.35E-06 3.09E+01 9.33E+00 3.72E-05
1.65E-02 3.50E-07 1.56E+01 5.36E-06 5.36E-06 2.67E+01 1.13E+01 3.55E-05
1.67E-02 8.83E-07 2.07E+01 5.36E-06 5.36E-06 1.98E+01 1.43E+01 1.00E-04
1.67E-02 5.01E-06 2.59E+01 5.07E-06 5.07E-06 1.14E+01 1.79E+01 4.07E-02
1.70E-02 5.06E-06 3.14E+01 5.07E-06 5.07E-06 9.58E+00 1.91E+01 2.45E-01
1.80E-02 4.93E-06 3.67E+01 4.93E-06 4.93E-06 5.65E+00 2.09E+01 6.17E-01
2.15E-02 5.29E-06 4.18E+01 5.30E-06 5.30E-06 3.32E+00 2.21E+01 1.32E+00
2.68E-02 5.30E-06 4.66E+01 5.30E-06 5.30E-06 0.00E+00 2.37E+01 2.29E+00
3.07E-02 5.30E-06 5.09E+01 5.30E-06 5.30E-06 0.00E+00 2.40E+01 2.34E+00
3.33E-02 9.45E-06 5.51E+01 9.46E-06 9.46E-06 0.00E+00 2.43E+01 2.40E+00
3.68E-02 9.45E-06 5.90E+01 9.46E-06 9.46E-06 0.00E+00 2.47E+01 1.82E+00
4.14E-02 1.13E-05 6.23E+01 1.13E-05 1.13E-05 0.00E+00 2.50E+01 1.38E+00
4.48E-02 1.15E-05 6.56E+01 1.15E-05 1.15E-05 0.00E+00 2.54E+01 2.09E+00
4.87E-02 1.10E-05 6.87E+01 1.10E-05 1.10E-05 0.00E+00 2.58E+01 1.82E+00
5.19E-02 1.10E-05 7.12E+01 1.10E-05 1.10E-05 0.00E+00 2.62E+01 1.58E+00
5.47E-02 1.04E-05 7.36E+01 1.04E-05 1.04E-05 0.00E+00 2.67E+01 2.29E+00
5.75E-02 1.04E-05 7.59E+01 1.04E-05 1.04E-05 0.00E+00 2.71E+01 1.62E+00
5.97E-02 1.04E-05 7.78E+01 1.04E-05 1.04E-05 0.00E+00 2.76E+01 1.12E+00
6.17E-02 1.03E-05 7.95E+01 1.03E-05 1.03E-05 0.00E+00 2.81E+01 8.91E-01
6.33E-02 1.04E-05 8.11E+01 1.04E-05 1.04E-05 0.00E+00 2.86E+01 8.51E-01
6.44E-02 1.05E-05 8.24E+01 1.06E-05 1.06E-05 0.00E+00 2.90E+01 7.59E-01
6.48E-02 1.13E-05 8.24E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 8.71E-01
6.53E-02 1.13E-05 8.23E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 1.12E+00
6.57E-02 1.13E-05 8.21E+01 1.13E-05 1.13E-05 0.00E+00 2.90E+01 1.00E+00
6.59E-02 1.03E-05 8.20E+01 1.03E-05 1.03E-05 0.00E+00 2.90E+01 8.51E-01];
%-------------------------------------------------------------------------------------------------------------------------------------
syms CSO2_gas(t) CCO2_gas(t) S_total(t) C_total(t) Ca2_total(t) CCaCO3(t) CCaSO3(t) CH(t) V_Headspace F_rate CSO2_in CCO2_in R T HSO2 kga DCa2 DSO2 DHSO3 DSO32 kLa_SO2 kLa_CO2 HCO2 DCO2 DHCO3 DCO32 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad KCO2 KHCO3 KSO2 KHSO3 Kw CCa2
eqn1 = diff(CSO2_gas(t),t)== (F_rate/ V_Headspace) * ( CSO2_in - CSO2_gas(t)) - (( CSO2_gas(t) * R * T - HSO2 * ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / ((1/kga) + 1/( kLa_SO2 * ((( DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 * CSO2_gas(t) * R * T/CH(t)) - ((S_total(t) * KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))) + DSO32 * (((( KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / ( DSO2 * ( HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))))));
eqn2 = diff(CCO2_gas(t) ,t)== (F_rate/ V_Headspace) * ( CCO2_in - CCO2_gas(t)) - ( kLa_CO2 * ((((DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))) + DHCO3 * (((KCO2 * HCO2 * CCO2_gas(t) * R * T/CH(t)) - ((C_total (t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2*KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * CCO2_gas(t) * R * T/CH(t)^2) - ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total (t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))))))) * (CCO2_gas(t) * R * T / HCO2 - (( C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))));
eqn3 = diff(S_total(t),t)== (( CSO2_gas(t) * R * T - HSO2 * ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / ((1/kga) + 1/(kLa_SO2 * ((( DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) + DHSO3 * ((( KSO2 * HSO2 * CSO2_gas(t) * R * T/CH(t)) - ((S_total(t) * KSO2 * CH(t))/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))) + DSO32 * (((( KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / ( DSO2 * ( HSO2 * CSO2_gas(t) * R * T - (( S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))))) - ( 0.162 * exp(-5153 / T) * ((( CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * ((((( KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3);
eqn4 = diff(C_total(t),t)== ( kLa_CO2 * ((((DCO2 * ( HCO2 * CCO2_gas(t) * R * T - ((C_total(t) * CH(t)^2) / (CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))) + DHCO3 * (((KCO2 * HCO2 * CCO2_gas(t) * R * T / CH(t)) - ((C_total (t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * CCO2_gas(t) * R * T / CH(t)^2) - (( C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))))) / ( DCO2 * (HCO2 * CCO2_gas(t) * R * T - ((C_total (t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))))))) * (CCO2_gas(t) * R * T / HCO2 - ((C_total(t) * CH(t)^2)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)))) + (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) *(1 - (Kad * CH(t))/(1 + Kad * CH(t))));
eqn5 = diff(Ca2_total(t),t)== (-1) * ( ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) * (1 - (Kad * CH(t))/(1 + Kad * CH(t)))) - (0.162 * exp(-5153 / T) * ((( CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(CCaCO3(t),t)== (-1) * (ktot * BETCaCO3 * MWCaCO3 * CCaCO3(t) * CH(t) * (1 - (Kad * CH(t)) / (1 + Kad * CH(t))));
eqn7 = diff(CCaSO3(t),t)== 0.162 * exp(-5153 / T) * ((((CCa2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3 / ((CCa2 * (((((KSO2 * KHSO3 * HSO2 * CSO2_gas(t) * R * T/CH(t)^2) - ((S_total(t) * KSO2 * KHSO3) / (CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)))))) / (DSO2 * (HSO2 * CSO2_gas(t) * R * T - ((S_total(t) * CH(t)^2)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(CH(t),t)== CH(t) + 2 * CCa2 - ((S_total(t) * KSO2 * CH(t)) / (CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)) - 2 * ((S_total(t) * KSO2 * KHSO3)/(CH(t)^2 + KSO2 * CH(t) + KSO2 * KHSO3)) - ((C_total(t) * KCO2 * CH(t))/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3)) - 2 * ((C_total(t) * KCO2 * KHCO3)/(CH(t)^2 + KCO2 * CH(t) + KCO2 * KHCO3))- Kw / CH(t);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [CSO2_gas(t); CCO2_gas(t); S_total(t); C_total(t); Ca2_total(t); CCaCO3(t); CCaSO3(t); CH(t)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHCO3, DCO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
V_Headspace = 1.5e-6; F_rate = 1.66667e-5; CSO2_in = 6.51332e-2; CCO2_in = 0; CCa2 = 10; R = 8.314; T = 323.15; HSO2 = 149; kga = 4.14e-6; DCa2 = 1.39e-9; DSO2 = 2.89e-9; DHSO3 = 2.89e-9; DSO32 = 2.89e-9; kLa_SO2 = 8.4e-4; kLa_CO2 = 9.598e-4; HCO2 = 5.15e3; DCO2 = 3.53e-9; DHCO3 = 2.89e-9; DCO32 = 2.89e-9; KSPCaSO3 = 1.07e-7; BETCaSO3 = 10; ktot = 8.825e-3; BETCaCO3 = 12.54; MWCaCO3 = 100.0869; Kad = 0.84; KCO2 = 1.7e-3; KHCO3 = 6.55e-8; KSO2 = 6.24; KHSO3 = 5.68e-5; Kw = 5.3e-8;
F = @(t, Y, YP) f(t, Y, YP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, kga, DCa2, DSO2, DHSO3, DSO32, kLa_SO2, kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, ktot, BETCaCO3, MWCaCO3, Kad, KCO2, KHCO3, KSO2, KHSO3, Kw);
vars;
y0est = [0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
yp0est = zeros(8,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [600, 183000], y0, yp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
figure(1)
plot(tSol,ySol)
hold on
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
hold off
end
%eqns0 = [4.14e-6; 8.4e-4; 8.825e-3; 0.84; 0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
%B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
%fprintf(1, '\n\tParameters:\n\t\tkga = %11.3E\n\t\tkLa = %11.3E\n\t\tktot = %11.3E\n\t\tKad = %11.3E\n\tInitial Conditions:\n\t\tCSO2gas = %11.3E\n\t\tCCO2gas = %11.3E\n\t\tStotal = %11.3E\n\t\tCtotal = %11.3E\n\t\tCatotal = %11.3E\n\t\tCCaCO3 = %11.3E\n\t\tCCaSO3 = %11.3E\n\t\tCH = %11.3E\n', B);
%F = @(t, Y, YP) f(t, B, BP, V_Headspace, F_rate, CSO2_in, CCO2_in, CCa2, R, T, HSO2, B(1), DCa2, DSO2, DHSO3, DSO32, B(2), kLa_CO2, HCO2, DCO2, DHSO3, DSO32, KSPCaSO3, BETCaSO3, B(3), BETCaCO3, MWCaCO3, B(4), KCO2, KHCO3, KSO2, KHSO3, Kw);
%[tSol,eqns] = ode15i(F, ExperimentalTime, B(5:12), bp0, opt);
%figure(2)
%plot(ExperimentalTime, ExperimentalConcentrations, 'p')
%hold on
%plot(tSol, ySol, '--')
%hold off
%grid
%legend('Exp-SO_{2-gas}','Exp-CO_{2-gas}','Exp-S_{total}','Exp-C_{total}','Exp-Ca^{2+}_{total}','Exp-CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-H^{+}', 'Mod-SO_{2-gas}','Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^{+}', 'Location','best')
%end
Answers (1)
Steven Lord
on 1 Oct 2018
The relevant line of your code included in the error message:
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT3 (line 155)
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations)
What requirements does lsqcurvefit place on the function handle you pass into it as the first input?
"Function you want to fit, specified as a function handle or the name of a function. fun is a function that takes two inputs: a vector or matrix x, and a vector or matrix xdata. fun returns a vector or matrix F, the objective function evaluated at x and xdata."
Does your function ode15ifun satisfy those requirements? Here's how it was defined:
function ode15ifun
How many input arguments does ode15ifun accept? Does it accept two?
How many output arguments does ode15ifun return? Does it return one?
Even if you corrected the signature of your ode15ifun, if you're calling lsqcurvefit inside ode15ifun with ode15ifun as input you'll get into an infinite loop unless you're very careful.
ode15ifun calls lsqcurvefit which calls ode15ifun which calls lsqcurvefit which calls ode15ifun ...
2 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)