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

Question

Hussam 2024-3-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098386-optimization-problem-for-3-element-windkessel-model-need-help-with-improving-my-optimization-and-lo

回答： Sam Chak 2024-3-28，15:04

Ref_P_Q_WK_No-Osc.csv

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Sam Chak 2024-3-25

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 2024-3-25

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 2024-3-28，15:04

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098386-optimization-problem-for-3-element-windkessel-model-need-help-with-improving-my-optimization-and-lo#answer_1432971

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.

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

请先登录，再进行评论。

Answer 2

Sam Chak 2024-3-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098386-optimization-problem-for-3-element-windkessel-model-need-help-with-improving-my-optimization-and-lo#answer_1430691

在 MATLAB Online 中打开

Ref_P_Q_WK_No-Osc.csv

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Sam Chak 2024-3-25

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 2024-3-28，13:56

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.

请先登录，再进行评论。