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

조회 수: 15 (최근 30일)
Hussam 2024년 3월 25일
답변: Sam Chak 2024년 3월 28일
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
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.

댓글을 달려면 로그인하십시오.

채택된 답변

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

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Sam Chak 2024년 3월 25일
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
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일
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일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Chemistry에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by