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

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?