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

25 次查看（过去 30 天）
Hussam
Hussam 2024-3-25
回答： Rykelle 2024-10-24
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更早的评论
Sam Chak
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
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.

请先登录，再进行评论。

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

采纳的回答

Sam Chak
Sam Chak 2024-3-28
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.

更多回答（2 个）

Sam Chak
Sam Chak 2024-3-25
Ran in:
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更早的评论
Sam Chak
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
Hussam 2024-3-28
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.

请先登录，再进行评论。

Rykelle
Rykelle 2024-10-24
Hello! I am also building a Windkessel right now and found this source helpful. It has code included :>

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by