esitimate parameter and curve fitting kinetic model

20 visualizaciones (últimos 30 días)
Alfiyya Isfahani
Alfiyya Isfahani el 28 de Sept. de 2022
Comentada: Alfiyya Isfahani el 29 de Sept. de 2022
Hi! I want to ask some question.
So, I want to estimate parameter from this kinetic model. The reaction is isomerization glucose to fructose with net rate constant (reversible reaction). The parameter that I want to estimate is 'k' (net rate constant) and 'K' (equilibrium constant) from experimental data.
I make this code with fminsearch and ODE45 for solve the problem, but when I change the initial guess, the result is change too. Is there anyone can help me with the code?
Thanks in advance!
function estimpara
clc
clear all
function dcdt = subpro(t,C,p)
%Unknown Parameters
K = p(1); %equilibrium constant
k = p(2); %net rate constant
%Model
dcdt = zeros(2,1);
dcdt(1) = -k.*(C(1).*(1-(1./K)));
dcdt(2) = k.*(C(1).*(1-(1./K)));
dcdt = dcdt;
end
function obj = objective(p)
%initial concentration at t=0
C0 = [0.068;0.015];
%spesific time points
ts = [1, 2, 3, 5, 8, 10, 12, 15, 20, 30];
%function handle to pass p through substrate function and obtain numerical
%model by solving nonlinear differential equation
[t,C] = ode45(@(t,C)subpro(t,C,p),ts,C0);
%experimental data at each time point ts
Cmeasured = [0.068 0.015
0.065 0.016
0.062 0.017
0.06 0.018
0.059 0.02
0.058 0.0206
0.057 0.021
0.056 0.022
0.0553 0.023
0.055 0.024];
%objective function to minimise the square of the difference
%between the numerical model and experimental data
A = rms((Cmeasured(:)-C(:))./Cmeasured(:));
obj = sum(A);
end
%Parameter Initial Guess
K = 2;
k = 1;
p0 = [K;k];
fun = @objective;
options = optimset('Display','iter','PlotFcns',@optimplotfval);
[p,fval] = fminsearch(fun,p0,options);
%optimized parameter values (untuk katalis Amine)
K = p(1);
disp(['K : ' num2str(K)])
k = p(2);
disp(['k : ' num2str(k)])
%calculate model with updated parameter
p1 = K;
p2 = k;
ts = [1, 2, 3, 5, 8, 10, 12, 15, 20, 30];
C0 = [0.068;0.015];
[t,C] = ode45(@(t,C)subpro(t,C,p),ts,C0);
Cmeasured = [0.068 0.015
0.065 0.016
0.062 0.017
0.06 0.018
0.059 0.02
0.058 0.0206
0.057 0.021
0.056 0.022
0.0553 0.023
0.055 0.024];
%plot of updated parameter values
figure (1)
plot(ts,Cmeasured(:,1),'o',ts,C(:,1),'b')
hold on
plot(ts,Cmeasured(:,2),'o',ts,C(:,2),'b')
xlabel('time')
ylabel('concentration')
end

Respuestas (1)

Bjorn Gustavsson
Bjorn Gustavsson el 28 de Sept. de 2022
Is your model correct? To my (very much non-biology-expert) eyes the ODE for c(1) seem to reduce to:
dcdt(1) = C(1)*k*(1/K - 1);
That can be re-parameterized with only one reaction-konstant, lets say B:
B = k*(1/K-1);
Fminsearch can find you best values for B but any pair of k and K that produce that single-parameter B will give the same fit to your data. That's why you get different values from different starting-points.
Since this model leads to an exponential decrease of c(1) it doesn't make much sense to talk about equilibrium and reversible reaction - unless there should also be some transition from c(2) to c(1) as well.
HTH

Productos


Versión

R2016a

Community Treasure Hunt

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

Start Hunting!

Translated by