Optimisation for 3-element Windkessel Model Not Working - Using Least Squares Method and fminsearch. Advice? Should I use other optimisation algorithms?

So I am trying to optimise the parameters (R1, R2, C) for the 3-element windkessel model for the aorta. I am supplying it with tabular data for the Pressure and Flow (Q), but neither of the models are doing anything, and not giving any errors either. Kindly advise.

My code is as follows:

clear all
close all
clc

%% Inputs
R1 = 4.4e+6;
R2 = 1.05e+8;
C = 1.31e-8;
% L = 6.799e+5;
TimeStep = 0.001;
Beta = TimeStep / (C*R2);

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

%% Importing flow and pressure data
RefPQWK = readtable('Ref_P_Q_WK_No-Osc_1-Cycle.csv');

t=RefPQWK.t;
Q=RefPQWK.Q;
P=RefPQWK.P;
P0=RefPQWK.P(1);
Q0=RefPQWK.Q(1);
dt = TimeStep;

t_0 = 0:dt:t(end);
Q = smooth(interp1(t, Q, t_0), 50);
P = smooth(interp1(t, P, t_0), 50);

%Initial conditions
IC_Q = Q(1);
IC_P = P(1);

dP = diff(P)./dt;
dP(end+1) = dP(end);

dQ = diff(Q)./dt;
dQ(end+1) = dQ(end);

%dqdt = derivative(Q,t_0);

Unrecognized function or variable 'derivative'.

%% Curve fitting using the method of least squares

LB = [1e-10 1 1]; UB = [1 1e10 1e10]; % Lower and Upper bounds

options = optimoptions('lsqcurvefit',  'Algorithm', 'levenberg-marquardt', 'StepTolerance', 1e-100,'FunctionTolerance', 1e-10,'Display','iter');
%options = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations', 1e100, 'MaxIterations', 1e4, 'StepTolerance',1e-300, 'Display','iter');
%options =optimoptions('lsqcurvefit','MaxFunctionEvaluations', 1e100, 'MaxIterations', 1e100,'Display','iter');

f = @(x, t_0)fun_lsq(x, t_0, dt, P, Q, dQ, IC_P);

 %ydata=zeros(1, height(P))';
 sol = lsqcurvefit(f, IC, t_0, P, LB, UB, options);
 C_optim = sol(1); Rp_optim = sol(2); Rd_optim = sol(3);

%% Curve fitting using the fminsearch
% options = optimset('Display','iter','PlotFcns',@optimplotfval);
% 
% f = @(x)fun_fmin(x, t_0, dt, P, Q, dQ, IC_P);
% sol = fminsearch(f,IC);

%% Solving and Plotting ODE for P

tspan = [t_0(1):dt:t_0(end)]; %Time interval of the integration

funP = @(t,y_P) P_odefcn(t, y_P, C, R1, R2, t_0, Q, dQ);
[t,y_P] = ode78(funP, tspan, IC_P);

funP_opt = @(t,y_P_opt) P_odefcn(t, y_P_opt, C, R1, R2, t_0, Q, dQ);
[t,y_P_opt] = ode78(funP, tspan, IC_P);

% y_Q = smooth(interp1(t, y_Q, t_0), 50);
y_P = smooth(interp1(t, y_P, t_0), 50);
y_P_opt = smooth(interp1(t, y_P, t_0), 50);

figure(1), plot(t_0, P), hold on, plot(t_0, y_P); hold on, plot(t_0, y_P_opt); 
xlabel('Time'); ylabel('P')
legend('Actual', 'Stergiopoulos', 'Optimised')

%% Function for lsq optimisation

function f = fun_lsq(x, time, dt, P_in, Q_out, dQ_out, IC_p) 

    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);

    tspan = [0:dt:time(end)]; %Time interval of the integration

    %Solving the ODEs
    fun2 = @(t,y) P_odefcn(t, y, C, R1, R2, time, Q_out, dQ_out);
    [t,y] = ode78(fun2, tspan, IC_p);

    P = y(:,1);

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

    f = P-P_in;

end

%% Function for fminsearch optimisation

