MATLAB Answers

Jhon MT
0

Parameter Estimation for a System of Differential Equations and DAEs

Jhon MT
さんによって質問されました 2019 年 1 月 25 日
最新アクティビティ Alan Weiss
さんによって 回答されました 2019 年 2 月 5 日
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;
end
D = 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^3
delta_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 on
plot (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

  0 件のコメント

サインイン to comment.

2 件の回答

回答者: Alan Weiss
2019 年 2 月 5 日
 採用された回答

It is possible that the results are too noisy for direct estimation by lsqcurvefit, which is a derivative-based solver. I have two suggestions:
  • As suggested in Optimizing a Simulation or ODE, take larger-than-default finite difference steps.
  • Also as suggested in Optimizing a Simulation or ODE, try using patternsearch to do the optimization. patternsearch is much less sensitive to noise than a derivative-based optimizer. Be sure to put appropriate bounds on your control variables, as patternsearch can wander into regions of parameter spacee that you might not want.
Good luck,
Alan Weiss
MATLAB mathematical toolbox documentation

  0 件のコメント

サインイン to comment.


回答者: Jhon MT
2019 年 2 月 5 日

Anyone that can help with any idea????

  0 件のコメント

サインイン to comment.



Translated by