이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Parameter estimation when solving DAE using ODE15i solver
조회 수: 3 (최근 30일)
이전 댓글 표시
Hi everyone,
Is it possible to estimate parameters/curve fitting when solving solving DAE using ODE15i solver? Has anyone done it?
Regards, Dursman
채택된 답변
Torsten
2018년 9월 26일
Here is a link to the reference application:
https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting
Best wishes
Torsten.
댓글 수: 2
Dursman Mchabe
2018년 9월 27일
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?
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Data Preprocessing에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- 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)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)