Is it possible to use lsqcurvefit for parameter estimation while using ODE15i solver?

조회 수: 1 (최근 30일)
Dursman Mchabe
Dursman Mchabe 2018년 9월 2일
Hi everyone,
I am trying to estimate 14 parameters. The code below is with guessed parameters, and it does not give errors, however the model does not fit the experimental results:
clear all
close all
clc
syms x(t) y(t) z(t) A B C D E FF G H I J K LL MM N O P Q R SS
eqn1 = diff(x(t),t) == 1/ A * ((B * C - B * x(t))) - (x(t) * D * E - FF*((y(t) * z(t)^2)/(z(t)^2 + G * z(t) + G * H))) / ((1/I) + 1/(J * (((K*(FF* x(t) * D * E - ((y(t) * z(t)^2)/(z(t)^2 + G * z(t) + G*H))))) + LL * (((G * FF * x(t) * D * E/z(t)) - ((y(t) * G * z(t))/(z(t)^2 + G * z(t) + G*H)))) + MM * ((((G * H *FF* x(t) * D * E/z(t)^2) - ((y(t) * G * H)/(z(t)^2 + G * z(t) + G*H))))))) / (K * (FF * x(t) * D * E - ((y(t) * z(t)^2)/(z(t)^2 + G * z(t) + G*H)))));
eqn2 = diff(y(t),t) == (x(t) * D * E - FF*((y(t) * z(t)^2)/(z(t)^2 + G * z(t) + G*H))) / ((1/I) + 1/(J * (((K*(FF* x(t) * D * E - ((y(t) * z(t)^2)/(z(t)^2 + G * z(t) + G * H)) ) + LL * (((G * FF * x(t) * D * E/z(t)) - (( y(t) * G * z(t))/(z(t)^2 + G * z(t) + G*H)))) + MM * ((((G * H * FF * x(t) * D * E / z(t)^2) - ((y(t) * G * H)/(z(t)^2 + G * z(t) + G * H))))))/ (K * (FF * x(t) * D * E - (( y(t) * z(t)^2)/(z(t)^2 + G * z(t) + G * H))))))) - 0.162 * exp(-5153/E) * (((N * (( y(t) * G * H)/(z(t)^2 + G * z(t) + G*H)))/O) - 1)^3 * (P / ((N * (((y(t) * G * H)/(z(t)^2 + G * z(t) + G * H))))) / O));
eqn3 = z(t) + 2 * N - ((y(t) * G * z(t))/(z(t)^2 + G * z(t) + G*H)) - 2 * (( y(t) * G * H)/(z(t)^2 + G * z(t) + G*H)) - ((N * Q * z(t))/(z(t)^2 + Q * z(t) + Q*R)) - 2 * ((N * Q * R)/(z(t)^2 + Q * z(t) + Q*R)) - SS/z(t) == 0;
eqns = [eqn1 eqn2 eqn3];
vars = [x(t); y(t); z(t)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H, I, J, K, LL, MM, N, O, P, Q, R, SS);
A = 1.5e-6;
B = 1.66667e-5;
C = 6.51332e-2;
D = 8.314;
E = 323.15;
FF = 149;
G = 6.24;
H = 5.68e-5;
I = 4.14e-6;
J = 7.25E-2;
K = 2.98e-9;
LL = 2.35e-9;
MM = 1.69e-9;
N = 8;
O = 1.07e-7;
P = 10;
Q = 1.7e-3 ;
R = 6.55e-8;
SS = 5.3e-8 ;
F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H, I, J, K, LL, MM, N,O, P, Q, R, SS);
vars;
y0est = [1.65e-02; 0.34; 7.4131e-06];
yp0est = zeros(3,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
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_SO2_Headspace = [
1.65e-02
1.67e-02
1.67e-02
1.68e-02
1.70e-02
1.70e-02
1.71e-02
1.72e-02
1.74e-02
1.76e-02
1.80e-02
2.15e-02
2.68e-02
3.07e-02
3.33e-02
3.68e-02
4.14e-02
4.48e-02
4.87e-02
5.19e-02
5.47e-02
5.75e-02
5.97e-02
6.17e-02
6.33e-02
6.44e-02
6.48e-02
6.53e-02
6.57e-02
6.59e-02
];
Exp_S_total = [
0.34
0.69
1.05
1.43
6.00
10.76
15.61
20.61
25.77
31.19
36.51
41.59
46.36
50.66
54.84
58.79
62.10
65.33
68.47
71.00
73.40
75.61
77.55
79.29
80.80
82.13
82.10
82.01
81.88
81.72
];
Exp_H = [
7.4131E-06
9.33254E-06
1.20226E-05
1.04713E-05
1.7378E-05
3.71535E-05
3.54813E-05
0.0001
0.040738028
0.245470892
0.616595002
1.318256739
2.290867653
2.344228815
2.398832919
1.819700859
1.380384265
2.089296131
1.819700859
1.584893192
2.290867653
1.621810097
1.122018454
0.891250938
0.851138038
0.758577575
0.87096359
1.122018454
1
0.851138038
];
yyaxis left
plot(Exp_time, Exp_S_total,'bo')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','b')
hold on
plot(tSol,ySol(:,2),'k-.',tSol,ySol(:,3),'b--')
yyaxis right
plot(Exp_time,Exp_SO2_Headspace,'kd')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','k')
hold off
legend('Exp-S(IV)','Mod-SO_2','Mod-S(IV)', 'Exp-SO_2','location','Best')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','k')
clear all
On trying to make the model to fit the experimental results, I made the following changes to the code:
function SlurryFitted02Sept2018
syms x(t) y(t) z(t) A B C D E k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8) k(9) k(10) k(11) k(12) k(13) k(14)
eqn1 = diff(x(t),t) == 1/ A * ((B * C - B * x(t))) - (x(t) * D * E - k(1)*((y(t) * z(t)^2)/(z(t)^2 + k(2) * z(t) + k(2) * k(3)))) / ((1/k(4)) + 1/(k(5) * (((k(6)*(k(1)* x(t) * D * E - ((y(t) * z(t)^2)/(z(t)^2 + k(2) * z(t) + k(2)*k(3)))))) + k(7) * (((k(2) * k(1) * x(t) * D * E/z(t)) - ((y(t) * k(2) * z(t))/(z(t)^2 + k(2) * z(t) + k(2)*k(3))))) + k(8) * ((((k(2) * k(3) *k(1)* x(t) * D * E/z(t)^2) - ((y(t) * k(2) * k(3))/(z(t)^2 + k(2) * z(t) + k(2)*k(3)))))))) / (k(6) * (k(1) * x(t) * D * E - ((y(t) * z(t)^2)/(z(t)^2 + k(2) * z(t) + k(2)*k(3))))));
eqn2 = diff(y(t),t) == (x(t) * D * E - k(1)*((y(t) * z(t)^2)/(z(t)^2 + k(2) * z(t) + k(2)*k(3)))) / ((1/k(4)) + 1/(k(5) * (((k(6)*(k(1)* x(t) * D * E - ((y(t) * z(t)^2)/(z(t)^2 + k(2) * z(t) + k(2) * k(3))) ) + k(7) * (((k(2) * k(1) * x(t) * D * E/z(t)) - (( y(t) * k(2) * z(t))/(z(t)^2 + k(2) * z(t) + k(2)*k(3))))) + k(8) * ((((k(2) * k(3) * k(1) * x(t) * D * E / z(t)^2) - ((y(t) * k(2) * k(3))/(z(t)^2 + k(2) * z(t) + k(2) * k(3)))))))/ (k(6) * (k(1) * x(t) * D * E - (( y(t) * z(t)^2)/(z(t)^2 + k(2) * z(t) + k(2) * k(3)))))))) - 0.162 * exp(-5153/E) * (((k(9) * (( y(t) * k(2) * k(3))/(z(t)^2 + k(2) * z(t) + k(2)*k(3))))/k(10)) - 1)^2 * (k(11) / ((k(9) * (((y(t) * k(2) * k(3))/(z(t)^2 + k(2) * z(t) + k(2) * k(3)))))) / k(10)));
eqn3 = z(t) + 2 * k(9) - ((y(t) * k(2) * z(t))/(z(t)^2 + k(2) * z(t) + k(2)*k(3))) - 2 * (( y(t) * k(2) * k(3))/(z(t)^2 + k(2) * z(t) + k(2)*k(3))) - ((k(9) * k(12) * z(t))/(z(t)^2 + k(12) * z(t) + k(12)*k(13))) - 2 * ((k(9) * k(12) * k(13))/(z(t)^2 + k(12) * z(t) + k(12)*k(13))) - k(14)/z(t) == 0;
eqns = [eqn1 eqn2 eqn3];
vars = [x(t); y(t); z(t)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, A, B, C, D, E, k(1), k(2), k(3), k(4), k(5), k(6), k(7), k(8), k(9), k(10), k(11), k(12), k(13), k(14));
A = 1.5e-6;
B = 1.66667e-5;
C = 6.51332e-2;
D = 8.314;
E = 323.15;
F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, k(1), k(2), k(3), k(4), k(5), k(6), k(7), k(8), k(9),k(10), k(11), k(12), k(13), k(14));
vars;
y0est = [1.65e-02; 0.34; 7.4131e-06];
yp0est = zeros(3,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
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_SO2_Headspace = [
1.65e-02
1.67e-02
1.67e-02
1.68e-02
1.70e-02
1.70e-02
1.71e-02
1.72e-02
1.74e-02
1.76e-02
1.80e-02
2.15e-02
2.68e-02
3.07e-02
3.33e-02
3.68e-02
4.14e-02
4.48e-02
4.87e-02
5.19e-02
5.47e-02
5.75e-02
5.97e-02
6.17e-02
6.33e-02
6.44e-02
6.48e-02
6.53e-02
6.57e-02
6.59e-02
];
Exp_S_total = [
0.34
0.69
1.05
1.43
6.00
10.76
15.61
20.61
25.77
31.19
36.51
41.59
46.36
50.66
54.84
58.79
62.10
65.33
68.47
71.00
73.40
75.61
77.55
79.29
80.80
82.13
82.10
82.01
81.88
81.72
];
Exp_H = [
7.4131E-06
9.33254E-06
1.20226E-05
1.04713E-05
1.7378E-05
3.71535E-05
3.54813E-05
0.0001
0.040738028
0.245470892
0.616595002
1.318256739
2.290867653
2.344228815
2.398832919
1.819700859
1.380384265
2.089296131
1.819700859
1.584893192
2.290867653
1.621810097
1.122018454
0.891250938
0.851138038
0.758577575
0.87096359
1.122018454
1
0.851138038
];
k0 = [149,6.24,5.68e-5,4.14e-6,7.25E-2,2.98e-9,2.35e-9,1.69e-9,8,1.07e-7,10,1.7e-3,6.55e-8,5.3e-8];
k = lsqcurvefit(F,k0,Exp_time,ydata,Exp_SO2_Headspace,Exp_S_total,Exp_H);
times = linspace(Exp_time(1),Exp_time(end));
plot(xdata,ydata,'ko',times,fun(k,times),'b-')
yyaxis left
plot(Exp_time, Exp_S_total,'bo',times,fun(k,times),'b-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','b')
hold on
plot(tSol,ySol(:,2),'k-.',tSol,ySol(:,3),'b--')
yyaxis right
plot(Exp_time,Exp_SO2_Headspace,'kd')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','k')
hold off
legend('Exp-S(IV)','Mod-SO_2','Mod-S(IV)', 'Exp-SO_2','location','Best')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
set(gca,'YColor','k')
end
Am I in the right direction?

이 질문에 답변하려면 로그인하십시오.

답변 (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