Azzera filtri
Azzera filtri

Parameter estimation when solving DAE using ODE15i solver

1 visualizzazione (ultimi 30 giorni)
Dursman Mchabe
Dursman Mchabe il 26 Set 2018
Commentato: Dursman Mchabe il 27 Set 2018
Hi everyone,
Is it possible to estimate parameters/curve fitting when solving solving DAE using ODE15i solver? Has anyone done it?
Regards, Dursman

Accedi per rispondere a questa domanda.

Risposta accettata

Torsten
Torsten il 26 Set 2018
Here is a link to the reference application:
https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting
Best wishes
Torsten.
  2 Commenti
Dursman Mchabe
Dursman Mchabe il 27 Set 2018
Hi Torsten, I hope you will see this comment. While trying out the monod-kinetics-and-curve-fitting example, I came across a more suitable example where you were helping Dave Black. The link is:
https://www.mathworks.com/matlabcentral/answers/329901-solving-coupled-differential-equations
I decided to follow it instead. My code runs well with the experimentally-measured parameters, however the model does not fit the experimental results. I think some of the experimentally parameters might be erroneous, hence I need to estimate them from curve-fitting. Below is the code with experimentally-measured parameters:
%%MODEL
clear all
close all
clc
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
%plot(tSol,ySol)
%legend('CSO2-gas', 'CCO2-gas','S-total','C-total','Ca2-total','CCaCO3','CCaSO3','CH')
%%EXPERIMENTAL RESULTS
Exp_time = [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
];
Exp_CSO2_gas = [0.000
3.45e-5
3.46e-5
3.47e-5
3.46e-5
3.63e-5
3.47e-5
3.55e-5
3.58e-5
3.69e-5
3.94e-5
4.74e-5
5.96e-5
6.88e-5
7.52e-5
8.37e-5
9.47e-5
1.03e-4
1.13e-4
1.22e-4
1.29e-4
1.37e-4
1.43e-4
1.49e-4
1.54e-4
1.58e-4
1.59e-4
1.60e-4
1.61e-4
1.61e-4
];
Exp_S_total = [3.34e-1
6.79e-1
1.040
1.410
6.000
1.07e+1
1.56e+1
2.07e+1
2.59e+1
3.14e+1
3.67e+1
4.18e+1
4.66e+1
5.09e+1
5.51e+1
5.90e+1
6.23e+1
6.56e+1
6.87e+1
7.12e+1
7.36e+1
7.59e+1
7.78e+1
7.95e+1
8.11e+1
8.24e+1
8.24e+1
8.23e+1
8.21e+1
8.20e+1
];
Exp_Ca2_total = [4.88e-6
4.88e-6
5.11e-6
5.11e-6
5.35e-6
5.35e-6
5.36e-6
5.36e-6
5.07e-6
5.07e-6
4.93e-6
5.30e-6
5.30e-6
5.30e-6
9.46e-6
9.46e-6
1.13e-5
1.15e-5
1.10e-5
1.10e-5
1.04e-5
1.04e-5
1.04e-5
1.03e-5
1.04e-5
1.06e-5
1.13e-5
1.13e-5
1.13e-5
1.03e-5
];
Exp_CCaCO3 = [4.86e+1
4.92e+1
4.94e+1
4.96e+1
3.78e+1
3.09e+1
2.67e+1
1.98e+1
1.14e+1
9.580
5.650
3.320
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
];
Exp_CCaSO3 = [0.000
0.000
0.000
0.000
6.300
9.330
1.13e+1
1.43e+1
1.79e+1
1.91e+1
2.09e+1
2.21e+1
2.37e+1
2.40e+1
2.43e+1
2.47e+1
2.50e+1
2.54e+1
2.58e+1
2.62e+1
2.67e+1
2.71e+1
2.76e+1
2.81e+1
2.86e+1
2.90e+1
2.90e+1
2.90e+1
2.90e+1
2.90e+1
];
Exp_CH = [7.41e-6
9.33e-6
1.20e-5
1.05e-5
1.74e-5
3.72e-5
3.55e-5
1.00e-4
4.07e-2
2.45e-1
6.17e-1
1.320
2.290
2.340
2.400
1.820
1.380
2.090
1.820
1.580
2.290
1.620
1.120
8.91e-1
8.51e-1
7.59e-1
8.71e-1
1.120
1.000
8.51e-1
];
%%POST PROCESSING
figure
yyaxis left
plot(Exp_time,Exp_CCaCO3,'rd', tSol,ySol(:,6),'r--', Exp_time,Exp_CCaSO3,'bv', tSol,ySol(:,7),'b--',Exp_time, Exp_S_total,'ks', tSol,ySol(:,7),'k--')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','k');
hold on
plot(tSol,ySol(:,6),tSol,ySol(:,7), tSol,ySol(:,7))
yyaxis right
plot(Exp_time,Exp_CSO2_gas,'co', tSol,ySol(:,1),'c--',Exp_time,Exp_CH,'gh', tSol,ySol(:,8),'g--')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','k');
legend('Exp-CaCO_{3}','Mod CaCO_{3}','Exp-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Exp-S_{total}','Mod-S_{total}',...
'Mod-H^{+}','Exp-H^{+}','Mod-SO_{2-gas}','Exp-SO_{2-gas}','location','best')
hold off
When I try to use lsqcurvefitting to estimate the parameters that I suspect to be erroneous, I get this error message:
Error using symfun.parseString (line 50)
Invalid variable name.
Error in syms (line 225)
[name, vars] = symfun.parseString(x);
Error in SlurryAbsorptionModel2/odefit (line 69)
syms x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) 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 SlurryAbsorptionModel2 (line 151)
B = lsqcurvefit(@odefit, X0, ExperimentalTime, ExperimentalConcentrations, zeros(size(X0)), inf(size(X0)))
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
I suspect that the invalid variable names are x(1) ... x(8), however, I don't know what else to use. Below is the code:
function SlurryAbsorptionModel2
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 = odefit(b,t)
% b(1:4) = Parameters, b(5:12) = Initial Conditions
% MAPPING: x(1)= CSO2_gas(t), x(2)= CCO2_gas(t), x(3)= S_total(t), ...
% x(4)= C_total(t), x(5)= Ca2_total(t),x(6)= CCaCO3(t),x(7)= CCaSO3(t),...
% x(8)= CH(t), b(1) = kga, b(2) = kLa, b(3) = ktot, b(4) = Kad
syms x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) 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
% CONSTANT PARAMETERS
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; Introduced as b(1)
DCa2 = 1.39e-9;
DSO2 = 2.89e-9;
DHSO3 = 2.89e-9;
DSO32 = 2.89e-9;
%kLa_SO2 = 8.4e-4; Introduced as b(2)
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; Introduced as b(3)
BETCaCO3 = 12.54;
MWCaCO3 = 100.0869;
%Kad = 0.84; Introduced as b(4)
KCO2 = 1.7e-3;
KHCO3 = 6.55e-8;
KSO2 = 6.24;
KHSO3 = 5.68e-5;
Kw = 5.3e-8;
eqn1 = diff(x(1),t)== F_rate/ V_Headspace * ( CSO2_in - x(1)) - (x(1) * R * T - HSO2 * ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2 * (HSO2 * x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))) + DHSO3 * (((KSO2 * HSO2 * x(1) * R * T/x(8)) - ((x(3) * KSO2 * x(8))/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* x(1) * R * T/x(8)^2) - ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3))))))));
eqn2 = diff(x(2) ,t)== F_rate/ V_Headspace * ( CCO2_in - x(2)) - (kLa_CO2 * ((((DCO2*(HCO2 * x(2) * R * T - ((x(4) * x(8)^2)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3))) + DHCO3*(((KCO2 * HCO2 * x(2) * R * T/x(8)) - ((C_total (t) * KCO2 * x(8))/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * x(2) * R * T/x(8)^2) - ((x(4) * KCO2 * KHCO3)/(x(8)^2 + KCO2 * x(8) + KCO2 * KHCO3)))))) / (DCO2 * (HCO2 * x(2) * R * T - ((C_total (t) * x(8)^2)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3)))))))*(x(2) * R * T/HCO2 - ((x(4) * x(8)^2)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3))));
eqn3 = diff(x(3),t)== (x(1) * R * T - HSO2*((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3))) / ((1/b(1)) + 1/(b(2) * (((DSO2*(HSO2* x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3))) + DHSO3*(((KSO2 * HSO2 * x(1) * R * T/x(8)) - ((x(3) * KSO2 * x(8))/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))) + DSO32 * ((((KSO2 * KHSO3 *HSO2* x(1) * R * T/x(8)^2) - ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))))))) - 0.162 * exp(-5153/T) * (((CCa2 * ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* x(1) * R * T/x(8)^2) - ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn4 = diff(x(4),t)== (kLa_CO2 * ((((DCO2*(HCO2 * x(2) * R * T - ((x(4) * x(8)^2)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3))) + DHCO3*(((KCO2 * HCO2 * x(2) * R * T/x(8)) - ((C_total (t) * KCO2 * x(8))/(x(8)^2 + KCO2 * x(8) + KCO2 * KHCO3)))) + DCO32 * ((((KCO2 * KHCO3 * HCO2 * x(2) * R * T/x(8)^2) - ((x(4) * KCO2 * KHCO3)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3)))))) / (DCO2 * (HCO2 * x(2) * R * T - ((C_total (t) * x(8)^2)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3)))))))*(x(2) * R * T /HCO2 - ((x(4) * x(8)^2)/(x(8)^2 + KCO2 * x(8) + KCO2*KHCO3)))) + (b(3) * BETCaCO3 * MWCaCO3 * x(6) * x(8) *(1 - (b(4) * x(8))/(1 + b(4) * x(8))));
eqn5 = diff(x(5),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * x(6) * x(8) *(1 - (b(4) * x(8))/(1 + b(4) * x(8)))) - (0.162 * exp(-5153/T) * (((CCa2 * ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* x(1) * R * T/x(8)^2) - ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn6 = diff(x(6),t)== (-1) * (b(3) * BETCaCO3 * MWCaCO3 * x(6) * x(8) * (1 - (b(4) * x(8))/(1 + b(4) * x(8))));
eqn7 = diff(x(7),t)== 0.162 * exp(-5153/T) * ((((CCa2 * ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))) / KSPCaSO3) - 1)^2 * (BETCaSO3/((CCa2 * (((((KSO2 * KHSO3 *HSO2* x(1) * R * T/x(8)^2) - ((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2*KHSO3)))))) / (DSO2 * (HSO2 * x(1) * R * T - ((x(3) * x(8)^2)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))))))))/ KSPCaSO3;
eqn8 = diff(x(8),t)== x(8) + 2 * CCa2 - ((x(3) * KSO2 * x(8))/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3))- 2*((x(3) * KSO2 * KHSO3)/(x(8)^2 + KSO2 * x(8) + KSO2 * KHSO3)) - ((x(4) * KCO2 * x(8))/(x(8)^2 + KCO2 * x(8) + KCO2 * KHCO3)) - 2 * ((x(4) * KCO2 * KHCO3)/(x(8)^2 + KCO2 * x(8) + KCO2 * KHCO3))- Kw / x(8);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8];
vars = [x(1); x(2); x(3); x(4); x(5); x(6); x(7); x(8)];
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);
F = @(t, Y, YP) f(t, Y, YP,B(1:4), 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);
vars;
x0est = [0; 0; 3.3351e-1; 4.879e-6; 4.879e-6; 4.862e+1; 0; 7.413e-6];
xp0est = zeros(8,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[x0, xp0] = decic(F, 0, x0est, [], xp0est, [], opt);
[tSol,xSol] = ode15i(F, ExperimentalTime, x0, xp0, 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(@odefit, 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,B(1:4), 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,xSol] = ode15i(F, ExperimentalTime, x0, xp0, 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
What can I do?

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Data Preprocessing in Help Center e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by