Estimation of parameters and initial conditions using lsqcurvefit

2 次查看（过去 30 天）
Dursman Mchabe
Dursman Mchabe 2018-9-28
评论： Dursman Mchabe 2018-9-30
Hi everyone, I am trying to follow the example from the link:
https://www.mathworks.com/matlabcentral/answers/329901-solving-coupled-differential-equations
to estimate parameters and initial conditions. However, I get the error message:
Error using symfun.parseString (line 50)
Invalid variable name.
Error in syms (line 225)
[name, vars] = symfun.parseString(x);
Error in EstimationOfParametersAndInitialConditions/ode15ifun (line 65)
syms b(5) b(6) b(7) b(8) b(9) b(10) b(11) b(12) V_Headspace F_rate CSO2_in CCO2_in 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 CCa2
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in EstimationOfParametersAndInitialConditions (line 149)
B = lsqcurvefit(@ode15ifun, X0, ExperimentalTime, ExperimentalConcentrations, zeros(size(X0)), inf(size(X0)))
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
My code is as follows:
function EstimationOfParametersAndInitialConditions
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];
%-------------------------------------------------------------------------------------------------------------------------------
function X = ode15ifun(b,t)
syms b(5) b(6) b(7) b(8) b(9) b(10) b(11) b(12) V_Headspace F_rate CSO2_in CCO2_in 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 CCa2
% 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(t), b(6)= CCO2_gas(t), b(7)= S_total(t), b(8)= C_total(t), b(9)= Ca2_total(t),b(10)= CCaCO3(t),b(11)= CCaSO3(t),b(12)= CH(t)
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;
%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, b(5:12), [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [600, 183000], y0, yp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
end
%-------------------------------------------------------------------------------------------------------------------------------
X0 = [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, X0, ExperimentalTime, ExperimentalConcentrations, zeros(size(X0)), inf(size(X0)))
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, 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);
[tSol,ySol] = ode15i(F, [600, 183000], y0, yp0, opt);
figure(1)
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
hold on
plot(tSol, X, '--')
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 my mistake. Where am I going wrong?
  2 个评论
Torsten
Torsten 2018-9-28
Your function "ode15ifun" does not seem to return anything to lsqcurvefit. At least I don't see a line where you assign a value to "X".
Dursman Mchabe
Dursman Mchabe 2018-9-30
Hi Torsten, over the past two days I tried different approaches to emulate Dave Black's method:
https://www.mathworks.com/matlabcentral/answers/329901-solving-coupled-differential-equations
But now, I get the error message:
Index exceeds array bounds.
Error in sym/subsref (line 859)
R_tilde = builtin('subsref',L_tilde,Idx);
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT1/ode15ifun (line 75)
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))))))));
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT1 (line 155)
B = lsqcurvefit(@ode15ifun, eqns0, ExperimentalTime, ExperimentalConcentrations, zeros(size(eqns0)), inf(size(eqns0)))
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
My current code is:
function EstimationOfParametersAndInitialConditionsUsingLSQCURVEFIT1
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];
%-------------------------------------------------------------------------------------------------------------------------------------
function eqns = ode15ifun(b,t,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)
% 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(t), b(6)= CCO2_gas(t), b(7)= S_total(t), b(8)= C_total(t), b(9)= Ca2_total(t),b(10)= CCaCO3(t),b(11)= CCaSO3(t),b(12)= CH(t)
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, [600, 183000], 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)))
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, [600, 183000], B(5:12), bp0, opt);
figure(1)
plot(ExperimentalTime, ExperimentalConcentrations, 'p')
hold on
plot(tSol, eqns, '--')
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
What can I try?

请先登录，再进行评论。

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by