Global min optimisation for aorta 3-element windkessel model

Question

Hussam 2024-3-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2100156-global-min-optimisation-for-aorta-3-element-windkessel-model

评论： William Rose 2024-4-1，2:09

So, I'm currently trying to optimise a 3-element windkessel model for the aorta using lsqcurvefit (only optimises 1 parameter) and fminsearch (okay, but very dependent on the initial parameters and does not match the pressure or flow exactly or very close). So I am thinking of implementing a global min optimiser, such as the Monte Carlo method but without the large computational cost. Any recommendations?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2024-3-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2100156-global-min-optimisation-for-aorta-3-element-windkessel-model#answer_1433001

@Hussam,

Use fmincon(). I have used it to optimize a model of the circulation with up to 9 parameters, in order to match presure and flow traces recorded from the radial, carotid, and femoral arteries.

A windkessel model, by itself, cannot predict both pressure and flow. It can predict pressure from flow, or flow from pressure, but not both. To get both pressure and flow, you also must have a model of the heart itself, such as a time-varying elastance model. See, for example, this recent article by me and colleagues:

https://onlinelibrary.wiley.com/doi/10.1111/aor.14624

9 个评论
显示 7更早的评论隐藏 7更早的评论

Hussam 2024-3-29，2:05

在 MATLAB Online 中打开

I'm trying to use fmincon() as follows, but I get the following error. Kindly advise.

fmincon error nonlcon must be a function

% fmincon inputs
WK_P_ini_fmincon = [5e-10, 1e2, 1e2];
fmincon_lb = [0,0,0];
fmincon_ub = [5e-9, 1e10, 1e10];
Q_peak = 10e-4;  % Peak flow
t_s = 0.28;     % Systolic Time
t_c = 0.84;     % Cycle Time
dt=0.0001;      % (CFD/FSI) Simulation Time Step
%% Importing flow and pressure data
RefQWK = readtable('Ref_P_Q_WK_No-Osc.csv');
RefPWK = readtable('Ref_P_Q_WK_No-Osc.csv');
P_ref=RefPWK.P;                  % multiply by 0.00750062 for mmHg
P0_ref=RefPWK.P(1);              % multiply by 0.00750062 for mmHg
t_ref=RefPWK.t;
t=0:dt:t_ref(end);
P_smooth = smooth(interp1(t_ref, P_ref, t), 50);
P_spline = spline(t,P_smooth);
P_spline_dot = fnder(P_spline);
P=@(t) ppval(P_spline,t);
dPdt=@(t) ppval(P_spline_dot,t);
%Initial P condition for P equation
P0 = P(0);
%Flow equation:
Q_eq = @(t) Q_peak*sin((pi*t)./t_s).^2.*(t>=0).*(t<=t_s);
dQdt_eq = @(t) Q_peak*2*sin(pi*t./t_s).*cos(pi*t./t_s)*pi./t_s.*(t>=0).*(t<=t_s); 
Q = @(t) Q_eq(mod(t,t_c));
dQdt = @(t) dQdt_eq(mod(t,t_c));
%Initial Q condition
Q0 = Q(0);          
%% Curve fitting using the fmincon
 