function f = fun_fmin(x, time, dt, P_in, Q_out, dQ_out, IC_p) 

    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);

    tspan = [0:dt:time(end)]; %Time interval of the integration

    %Solving the ODEs
    fun2 = @(t,y) P_odefcn(t, y, C, R1, R2, time, Q_out, dQ_out);
    [t,y] = ode78(fun2, tspan, IC_p);

    P = y(:,1);

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

    f = sqrt(sum((P-P_in).^2));

end

%% Function for ODE

function dPdt = P_odefcn(t, y, C, Rp, Rd, time, Q_in, dQ_in)

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

%Determine the closest point to the current time
    [d, closestIndex]=min(abs(time-t));
    Q = Q_in(closestIndex);
    dQ = dQ_in(closestIndex);

   dPdt(1) =(Q/C)*(1+Rp/Rd) + Rp*dQ - y(1)/(C*Rd);
      
end

Both results look like this:

I do get the following for the LSQ method though:

[Torsten]

We miss your .csv file.

[Hussam]

Attached it to the file. So they are results from FSI simulations with oscillations. I manually cleaned it. Both files are attached. I also reduced it to 1 and 2 cycles, attached to this comment. Thanks!

[Torsten]

You didn't include the function "derivative". Please test beforehand under MATLAB online if the code works.

[Hussam]

Hi, I'm not using the derivative function (I was before). I have commented it. The code should be working. Thank you.

[Torsten]

If you use "lsqcurvefit", you must return

f = P;

not

f = P-P_in;

Don't you have a separate equation for Q so that you wouldn't need to differentiate your measurement data ?

[d, closestIndex]=min(abs(time-t));
Q = Q_in(closestIndex);
dQ = dQ_in(closestIndex);

[Hussam]

Hi, thank you for your feedback.

Fixed this, but still getting the same result.

f = P-P_in;

to

f = P;

Unfortunately, I do not have an equation for Q, since it is the output of a simulation, I can try to formulate one though - would that make a difference?

[d, closestIndex]=min(abs(time-t));
Q = Q_in(closestIndex);
dQ = dQ_in(closestIndex);

Actually, I think this piece of code above is redundant since I smooth out my data according to dt.

Answer 1

[Torsten]

I'd try to continue with the code below.

It seems that a variation of your initial parameters does not cause a change in P such that the initial values for C, R1 and R2 are returned by the solver.

Unfortunately, I do not have an equation for Q, since it is the output of a simulation, I can try to formulate one though - would that make a difference?

The difference is that the inputs to the differential equation for P would become smooth which is usually necessary for the integrators to succeed.

%% Inputs
R1 = 4.4e+6;
R2 = 1.05e+8;
C = 1.31e-8;

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

%% Importing flow and pressure data
RefPQWK = readtable('Ref_P_Q_WK_No-Osc_1-Cycle.csv');

t=RefPQWK.t;
Q=RefPQWK.Q;
P=RefPQWK.P;
P0=RefPQWK.P(1);
Q0=RefPQWK.Q(1);

%Initial conditions
IC_P = P(1);

%% Curve fitting using the method of least squares

LB = [1e-10 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, Q, IC_P);

sol = lsqcurvefit(f, IC, t, P, LB, UB, options);

                                        First-order                     Norm of 
 Iteration  Func-count      Resnorm      optimality       Lambda           step
     0           4      2.18762e+10          0.0747         0.01
     1          18                                                     32.8184

Local minimum possible.
lsqcurvefit stopped because the relative size of the current step is less than
the value of the step size tolerance.

C_optim = sol(1) 

C_optim = 1.3100e-08

Rp_optim = sol(2) 

Rp_optim = 4400000

Rd_optim = sol(3)

Rd_optim = 105000000

%% Solving and Plotting ODE for P
[t,Psim] = integration(C_optim,Rp_optim,Rd_optim,t,Q,IC_P);

hold on
plot(t,P)
plot(t,Psim)
hold off

%% Function for lsq optimisation

function f = fun_lsq(x, time, P_in, Q_out, IC_p) 

    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);

    [t,y] = integration(C,R1,R2,time,Q_out,IC_p);

    P = y(:,1);

    f = P;

end

