Optimization Problem for 3-Element Windkessel Model. Need help with improving my optimization and looking for better optimization alogrithims.

Hi everyone, so I am trying to optimize a standard 3-element windkessel model for the aorta based on pressure and flow data (attached). I have tested the least square method (LSQ) and fminsearch. I think LSQ does not do a good job at all, it does not change any of my parameters, and even though fminsearch is better, it stil has a lot of variance/error. Kindly advise. The circuit and equations are as below, and I have also added my code and results for fminsearch. Cheers!

clear all
close all
clc

%% Inputs
R1 = 9.3896e+6;
R2 = 2.4023e+7;
C = 8.0187e-9;
%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.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-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, 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_sterg] = integration(1.31e-8,4.4e+6,1.05e+8,t,Q,I0,Ts,IC_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);
% [t,Psim_fmin_manual] = integration(8e-9,6e+6,2.4023e+7,t,Q,I0,Ts,IC_P);

plot(t,P)
hold on
% plot(t,Psim_sterg)
% plot(t,Psim_lsq)
plot(t,Psim_fmin)
% plot(t,Psim_fmin_manual)
hold off

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

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

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

  %input current flow for 4 cycle
  I=@(t)I0*sin((pi*t)./Ts).^2.*(t<=Ts | 0.84<=t & t<=0.84+Ts | 2*0.84<=t & t<=2*0.84+Ts | 3*0.84<=t & t<=3*0.84+Ts); %input current flow
  Idot = @(t)I0*2*sin(pi*t./Ts).*cos(pi*t./Ts)*pi./Ts.*(t<=Ts | 0.84<=t & t<=0.84+Ts | 2*0.84<=t & t<=2*0.84+Ts | 3*0.84<=t & t<=3*0.84+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

[Sam Chak]

Hi @Hussam

Just out of curiosity, I'm wondering how you managed to express this elegant two-state coupled differential equation of Windkessel into a concise single-line equation as described in 'fun2'?

It is crucial to acknowledge that if the model is inaccurately described from the outset, any optimization efforts invested may end up being futile.

Additionally, it is vital to recognize that the mathematical notation of I as a function of time, , does not directly translate to 'I(t)' in MATLAB code. For example, 'y(1)' in the code does not mean the value of at .

I'm uncertain about the best way to explain this, as it may require the expertise of an MVP-level user to provide clarification.

[Hussam]

Hi @Sam Chak

Thank you for your feedback. I am not very sure I am following. What is the problem with the ODE function 'fun2'?

Also, if I(t) is not the same as it is mathematically, then what is it? I understand that y(1) in my code represents y(t) in general, for all t.

May you please make some changes to my code so that I can understand better?

Cheers!

[Sam Chak]

My bad, @Hussam. I overlooked the anonymous functions that declared and . I'm unfamilar with the 3-Element Windkessel Model. I believe is although you didn't explicitly mention that. Isn't time series data Q(t) provided in the csv file?

[Hussam]

No worries. Oh yes, sorry I didn't make that clear. Yes, I was using tabular Q(t) data before, but was recommended to move to an equation for flow instead, I(t). I think the optimization results were better when I did that.

Answer 1

[Sam Chak]

Hi @Hussam

In my solution, I suggested utilizing the dI/dt-free model:

This approach enables me to directly utilize the raw data points {t, , and } from the CSV file, without the need to inaccurately estimate . Here is the proof demonstrating that both models are equivalent. In fact, I employed basic differentiation skills and a change of coordinates to achieve this.

Rearranging the equation to get

I hope that this alternative model offers you a viable approach to identify the parameters in the 3-element Windkessel Model.

Answer 2

[Sam Chak]

Hi @Hussam

I also haven't been able to obtain a better identified system for the 3-element Windkessel Model. However, it would be beneficial for you to review if my modeling approach aligns with your understanding. Sometimes, when the results don't match as expected, it might indicate the need to fit the data using the 4-element Windkessel Model.

%% Extract data from csv
RefPQWK = readtable('Ref_P_Q_WK_No-Osc.csv');
tout    = RefPQWK.t;
Iout    = RefPQWK.Q;
Pout    = RefPQWK.P;

%% Identified parameters
C       = 1.22714690629244e-08;
R1      = 10203268.9540008; 
R2      = 15476330.6281995;

%% Call ode45 to solve identified 3-element Windkessel Model
tspan   = [tout(1), tout(end)];
P0      = Pout(1);
I0      = Iout(1);
x0      = P0 - R1*I0;
fun     = @(t, x) odefcn(t, x, C, R1, R2, P0, tout, Iout);
sol     = ode45(fun, tspan, x0);
x       = deval(sol, tout);
P       = 1*x' + R1*Iout;   % Output of 3-element Windkessel Model

%% Plot result
plot(tout, P),    hold on   % Simulated Pressure
plot(tout, Pout), grid on   % Measured  Pressure
title('3-element Windkessel Model')
xlabel('Time (s)')
ylabel('Pressure (Pa)')
legend('3eWM Response', 'Measured Data')

%% 3-element Windkessel Model
function dx = odefcn(t, x, C, R1, R2, P0, tout, Iout)
    % Parameters
    a1  = 1/(C*R2);
    b0  = R1;
    b1  = (1/C)*(1 + R1/R2);
    A   = - a1;                     % state coefficient
    B   = b1 - a1*b0;               % input coefficient
    
    % Interpolate the data set (tout, Iout) at time t
    I   = interp1(tout, Iout, t);
    
    % state equation of 3-element Windkessel Model
    dx  = A*(x - P0) +  B*I;
end

[Sam Chak]

Hi @Hussam

By the way, I used the following state-space model to represent the 3-element Windkessel Model:

State equation

Output equation

[Hussam]

Ah great, thank you for this! I'm not sure, but is there any optimization taking place? Might have missed it. Thanks!

[Sam Chak]

Hi @Hussam

The system identification process was conducted offscreen, and while the optimization methods are important, the crucial factors are the mathematical model and the chosen cost function for the optimization process.

If the appropriate Windkessel model is employed, any effective optimization method will yield more or less the same result because the computed cost functions will be identical.

One key distinction to note is that my approach to modeling the 3-element Windkessel effect does NOT involve differentiating the flow signal, . This avoids amplifying any noise that may be present in the flow signal.

[Hussam]

Hi @Sam Chak, thank you for your answer. I think I might be a little lost with the math itself, how did you substitute or replace the differentiation of I(t)? Thanks.