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에서 Conversion Between Symbolic and Numeric에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by