function [t,y] = integration(C,R1,R2,time,Q_out,IC_p)

  fun2 = @(t,y) P_odefcn(t, y, C, R1, R2, time, Q_out);
  [t,y] = ode45(fun2, time, IC_p);

end

%% Function for ODE

function dPdt = P_odefcn(t, y, C, Rp, Rd, time, Q_in)

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

%Determine the closest point to the current time
    %[d, closestIndex]=min(abs(time-t));
    Q = interp1(time,Q_in,t);
    for i = 1:numel(time)
      if time(i+1) >= t
        dQ = (Q_in(i+1)-Q_in(i))/(time(i+1)-time(i));
        break
      end
    end
    %dQ = dQ_in(closestIndex);

   dPdt(1) =(Q/C)*(1+Rp/Rd) + Rp*dQ - y(1)/(C*Rd);
      
end

[Hussam]

Thank you for your feedback and modifications. I will continue with that. I will also try to generate an equation for Q. Out of curiosity, do you think it would make a different to optimize against Q instead of P?

Also, what do you think of using other algorithm methods, such as fminsearch, or such?

Cheers!

[Torsten]

According to your curve for P, P returns to the initial pressure P(1) as Q and dQ become zero. But your differential equation for P indicates that as Q and dQ become 0, P goes to 0 as you set p_out = 0. Is that correctly set by you ?

[Hussam]

Well, if I understand correctly, the pressure p_out (that is grounded after the 3 elements) should be zero, and therefore I have set it to zero. Please correct me if I am wrong. However, I might not have understood correctly. The pressure at the beginning of the circuit should return to P(1), as shown on the graph. That is when Q starts to increase again. Does this answer your question?

[Torsten]

When Q and dQ are zero in your measurement curves, your measurement pressure curve returns to the initial pressure P at time t = 0.

When Q and dQ are zero in your model equation, your simulated pressure will go to p_out which is set to 0.

So I think your model pressure and your measured pressure cannot be compared. Either p_out must be set to the measured pressure at time t = 0 or your P from the model equation and your measured P must be taken at the same locus.

[Hussam]

Also, I think I can approximate Q as follows:

Ts = 0.28;
t=0:0.0001:0.84;
Q0=3.5e-4;
I=Q0*sin((pi*t)/Ts).^2.*(t<=Ts);
Q=@(t)Q0*sin((pi*t)/Ts).^2.*(t<=Ts);

plot(t,I);

[Hussam]

So, I've played around with P_out, setting it to zero and to P(1) and I got the following curves. I also added fminsearch optimisation, which seems to be better (added code at the end). It seems like P_out=0 is capturing the pressure better. I think your 1st order discretization of the differentiation is better. I will implement the sin equation instead of the Q data and investigate further - will keep you updated! Thank you very much for you help! Would you recommend any other optimization alogrithims to try?

Seems to be getting worse with multiple cycles though:

%% Function for fminsearch optimisation
function f = fun_fmin(x, time, P_in, Q_out, IC_p) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [t,y] = integration(C,R1,R2,time,Q_out,IC_p);
    P = y(:,1);
    f = sqrt(sum((P-P_in).^2));
end

%% Curve fitting using the fminsearch

% % options = optimset('Display','iter','PlotFcns',@optimplotfval);
% % options = optimset('TolX', 1e-16, 'TolFun', 1e-16, 'MaxFunEvals', 1e100, 'MaxIter', 1e100, 'Display','iter','PlotFcns',@optimplotfval);
% 
options = optimset('Display','iter','PlotFcns',@optimplotfval);
f = @(x)fun_fmin(x, t, P, Q, IC_P);
sol = fminsearch(f,IC,options);

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

[Torsten]

In my opinion, setting p_out = 0 is simply wrong because measurements and model have different asymptotic behavior (measurements tend to P(1), model tends to 0).

In your objective function for fminsearch, you should replace

f = sqrt(sum((P-P_in).^2))

by

f = sum((P-P_in).^2);

The sqrt makes your objective function non-differentiable at 0.

[Hussam]

I think I agree with you regarding the p_out. I have implemented the sine function for Q and below are my results for 1 cycle. fminsearch is doing a better job so far. How can I improve this further? (My code is attached below).