options = optimoptions('fmincon','Display','iter');
A = [];
b = [];
Aeq = [];
beq = [];
f = @(x) fun_fmin(x, t, P, Q, dQdt, P0);
sol = fmincon(f,WK_P_ini_fmincon,A,b,Aeq,beq,fmincon_lb,fmincon_ub,options);
C_optim_fmincon = sol(1)
R1_optim_fmincon = sol(2)
R2_optim_fmincon = sol(3)
%% Function for fmin optimisation
function f = fun_fmin(x, time, P_in, Q, dQdt, P0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [t,y] = integration(C,R1,R2,time,Q,dQdt, P0);
    P = y(:,1);
    f = sum((P-P_in(t)).^2);
end
%% Function for ODE Integration using equation
function [t,y] = integration(C,R1,R2,time,Q,dQdt,P0)
  % For Q equation:
  ode_fun = @(t,y)((Q(t)/C)*(1+R1/R2) + R1*dQdt(t) - (y(1)-0)/(C*R2));
  
  [t,y] = ode45(ode_fun, time, P0);
  
end

William Rose 2024-3-29，2:35

在 MATLAB Online 中打开

@Hussam,

If you desire, you can explore the tools in the global optimization toolbox. Simulated annealing, particle swarm, and multiple starting point algoritms may be worth considering. The only method I have used is multiple starting points. There is now a toolbox item, multistart (here and here), to assist in this process. I started using multiple starting points before there was a global optimization toolbox, so I just wrote my own implementation (see below), and it seems to work fine. I have used it successfully to estimate parameters for cardiovascular models with up to nine parameters (resistors, capacitors, inertances, wave reflection).

You may decide not to use Matlab's multistart. There is a bit of a learning curve. It's not hard to do it yourself. Just call fmincon() with various initial guesses for the 3 parameters, and pick whichever solution gives the smallest error. For a three parameter model, the search space is a 3D cube, bounded by the lowest and highest reasonable values for each parameter. You could start at 9 points: the center of the cube, and 10% or 20% in from each of the 8 corners. Or you could start at 27 points: 10%, 50%, and 90% along each of the three axes.

For fmincon(), you don't need to use the equality or inequality constraints. You can just use [ ] for those arguments, when calling fmincon(), a shown:

x = fmincon(fun,x0,[],[],[],[],lb,ub);

Good luck.

William Rose 2024-3-30，4:36

编辑：William Rose 2024-3-30，4:54

@Hussam,

I have not used MultiStart so I cannot advise you about it.

I am not sure what you mean when you say "whenever I try to add 15 or more points...". Do you mean that the pressure and flow signals must be 14 points long or shorter? That would be strange. I have fitted signals that are hundreds or thousands of points.

When you say LB=0, I assume LB refers to lower bound, and I assume you mean LB=[0,0,0], because your model has three parameters. I would not use LB=[0,0,0], due to the possibilty of NaNs. I would specify very small positive values for the lower bound vector. For the resistances, I would specify a lower bound that is 1% of total peripheral resistance (TPR=mean pressure/mean flow). In your model data, the TPR is about 95 mmHg/(80 mL/s) ~= 1.2 mmHg-s/mL. (Convert to whatever units you are using.)

For C, notice that the time constant of decay of arterial pressure, when there is no flow, is on the order of 1 second. Therefore a very rough initial guess for C is C*TPR=1 s, i.e. C=(1 s)/TPR ~= 0.8 mL/mmHg. Choose a LB for C that is 1% of 0.8 mL/mmHg.

I plotted your P and F signals. At the end of diastole, the pressure is about 80 mmHg and is not changing at all. The fact that the "final" pressure asymptote is about 80 mmHg is not consistent with a standard 3 element windkessel model. In a standard 3 element windkessel, the pressure will decay to zero if you wait a long time. Therefore I do not expect a 3 element windkessel to fit your data well.

If you have more questions, please email me securely by clicking on the envelope icon in the pop-up window that appears when you click on the "WR" next to my posts. I may not answer for about two weeks, since I am very busy starting tomorrow.

Hussam 2024-3-30，11:10

Hi @William Rose, thank you very much for your feedback. With number of points, I meant the number of points for multistart. Sometimes fmincon() gets stuck with certain points.

Regarding C I did put down 1e-10 m^3/Pa (where your suggested calculcation is 6e-9, not very far off).

Regarding my pressure decay, would you recommend me to use a higher order model (4 element for example), or should I try to have a pressure waveform that decays to zero. What do you think? Thank you!

William Rose 2024-4-1，2:09

@Hussam,

You're welcome.

"Would you recommend me to use a higher order model (4 element for example), or should I try to have a pressure waveform that decays to zero. What do you think?"

To answer your question about how to proceed (for example, use a 4-4lement Windkessel), I would like to understand where the P and F data came from*, and what your goals are. Please email me securely by clicking on the envelope icon in the pop-up window that appears when you click on the "WR" next to my posts. I may not answer for a day or two, due to other obligations.

* I see "CFD/FSI" in a comment in your code. Did your data come from a CFD simulation involving fluid-structure interaction?

请先登录，再进行评论。