Parameter Estimation First Order ODE

조회 수: 7 (최근 30일)
Nora Rafael
Nora Rafael 2019년 10월 13일
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;

답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by