3-Element Windkessel Model (Optimizing against flow (Q)): Why is my code so slow? It does not reach the optimizations.

Hi, so I've wrote this code to optimize the parameters of a 3-element windkessel model. It was working when I was optimizing against the pressure (P), but when I switched it to flow (Q), it gets stuck calculating 'ppval' - it does not even reach the optimisation. Please advise. I have added the code and profiler results below. Is it the ODE function? I am suspecting that the ODE for Q is stiff, whereas for P it was not (even though it is the same equation)? How can I overcome this? Should I attempt the ODE solvers for stiff problems (e.g. ode15s) or for more accurate nonstiff problems (e.g. ode78)? Thanks!

The Windkessel model (P_out = 0):

%% Code to Optimize Parameters for Aorta 3-Element Windkessel Model
%%
clear all
close all
clc

%% Inputs

% For mmHg
% R1 * (0.000001/133.322);
% R2 * (0.000001/133.322);
% C  * (133.322/0.000001);

R1 = 4.40e+06;
R2 = 1.05e+08;
C  = 1.31e-08;

% Initial conditions for WK-3 parameters
IC = [C, R1, R2];

Q_peak = 5e-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
P0 = P(0);

%% Activate for flow data:
% Q_ref=RefQWK.Q;            % multiply by e+6 for mL/s
% Q0_ref=RefQWK.Q(1);
% Q_smooth = smooth(interp1(t_ref, Q_ref, t), 50);
% Q_spline = spline(t,Q_smooth);
% Q_spline_dot = fnder(Q_spline);
% 
% Q=@(t) ppval(Q_spline,t);
% dPdt=@(t) ppval(Q_spline_dot,t);
% 
%Initial P condition
% Q0 = P(0);

%Activate for flow equation:
Q_eq = @(t)Q_peak*sin((pi*t)./t_s).^2.*(t>=0).*(t<=t_s);
Q = @(t)Q_eq(mod(t,t_c));
%Initial P condition
Q0 = Q(0);          

%% Curve fitting using the method of least squares

LB = [1e-15 1 1]; UB = [1 1e10 1e10]; % Lower and Upper bounds
options = optimoptions('lsqcurvefit',  'Algorithm', 'levenberg-marquardt', 'MaxFunEvals',10000,'MaxIter',10000,'Display','iter');

f = @(x,t) fun_lsq(x, t, P, dPdt, Q0);
sol = lsqcurvefit(f, IC, t, Q(t), LB, UB, options);

C_optim_lsq = sol(1) 
R1_optim_lsq = sol(2) 
R2_optim_lsq = sol(3)

%% Curve fitting using the fminsearch
 
options = optimset('MaxFunEvals', 1e100, 'TolX', 1e-16, 'TolFun', 1e-16, 'Display','iter','PlotFcns',@optimplotfval);

f = @(x) fun_fmin(x, t, Q(t), P, dPdt, Q0);
sol = fminsearch(f,IC,options);

C_optim_fmin = sol(1)
R1_optim_fmin = sol(2)
R2_optim_fmin = sol(3)

%% Solving and Plotting ODE for Q

[t,Qsim_sterg] = integration(1.31e-8,4.4e+6,1.05e+8,t,P,dPdt,Q0);
[t,Qsim_lsq] = integration_data(C_optim_lsq,R1_optim_lsq,R2_optim_lsq,t,P,dPdt,Q0);
[t,Qsim_fmin] = integration_data(C_optim_fmin,R1_optim_fmin,R2_optim_fmin,t,P,dPdt,Q0);
[t,Qsim_fmin_manual] = integration_data(C,R1,R2,t,P,dPdt,Q0);

% PWK for FV
[t_FV_fmin, P_FV_fmin] = PWK_FV(C_optim_fmin, R1_optim_fmin, R2_optim_fmin, P0, Q(t), dt, t);
[t_FV_lsq, P_FV_lsq] = PWK_FV(C_optim_lsq, R1_optim_lsq, R2_optim_lsq, P0, Q(t), dt);

