Parameter Estimation First Order ODE

2 views (last 30 days)
Nora Rafael
Nora Rafael on 13 Oct 2019
I have an ODE with 3 unknown parameters: A, D and h. I have grouped them into one parameter k, where k = AD/h for simplification and trying to estimate this k parameter by fitting to experimental data. The ODE is dCb/dt=k*(Cs-Cb)/DV - this is essentially the Noyes-Whitney which simulates dissolution of a particle in a liquid. The experimental data are concentration (Cb) vs time (t) imported from excel. When I run the below code I get an output saying ans = 4.1855e-08, I am assuming this is the estimated k. But the graph makes no sense (see below). If I manually put in this value of k, it makes a bit more sense but the fitting is poor.
ODE-N-Y.JPG
function bestk = odeParamID_NR_4()
%Experimental Data FCT0002
FCT0002 = xlsread('FCT0002');
t_exp = FCT0002(:,1);
Cb_exp = FCT0002(:,2);
plot(t_exp, Cb_exp,'o');
% ODE Information
tSpan = [0 9E4];
r0 = 1378e-09/2; %m dv50
rho = 1370; %kg/m3
W=18/1000000; %kg
N= W/((4/3)*pi*r0^3*rho);
A=(4*pi*r0^2*N); %m2
D= 5E-10; %m2/s
h= 1.416E-5; %m
DV=0.0009; % m3, dosed volume
Cb0=0; %initial value for Cb
%k initial guess
k0=A*D/h; % m3/s
%solve ODE
ODE_Sol = ode45(@(t,Cb)updateStates(t, Cb, k0), tSpan, Cb0);
simCb = deval(ODE_Sol, t_exp); % Evaluate the solution at the experimental time steps
hold on
plot(t_exp, simCb, '-b')
hold off
m_exp=(Cb_exp*DV)*1E6;% Experimental mass, mg
sim_m=(simCb*DV)*1E6; %mg Simulated mass, mg
plot(t_exp, m_exp, 'ro')
hold on
plot (t_exp, sim_m, 'r-')
%Set up optimization
myObjective = @(x) objFcn(x, t_exp, Cb_exp,tSpan,Cb0);
lb = 1e-8;
ub = 1e-4;
bestk = lsqnonlin(myObjective, k0, lb, ub);
%Plot best result
ODE_Sol = ode45(@(t,Cb)updateStates(t,Cb,bestk), tSpan, Cb0);
bestCb = deval(ODE_Sol, t_exp);
plot(t_exp, bestCb, '-g')
legend('Exp Data','Initial Param','Best Param');
function f = updateStates(t, Cb, k)
DV=0.0009; % m3, dosed volume
Cs=0.032 %kg/m3;
f = k*(Cs-Cb)/DV;
function cost = objFcn(x,t_exp,Cb_exp,tSpan,Cb0)
ODE_Sol = ode45(@(t,Cb)updateStates(t,Cb,x), tSpan, Cb0);
simCb = deval(ODE_Sol, t_exp);
cost = simCb-Cb_exp;

Answers (0)

Community Treasure Hunt

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

Start Hunting!