Parameter Estimation First Order ODE
2 views (last 30 days)
Show older comments
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.
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;
0 Comments
Answers (0)
See Also
Categories
Find more on Linear Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!