figure(1)
subplot(2,1,1)
plot(t,P(t),'DisplayName',strcat('Actual'))
plot(t_FV_lsq, P_FV_lsq, 'DisplayName',strcat('FlowVision lsq'))
plot(t_FV_fmin, P_FV_fmin, 'DisplayName',strcat('FlowVision fmin'))
title('PWK')
xlabel('Time (s)')
ylabel('Pressure (Pa)')
legend show
legend('Location', 'southeast', 'FontSize', 7)

subplot(2,1,2)
plot(t,Q(t), 'DisplayName',strcat('Actual'))
hold on
plot(t,Qsim_sterg,'DisplayName',strcat('Stergiopulos, C = 1.31E-8, R1 = 4.4E+6, R2 = 1.05E+8'))
plot(t,Qsim_lsq,'DisplayName',strcat('Sim LSQ, C =', num2str(C_optim_lsq, '%.3E'), ', R1 = ', num2str(R1_optim_lsq, '%.3E'), ', R2 = ', num2str(R2_optim_lsq, '%.3E')))
plot(t,Qsim_fmin,'DisplayName',strcat('Sim fmin, C =', num2str(C_optim_fmin, '%.3E'), ', R1 = ', num2str(R1_optim_fmin, '%.3E'), ', R2 = ', num2str(R2_optim_fmin, '%.3E')))
plot(t,Qsim_fmin_manual,'DisplayName',strcat('Sim manual, C =', num2str(C, '%.3E'), ', R1 = ', num2str(R1, '%.3E'), ', R2 = ', num2str(R2, '%.3E')))
hold off
title('QWK')
xlabel('Time (s)')
ylabel('Flow (m^{3}/s)')
legend show
legend('Location', 'southeast', 'FontSize', 7)

%% Function for lsq optimisation

