Parameter Estimation for a System of Differential Equations and DAEs

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

Answer 1

[Alan Weiss]

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

Answer 2

[Jhon MT]

Anyone that can help with any idea????