clear all
close all
clc

%% Inputs
R1 = 4.4e+6;
R2 = 1.05e+8;
C = 1.31e-8;
%Initial conditions for WK-3 parameters
IC = [C, R1, R2];

I0 = 3.5e-4;
Ts = 0.28;

%% Importing flow and pressure data
RefPQWK = readtable('Ref_P_Q_WK_No-Osc_1-Cycle.csv');
t=RefPQWK.t;
dt=0.0001;
Q=RefPQWK.Q;
P=RefPQWK.P;
P0=RefPQWK.P(1);
Q0=RefPQWK.Q(1);
%Initial conditions
IC_P = P(1);

%% Curve fitting using the method of least squares
LB = [1e-10 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, Q, I0, Ts, IC_P);
sol = lsqcurvefit(f, IC, t, P, 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('Display','iter','PlotFcns',@optimplotfval);
% % options = optimset('TolX', 1e-16, 'TolFun', 1e-16, 'MaxFunEvals', 1e100, 'MaxIter', 1e100, 'Display','iter','PlotFcns',@optimplotfval);
% 
options = optimset('Display','iter','PlotFcns',@optimplotfval);
f = @(x)fun_fmin(x, t, P, Q, I0, Ts, IC_P);
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 P
[t,Psim_lsq] = integration(C_optim_lsq,R1_optim_lsq,R2_optim_lsq,t,Q,I0,Ts,IC_P);
[t,Psim_fmin] = integration(C_optim_fmin,R1_optim_fmin,R2_optim_fmin,t,Q,I0,Ts,IC_P);

plot(t,P)
hold on
plot(t,Psim_lsq)
plot(t,Psim_fmin)
hold off

title('PWK')
xlabel('Time (s)')
ylabel('Pressure (Pa)')
legend('Actual','Sim LSQ', 'Sim fmin')

%% Function for lsq optimisation
function f = fun_lsq(x, time, P_in, Q_out, I0, Ts, IC_p) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [t,y] = integration(C,R1,R2,time,Q_out, I0, Ts, IC_p);
    P = y(:,1);
    f = P;
end

%% Function for fminsearch optimisation
function f = fun_fmin(x, time, P_in, Q_out, I0, Ts, IC_p) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [t,y] = integration(C,R1,R2,time,Q_out,I0,Ts,IC_p);
    P = y(:,1);
    f = sum((P-P_in).^2);
end


%% Function for ODE Integration using flow tabular data
% function [t,y] = integration(C,R1,R2,time,Q_out,IC_p)
%   fun2 = @(t,y) P_odefcn(t, y, C, R1, R2, time, Q_out, IC_p);
%   [t,y] = ode45(fun2, time, IC_p);
% end

%% Function for ODE using flow tabular data
% function dPdt = P_odefcn(t, y, C, Rp, Rd, time, Q_in,p0)
% dPdt = zeros(1,1); %Initializing the output vector
% %Determine the closest point to the current time
%     %[d, closestIndex]=min(abs(time-t));
%     Q = interp1(time,Q_in,t);
%     for i = 1:numel(time)
%       if time(i+1) >= t
%         dQ = (Q_in(i+1)-Q_in(i))/(time(i+1)-time(i));
%         break
%       end
%     end
%     %dQ = dQ_in(closestIndex);
%    dPdt(1) =(Q/C)*(1+Rp/Rd) + Rp*dQ - (y(1)-0)/(C*Rd);
% 
% end

%% Function for ODE Integration using flow equation
function [t,y] = integration(C,R1,R2,time,Q_out,I0, Ts, IC_p)

  I = @(t)I0*sin((pi*t)./Ts).^2.*(t<=Ts); %input current flow
  Idot = @(t)I0*2*sin(pi*t./Ts).*cos(pi*t./Ts)*pi./Ts.*(t<=Ts);

  fun2 = @(t,y)[(I(t)/C)*(1+R1/R2) + R1*Idot(t) - (y(1)-IC_p)/(C*R2)];

  % tspan= [0 time(end)];
  [t,y] = ode45(fun2, time, IC_p);
end