function f = fun_lsq(x, time, P, dPdt, Q0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [~,y] = integration(C,R1,R2,time,P,dPdt,Q0);
    Q = y(:,1);
    f = Q;
end

%% Function for fminsearch optimisation

function f = fun_fmin(x, time, Q_in, P, dPdt, Q0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [~,y] = integration(C,R1,R2,time,P,dPdt,Q0);
    Q = y(:,1);
    f = sum((Q-Q_in).^2);
end

%% Function for ODE Integration using equation
function [t,y] = integration(C,R1,R2,time,P,dPdt,Q0)

  ode_fun = @(t,y)((P(t)/(C*R1*R2)) + (dPdt(t)/R1) - (y(1)*(1/C*R1 + 1/C*R2)));
  [t,y] = ode45(ode_fun, time, Q0);
  
end

%% Function for FV PWK

function [t_FV,P_FV] = PWK_FV(C_FV, R1_FV, R2_FV, P_FV_i, Q_FV, timestep, total_time)
    
    t_FV=0:timestep:total_time(end);
    P_FV=zeros(numel(t_FV):1);
    P_FV(1) = P_FV_i;

    for i=2:width(t_FV)

    Beta = timestep / (C_FV*R2_FV);
    Num1 = Q_FV(i)*( (Beta*R2_FV) + (Beta*R1_FV) + R1_FV );
    Num2 = R1_FV*Q_FV(i-1);
    Num3 = P_FV(i-1);
    Num = Num1 - Num2 + Num3;
    Den = 1 + Beta;
    P_FV(i) = Num / Den;

    end
end

I even tried to get rid of the equation for pressure and use tabular data (and reduced the data to 1 cycle), reduced the dt (timestep) to 0.01, as below, but it still gets stuck in the ODE function. Should I try another ode solver than ode45, maybe ode78?

%% Code to Optimize Parameters for Aorta 3-Element Windkessel Model
%%
clear all
close all
clc

%% Inputs

% For mmHg
% R1 * (0.000001/133.322);
% R2 * (0.000001/133.322);
% C  * (133.322/0.000001);

R1 = 4.40e+06;
R2 = 1.05e+08;
C  = 1.31e-08;

% Initial conditions for WK-3 parameters
IC = [C, R1, R2];

Q_peak = 5e-4;  % Peak flow
t_s = 0.28;     % Systolic Time
t_c = 0.84;     % Cycle Time

dt=0.01;      % (CFD/FSI) Simulation Time Step

%% Importing flow and pressure data

RefQWK = readtable('Ref_P_Q_WK_No-Osc_1-Cycle.csv');
RefPWK = readtable('Ref_P_Q_WK_No-Osc_1-Cycle.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);

P0=P0_ref; %Initial P condition
P = smooth(interp1(t_ref, P_ref, t), 50);
dPdt = diff(P)./dt;
dPdt(end+1) = dPdt(end);
P0=P(1);

% 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);

%% Activate for flow data:
% Q_ref=RefQWK.Q;            % multiply by e+6 for mL/s
% Q0_ref=RefQWK.Q(1);
% Q_smooth = smooth(interp1(t_ref, Q_ref, t), 50);
% Q_spline = spline(t,Q_smooth);
% Q_spline_dot = fnder(Q_spline);
% 
% Q=@(t) ppval(Q_spline,t);
% dQdt=@(t) ppval(Q_spline_dot,t);
% 
%Initial Q condition
% Q0 = Q(0);

%Activate for flow equation:
Q_eq = @(t)Q_peak*sin((pi*t)./t_s).^2.*(t>=0).*(t<=t_s);
Q = @(t)Q_eq(mod(t,t_c));
%Initial P condition
Q0 = Q(0);          

%% Curve fitting using the method of least squares

LB = [1e-15 1 1]; UB = [1 1e10 1e10]; % Lower and Upper bounds
options = optimoptions('lsqcurvefit',  'Algorithm', 'levenberg-marquardt', 'MaxFunEvals',10000,'MaxIter',10000,'Display','iter');

f = @(x,t) fun_lsq(x, t, P, dPdt, Q0);
sol = lsqcurvefit(f, IC, t_ref, Q(t), LB, UB, options);

C_optim_lsq = sol(1) 
R1_optim_lsq = sol(2) 
R2_optim_lsq = sol(3)

%% Curve fitting using the fminsearch
 
% options = optimset('MaxFunEvals', 1e100, 'TolX', 1e-16, 'TolFun', 1e-16, 'Display','iter','PlotFcns',@optimplotfval);
options = optimset('MaxFunEvals',10000,'MaxIter',10000, 'Display','iter','PlotFcns',@optimplotfval);

f = @(x) fun_fmin(x, t, Q(t), P, dPdt, Q0);
sol = fminsearch(f,IC,options);

C_optim_fmin = sol(1)
R1_optim_fmin = sol(2)
R2_optim_fmin = sol(3)

%% Solving and Plotting ODE for Q

[t,Qsim_sterg] = integration(1.31e-8,4.4e+6,1.05e+8,t,P,dPdt,Q0);
[t,Qsim_lsq] = integration_data(C_optim_lsq,R1_optim_lsq,R2_optim_lsq,t,P,dPdt,Q0);
[t,Qsim_fmin] = integration_data(C_optim_fmin,R1_optim_fmin,R2_optim_fmin,t,P,dPdt,Q0);
[t,Qsim_fmin_manual] = integration_data(C,R1,R2,t,P,dPdt,Q0);

% PWK for FV
[t_FV_fmin, P_FV_fmin] = PWK_FV(C_optim_fmin, R1_optim_fmin, R2_optim_fmin, P0, Q(t), dt, t);
[t_FV_lsq, P_FV_lsq] = PWK_FV(C_optim_lsq, R1_optim_lsq, R2_optim_lsq, P0, Q(t), dt);

figure(1)
subplot(2,1,1)
plot(t,P(t),'DisplayName',strcat('Actual'))
plot(t_FV_lsq, P_FV_lsq, 'DisplayName',strcat('FlowVision lsq'))
plot(t_FV_fmin, P_FV_fmin, 'DisplayName',strcat('FlowVision fmin'))
title('PWK')
xlabel('Time (s)')
ylabel('Pressure (Pa)')
legend show
legend('Location', 'southeast', 'FontSize', 7)

subplot(2,1,2)
plot(t,Q(t), 'DisplayName',strcat('Actual'))
hold on
plot(t,Qsim_sterg,'DisplayName',strcat('Stergiopulos, C = 1.31E-8, R1 = 4.4E+6, R2 = 1.05E+8'))
% plot(t,Qsim_lsq,'DisplayName',strcat('Sim LSQ, C =', num2str(C_optim_lsq, '%.3E'), ', R1 = ', num2str(R1_optim_lsq, '%.3E'), ', R2 = ', num2str(R2_optim_lsq, '%.3E')))
plot(t,Qsim_fmin,'DisplayName',strcat('Sim fmin, C =', num2str(C_optim_fmin, '%.3E'), ', R1 = ', num2str(R1_optim_fmin, '%.3E'), ', R2 = ', num2str(R2_optim_fmin, '%.3E')))
plot(t,Qsim_fmin_manual,'DisplayName',strcat('Sim manual, C =', num2str(C, '%.3E'), ', R1 = ', num2str(R1, '%.3E'), ', R2 = ', num2str(R2, '%.3E')))
hold off
title('QWK')
xlabel('Time (s)')
ylabel('Flow (m^{3}/s)')
legend show
legend('Location', 'southeast', 'FontSize', 7)

%% Function for lsq optimisation

function f = fun_lsq(x, time, P, dPdt, Q0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [~,y] = integration(C,R1,R2,time,P,dPdt,Q0);
    Q = y(:,1);
    f = Q;
end

%% Function for fminsearch optimisation

function f = fun_fmin(x, time, Q_in, P, dPdt, Q0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [~,y] = integration(C,R1,R2,time,P,dPdt,Q0);
    Q = y(:,1);
    f = sum((Q-Q_in).^2);
end

%% Function for ODE Integration using equation
function [t,y] = integration(C,R1,R2,time,P,dPdt,Q0)

  % For P data
  ode_fun = @(t,y) Q_odefcn(t, y, C, R1, R2, time, P, dPdt);

  % For P equation:
  % ode_fun = @(t,y)((P(t)/(C*R1*R2)) + (dPdt(t)/R1) - ((y(1)*(1 + R1/R2))/C*R1));
  
  [t,y] = ode45(ode_fun, time, Q0);
  
end

%% Function for ODE using flow tabular data
function dQdt = Q_odefcn(t, y, C, R1, R2, time, P_in, dPdt_in)
% dQdt = zeros(1,1); %Initializing the output vector
% Determine the closest point to the current time
%     [d, closestIndex]=min(abs(time-t));
%     P = interp1(time,P_in,t);
%     for i = 1:numel(time)
%       if time(i+1) >= t
%         dPdt = (P_in(i+1)-P_in(i))/(time(i+1)-time(i));
%         break
%       end
%     end
    %dQ = dQ_in(closestIndex);

   dQdt = zeros(1,1); %Initializing the output vector

%Determine the closest point to the current time
    [~, closestIndex]=min(abs(time-t));
    P = P_in(closestIndex);
    dPdt = dPdt_in(closestIndex);
  
   dQdt(1) = P/(C*R1*R2) + (dPdt/R1) - y(1)*(1/C*R1 + 1/C*R2);

end

%% Function for FV PWK

function [t_FV,P_FV] = PWK_FV(C_FV, R1_FV, R2_FV, P_FV_i, Q_FV, timestep, total_time)
    
    t_FV=0:timestep:total_time(end);
    P_FV=zeros(numel(t_FV):1);
    P_FV(1) = P_FV_i;

    for i=2:width(t_FV)

    Beta = timestep / (C_FV*R2_FV);
    Num1 = Q_FV(i)*( (Beta*R2_FV) + (Beta*R1_FV) + R1_FV );
    Num2 = R1_FV*Q_FV(i-1);
    Num3 = P_FV(i-1);
    Num = Num1 - Num2 + Num3;
    Den = 1 + Beta;
    P_FV(i) = Num / Den;

    end
end

[Matt J]

I think it is the fnder algorithm for differentiation that is really slowing it down.

Why do you think that? Did you profile it?

[Hussam]

Because everytime I stop it, it is always stuck somewhere here:

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);

This is the command window:

Operation terminated by user during ppval

In QWK_3E_I>@(t)ppval(P_spline,t) (line 42)

In QWK_3E_I>@(t,y)((P(t)/(C*R1*R2))+(dPdt(t)/R1)-((y(1)*(1+R1/R2))/C*R1)) (line 149)

In ode45 (line 308)

f6 = ode(t6, y6);

In QWK_3E_I>integration (line 150)

% Q = interp1(time,Q_in,t);

In QWK_3E_I>fun_lsq (line 131)

C = x(1); R1 = x(2); R2 = x(3);

In QWK_3E_I>@(x,t)fun_lsq(x,t,P,dPdt,Q0) (line 72)

In lsqcurvefit (line 251)

initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

In QWK_3E_I (line 73)

f = @(x,t) fun_lsq(x, t, P, dPdt, Q0);

I will rerun again with profiler open

[Hussam]

Looks like it is ppval. What would you recommend I change it to? I really would prefer to use an equation since the derivative is much clear/smoother. Thoughts?

[Catalytic]

Likely "clear all" is not helping. From the doc: "Calling clear all, clear classes, and clear functions decreases code performance, and is usually unnecessary."

[Hussam]

@Catalytic Interesting, did not know that, but I only call it in the beginning to clear everything for a fresh start, it is not called again, I think. Would it still be slowing down the code??

[Hussam]

I think I solved this problem by replacing ode45 with a stiff ode15s solver. But the results are bad (no optimization)? Please see below.

%% Code to Optimize Parameters for Aorta 3-Element Windkessel Model
%%
clear all
close all
clc

%% Inputs

% For mmHg
% R1 * (0.000001/133.322);
% R2 * (0.000001/133.322);
% C  * (133.322/0.000001);

R1 = 4.40e+06;
R2 = 1.05e+08;
C  = 1.31e-08;

% Initial conditions for WK-3 parameters
IC = [C, R1, R2];

Q_peak = 5e-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=P_ref;
% P0=P0_ref; %Initial P condition
% dPdt = 0; %dummy variable

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);

%% Activate for flow data:
% Q_ref=RefQWK.Q;            % multiply by e+6 for mL/s
% Q0_ref=RefQWK.Q(1);
% Q_smooth = smooth(interp1(t_ref, Q_ref, t), 50);
% Q_spline = spline(t,Q_smooth);
% Q_spline_dot = fnder(Q_spline);
% 
% Q=@(t) ppval(Q_spline,t);
% dQdt=@(t) ppval(Q_spline_dot,t);
% 
%Initial Q condition
% Q0 = Q(0);

%Activate for flow equation:
Q_eq = @(t)Q_peak*sin((pi*t)./t_s).^2.*(t>=0).*(t<=t_s);
Q = @(t)Q_eq(mod(t,t_c));
%Initial P condition
Q0 = Q(0);          

%% Curve fitting using the method of least squares

LB = [1e-15 1 1]; UB = [1 1e10 1e10]; % Lower and Upper bounds
options = optimoptions('lsqcurvefit',  'Algorithm', 'levenberg-marquardt', 'MaxFunEvals',10000,'MaxIter',10000,'Display','iter');

f = @(x,t) fun_lsq(x, t, P, dPdt, Q0);
sol = lsqcurvefit(f, IC, t, Q(t), LB, UB, options);

C_optim_lsq = sol(1) 
R1_optim_lsq = sol(2) 
R2_optim_lsq = sol(3)

%% Curve fitting using the fminsearch
 
options = optimset('MaxFunEvals', 1e100, 'TolX', 1e-16, 'TolFun', 1e-16, 'Display','iter','PlotFcns',@optimplotfval);
%options = optimset('MaxFunEvals',10000,'MaxIter',10000, 'Display','iter','PlotFcns',@optimplotfval);

f = @(x) fun_fmin(x, t, Q, P, dPdt, Q0);
sol = fminsearch(f,IC,options);

C_optim_fmin = sol(1)
R1_optim_fmin = sol(2)
R2_optim_fmin = sol(3)

%% Solving and Plotting ODE for Q

[t,Qsim_sterg] = integration(1.31e-8,4.4e+6,1.05e+8,t,P,dPdt,Q0);
[t,Qsim_lsq] = integration(C_optim_lsq,R1_optim_lsq,R2_optim_lsq,t,P,dPdt,Q0);
[t,Qsim_fmin] = integration(C_optim_fmin,R1_optim_fmin,R2_optim_fmin,t,P,dPdt,Q0);
[t,Qsim_fmin_manual] = integration(C,R1,R2,t,P,dPdt,Q0);

% PWK for FV
[t_FV_fmin, P_FV_fmin] = PWK_FV(C_optim_fmin, R1_optim_fmin, R2_optim_fmin, P0, Q(t), dt, t);
[t_FV_lsq, P_FV_lsq] = PWK_FV(C_optim_lsq, R1_optim_lsq, R2_optim_lsq, P0, Q(t), dt, t);

figure(1)
subplot(2,1,1)
plot(t,P(t),'DisplayName',strcat('Actual'))
plot(t_FV_lsq, P_FV_lsq, 'DisplayName',strcat('FlowVision lsq'))
plot(t_FV_fmin, P_FV_fmin, 'DisplayName',strcat('FlowVision fmin'))
title('PWK')
xlabel('Time (s)')
ylabel('Pressure (Pa)')
legend show
legend('Location', 'southeast', 'FontSize', 7)

subplot(2,1,2)
plot(t,Q(t), 'DisplayName',strcat('Actual'))
hold on
plot(t,Qsim_sterg,'DisplayName',strcat('Stergiopulos, C = 1.31E-8, R1 = 4.4E+6, R2 = 1.05E+8'))
plot(t,Qsim_lsq,'DisplayName',strcat('Sim LSQ, C =', num2str(C_optim_lsq, '%.3E'), ', R1 = ', num2str(R1_optim_lsq, '%.3E'), ', R2 = ', num2str(R2_optim_lsq, '%.3E')))
plot(t,Qsim_fmin,'DisplayName',strcat('Sim fmin, C =', num2str(C_optim_fmin, '%.3E'), ', R1 = ', num2str(R1_optim_fmin, '%.3E'), ', R2 = ', num2str(R2_optim_fmin, '%.3E')))
plot(t,Qsim_fmin_manual,'DisplayName',strcat('Sim manual, C =', num2str(C, '%.3E'), ', R1 = ', num2str(R1, '%.3E'), ', R2 = ', num2str(R2, '%.3E')))
hold off
title('QWK')
xlabel('Time (s)')
ylabel('Flow (m^{3}/s)')
legend show
legend('Location', 'southeast', 'FontSize', 7)

%% Function for lsq optimisation

function f = fun_lsq(x, time, P, dPdt, Q0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [~,y] = integration(C,R1,R2,time,P,dPdt,Q0);
    Q = y(:,1);
    f = Q';
end

%% Function for fminsearch optimisation

function f = fun_fmin(x, time, Q_in, P, dPdt, Q0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [t,y] = integration(C,R1,R2,time,P,dPdt,Q0);
    Q = y(:,1);
    Q_in = Q_in(t);
    f = sum((Q-Q_in).^2);
end

%% Function for ODE Integration using equation
function [t,y] = integration(C,R1,R2,time,P,dPdt,Q0)

  % For P data
  % ode_fun = @(t,y) Q_odefcn(t, y, C, R1, R2, time, P);

  % For P equation:
  ode_fun = @(t,y)((P(t)/(C*R1*R2)) + (dPdt(t)/R1) - (y(1)*(1/(C*R1) + 1/(C*R2))));
  
  [t,y] = ode15s(ode_fun, time, Q0);
  
end

%% Function for ODE using flow tabular data
function dQdt = Q_odefcn(t, y, C, R1, R2, time, P_in)
dQdt = zeros(1,1); %Initializing the output vector
%Determine the closest point to the current time
    %[d, closestIndex]=min(abs(time-t));
    P = interp1(time,P_in,t);
    for i = 1:numel(time)
      if time(i+1) >= t
        dPdt = (P_in(i+1)-P_in(i))/(time(i+1)-time(i));
        break
      end
    end
    %dQ = dQ_in(closestIndex);
   dQdt(1) = P/(C*R1*R2) + (dPdt/R1) - y(1)*(1/C*R1 + 1/C*R2);

end

%% Function for FV PWK

function [t_FV,P_FV] = PWK_FV(C_FV, R1_FV, R2_FV, P_FV_i, Q_FV, timestep, total_time)
    
    t_FV=0:timestep:total_time(end);
    P_FV=zeros(numel(t_FV):1);
    P_FV(1) = P_FV_i;

    for i=2:width(t_FV)

    Beta = timestep / (C_FV*R2_FV);
    Num1 = Q_FV(i)*( (Beta*R2_FV) + (Beta*R1_FV) + R1_FV );
    Num2 = R1_FV*Q_FV(i-1);
    Num3 = P_FV(i-1);
    Num = Num1 - Num2 + Num3;
    Den = 1 + Beta;
    P_FV(i) = Num / Den;

    end
end

[Torsten]

Your ode function is wrong.

It should read

ode_fun = @(t,y)((P(t)/(C*R1*R2)) + (dPdt(t)/R1) - (y(1)*(1/(C*R1) + 1/(C*R2))));

instead of

ode_fun = @(t,y)((P(t)/(C*R1*R2)) + (dPdt(t)/R1) - (y(1)*(1/C*R1 + 1/C*R2)));

[Hussam]

@Torsten Hi, thank you for that. I did realize that and fixed it, but still it does not work with ode45, and ode15s is not working properly I assume (as shown in the comments).

[Hussam]

Here are the results for the ode15s - very bad.

[Torsten]

Don't you see that the actual pressure curves always asymptotically converge towards atmospheric pressure while you still set pout =0 in your ode function ?

[Hussam]

@Torsten Yes, I do see that, even though mathematically it does make sense to set pout=p0, I think physiologically according to the windkessel model, the pressure cycle should go back to atmospheric (to close the loop correctly). If I try to implement a WK model that has pout=p0 in my CFD/FSI simulations, the pressure diverges (increases to infinity) as the flow changes and is depended on the pressure and the windkessel model keeps on adding p0 to the pressure at every timestep when calculating the pressure. Does this answer your question?

[Torsten]

Does this answer your question?

No. Your model variable "pout" has to be set to the asymptotic value of your measurement data for your parameter identification to make sense. Why don't you do that ?

[Hussam]

Hi @Torsten, I did try it, and it matches the diastolic pressure very well, however it causes my FSI and CFD simulations to diverge - I am assuming it is because the flow input is not fixed as it is in the matlab code. Ideas?

[Torsten]

I thought you try to reproduce the pressure or flow curves you supply in the .csv files. Here, pressure converges asymptotically not against 0, but against atmospheric pressure, and that's what your model equations should reflect. Or is there an iterative procedure between MATLAB and CFD ?

In CFD calculations, you often solve the equations for pressure above atmospheric pressure p'. The "real" pressure is then p = p' + p_atmospheric. Maybe you set p' = p_atmospheric when the solver diverges, thus an overpressure of 1.013 bar. But it's only speculation.

[Hussam]

So, I'm trying to target to output pressure to be the same as the one provided to the matlab code, but I am unaware of the flow. So I am assuming a certain flow (with the flow equation) for the code and optimize the WK parameters for that. Then I run the FSI simulation with WK at the outlet boundary condition and I end up with a different flow curve. I then take that flow curve and input it to the code, and iterate. However, since they are aortic valve FSI simulations, they take quite a long time, so I am trying to figure out a faster way to converge to the right parameters that I want. Any change in the WK parameters change my flow and cardiac output (so far the peak flow ranges from 5e-4 to 1e-3 m^3/s). Ideally, I would like to target the flow equation with a peak value of 5e-4 m^3/s.

In the FSI simulation, when I set p_out to the diastolic pressure (same as in the code), the pressure diverges upwards, keeps on increasing to infinity. It works well when I set p_out = 0. I understand that it should not be 0 (atmospheric), but it should not be diastolic pressure either (at least physiologically), possibly some small value. I have seen a lot of people in the community/literature use it as zero.

I hope this shows my thinking more. Please correct me if I am wrong. What do you think?