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

Question

Hussam 2024-3-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098256-optimisation-for-3-element-windkessel-model-not-working-using-least-squares-method-and-fminsearch

评论： Hussam 2024-3-25

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:

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

Torsten 2024-3-24

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098256-optimisation-for-3-element-windkessel-model-not-working-using-least-squares-method-and-fminsearch#answer_1430196

移动：Torsten 2024-3-24

在 MATLAB Online 中打开

Ref_P_Q_WK_No-Osc_1-Cycle.csv

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

%dQ = dQ_in(closestIndex);

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

end

8 个评论
显示 6更早的评论隐藏 6更早的评论

Hussam 2024-3-25

编辑：Hussam 2024-3-25

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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

请先登录，再进行评论。

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

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

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

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

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

8 个评论
显示 6更早的评论隐藏 6更早的评论