Hello. Ok, so I'm new to matlab and I've got a question regarding parameter estimation for a filtration model. I have 2 differential equations, that are related to an Equation that is not an ODE(ordinary differential equation), all the fitting parameters are present in all equations and all the equations are time dependent . The experimental data are related to the non ODE eq. at different times plus the initial condition. The kk's (fitting parameters) are the rate coefficients. I want to solve this system of DAE's using ode15 and then use the output to compute the experimental data minus the observed data and use these results to estimate the values of kk's using lsqcurvefit. I read that I must use ode15s for the DAE system, so I programed some codes, but the guessing parameters for kk don’t present stability, and in order to have some reasonable answer the initial guess most be really close to the answers, any help on how to solve this will be appreciated.
PS: the eq. are from “Transmembrane pressure profiles during constant flux microfiltration of bovine serum albumin” https://doi.org/10.1016/S0376-7388(02)00282-X
Thanks in advance.
function [D] = DE_dP_Ho_Zydney(kk, time )global Jo_m_s A_mem delta_P C_kg_m3 visc_agua por_mem Nt = length(time);Paso = Nt-1;Rm = delta_P/(Jo_m_s*visc_agua);F = 0.0003; IC_Aopen = A_mem ; IC_mp = 0; IC_dP_HZ = delta_P; IC = [IC_Aopen IC_mp IC_dP_HZ];[T,DV] = ode15s(@dDVdIV_Ho_Zydney,[0:time(Nt)*1.05/Paso:time(Nt)*1.05], IC); function [ dDVdIV ] = dDVdIV_Ho_Zydney( time,x ) % dAopen/dt = -alfa*Aopen*delta_P*C_kg_m3/(visc_agua*Rm*por_mem);
% dmp/dt = F*delta_P*C_kg_m3/(visc_agua*(Rm+kp*delta_P^Sz*mp+Rpo);
% delta_P = Jo_m_s*A_mem*visc_agua*Rm*(Rm+kp*delta_P^Sz*mp+Rpo)/(A_mem*Rm+Aopem*(kp*delta_P^Sz*mp+Rpo))
% with:
% Variables: x(1) = Aopen, x(2) = mp, x(3) = delta_P
% Parameters: alfa = kk(1), kp = kk(2), Rpo = kk(3) , Sz = kk(4)
xdot = zeros(3,1); xdot(1) = -kk(1).*x(1).*x(3)*C_kg_m3./(visc_agua*Rm*por_mem); xdot(2) = F*x(3)*C_kg_m3./(visc_agua*(Rm+kk(2)*x(3).^kk(4).*x(2)+kk(3))); xdot(3) = (Jo_m_s*A_mem*visc_agua*Rm*(Rm+kk(2)*x(3).^kk(4).*x(2)+kk(3))./(A_mem*Rm+x(1).*(kk(2)*x(3).^kk(4).*x(2)+kk(3))))-x(3); dDVdIV = xdot; endD = DV(:,3);D = D'; end
And I call like this:
por_mem = 0.10; A_mem = 4.1e-4; % m^2
L_mem = 10.0e-6; % m
visc_agua = 1e-3; % Pa*s
den_agua = 1000; % kg/m^3
den_par = 1350; % kg/m^3delta_P = 2400; % Pa
C_kg_m3 = 2.0; % concentracion de la solucion en kg/m^3
t_exp_s= [ 0.00 291.22 650.60 1009.98 1369.36 1734.94 2094.32 2304.99 2366.95 2428.92 2484.68 2608.61 2788.30... 2967.99 3153.87 3333.56 3519.45];delta_P_exp_pa = [ 2.344e3 3.172e3 4.344e3 6.343e3 9.515e3 1.579e4 2.689e4 3.599e4 3.916e4 4.192e4 4.544e4... 5.178e4 6.247e4 7.55e4 8.894e4 1.04e5 1.178e5];Jo_m_s = 0.00007; % Flujo inicial m/s
% Time in fuction of the experimental data, and increased by 5%
N_texp = length(t_exp_s);Paso_t = (t_exp_s(N_texp)*1.05*0.5/100);t = 0:Paso_t:t_exp_s(N_texp)*1.002; % el incremento es del 0.5% del tamanho total del tiempo
options = optimset('display','iter','DiffMinChange',1e-6,'MaxFunEvals',Inf,'TolX',1e-10); % Se puede aumentar 'TolFun',1e-10 para dar mas precision 'RobustWgtFun','bisquare', 'andrews', 'cauchy', 'fair', 'huber','logistic', 'talwar', or 'welsch
N_fig = 1;kk_alfa = 0.5; kk_kp = 1.5e13; kk_Rpo = 1.1e11; kk_Sz = 0.5;kk_ch = [kk_alfa kk_kp kk_Rpo kk_Sz];[k_HZ_rk1_aj_num,resnorm] = lsqcurvefit(@DE_dP_Ho_Zydney,kk_ch,t_exp_s,delta_P_exp_pa,[0,1.01e10,1.01e9,0],[1,5.1e13,1.1e12,1],options)L2_CB_CF = resnorm;hold onplot (t_exp_s,delta_P_exp_pa,'o','MarkerEdgeColor',[0.5 0.5 0.5],'MarkerSize',5);line (t,DE_dP_Ho_Zydney(k_HZ_rk1_aj_num,t),'linestyle','-.','color',[0 0 0]);grid on; ax = gca; ax.GridLineStyle = '--'; ax.GridColor = '[0.4 0.4 0.4]'; % Grid
xlabel ('t (s)'); ylabel ('\Delta P (psi)');text (t_exp_s(N_exp)*0.03,delta_P_exp_pa(N_exp)*0.5,'\it c_o = 10 [kg/m^{3}]')legend ('Experimental Data','Ho-Zydney Numerical Solution ','location','B');hold off
Best Answer