Fminsearch does not work after increasing maxiter

Hi everyone

I am new to Matlab. I come from Python which is slow for the type of problems that I am working on so I learned Matlab recently.

I am running an fminsearch to find the optimal paramters to fit a model in physiology. When I set the Maxiter to 10, the code ran all right but the result was not good. Naturally, I tried to increase the Maxiter number. The code did not run anymore and I got the following warning

Warning: Failure at t=8.127018e-01. Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (2.887297e-15) at time t.

I have attached below the code that I used. I can provide the data as well if it is neccesary.

Thank you very much.

Tung Nguyen

----

function BSF

clear all 
close all 
clc 

org_data = load("apnea_heart1.mat")
all_data = org_data.data;
ecg = all_data(1:12340);
[peaks, locs] = findpeaks(ecg, 'MinPeakHeight', 1.3*10^-3);
no_data_points = locs(2)-locs(1)+1;
time_frequency = 10^-3;
end_time = no_data_points*time_frequency;
time_span = linspace(0, no_data_points*time_frequency, no_data_points);
pressure = all_data(24681+locs(1):24681+locs(2));
flow = all_data(12341+locs(1):12341+locs(2));
tforward = time_span; 


initial_guess = [0.2  0.004 0];
r = initial_guess(1);
l = initial_guess(2);
ic = initial_guess(3);



                               
   function dydt = model(t,y, para_fit)
    P = interp1(time_span, pressure, t);
    R = para_fit(1);
    L = para_fit(2);
    dydt = [-R/L*y(1)+1/L * P];
    end

    function error_in_data = moder(k)
    % computing the error in the data
    
    
    % solves the ODE, output is written in t and y 
    IC = k(3);
    [t,y] = ode23s(@(t,y)model(t,y,k), tforward, IC);
    time_observed = 2000:500:6000; 
    error_in_data = sum((y(time_observed)-flow(time_observed)').^2);
    end 

k = [r l ic];
% k = [paras(1) paras(2)];
IC = k(3);
[t,y] = ode23s(@(t,y)model(t,y,k), tforward, IC);
figure(1)
subplot(1,2,1);
plot(tforward, flow, 'b');
hold on
plot(t, y, 'r'); % plot the solution with the initial data
xlabel('time');
ylabel('Blood flow');
legend("Measured", "Initial guess");
axis([0 end_time -40 200]);

options = optimset('Maxiter', 20);
[opt_paras, fval] = fminsearch(@moder, k, options);
disp(opt_paras);
IC = opt_paras(3);

% para_fit = [2.000792312628345, 0.04934547048449 0];
% k = [para_fit(1) para_fit(2)];
% IC = para_fit(3);

[t,y] = ode23s(@(t,y)model(t,y,opt_paras), tforward, IC);
figure(1)
subplot(1,2,2);
plot(tforward, flow, 'b');
hold on
plot(t, y, 'r'); % plot the solution with the initial data
xlabel('time');
ylabel('Blood flow');
axis([0 end_time -40 200]);
legend("Measured", "Optimized guess");
end 

[Star Strider]

‘Warning: Failure at t=8.127018e-01. Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (2.887297e-15) at time t.’

This is not a problem with fminsearch. Shortly after the integration begins (at about 0.8 time units) at some point in your code, the ‘model’ function encounters a singularity (±Inf), and stops. The most likely explanation is that ‘L’ becomes zero (or close to it).

[Tung]

Thank you very much for your help. That makes sense!

Answer 1

[Star Strider]

Following up on my Comment, see if this solves the problem —

function dydt = model(t,y, para_fit)

L = max(1,L);

P = interp1(time_span, pressure, t);

R = para_fit(1);

L = para_fit(2);

dydt = [-R/L*y(1)+1/L * P];

end

That essentially constrains ‘L’ to be no less than 1. If it can be less than 1 (and greater than 0), replace that value for 1 in the max call.

[Tung]

Thank you a lot for your answer. Can you please explain why the equation essential contrains 'L' to be no less than 1?

I have a follow up question. Instead of requiring L>1 (which does not happen in my problem because L is often small), can I set A = -R/L, B = 1/L, solve the optimization for A,B and then derive R and L from A and B? Is it a common practice?

[Star Strider]

As always, my pleasure!

It constrains ‘L’ because I added this assignment:

L = max(1,L);

so whenever ‘L’ is less than 1, the max function returns whatever value is larger of either ‘L’ or 1. (Ideally, you would use Lagrange multipliers to constrain ‘L’, so this is a relatively priimitive approach, although it should work. The idea is to prevent ‘L’ from becoming 0 and forcing ‘dydt’ to become infiinite, stopping the integration.)

To change the default value for ‘L’ to perhaps , that assignment becomes:

L = max(1E-4,L);

I would keep ‘L’ and ‘R’ separate, as they currently are. The reason is that they might very well diverge, and keeping them separate allows them to be estimated separately and individually. Combining parameters by addition or multiplication is not generally considered to be good programming practice, specifically because the confidence limits of the parameters need to be calculated separately so they can be assessed separetely. Their statistical characteristics are important.

You can certainly run your code combining them as ‘A’ and ‘B’, however I believe your results will be better if they are kept as they currently are, since combining them is not common practice and is generally discouraged. You can certainly try that experiment (run them separately and combined) and then compare the results. I believe you will find that keeping them separate is preferable.

.

[Tung]

Thank you! That makes perfect sense.

[Star Strider]

As always, my pleasure!

Answer 2

[Matt J]

The warning is not coming from fminsearch. It is coming from the ODE solver. As fminsearch iteratively explores different choices for k, it sometimes lands on one for which the ODE becomes numerically pathological. Only you know what this means for the problem at hand, and how it should be handled. Perhaps you should abort moder() with error_in_data=inf when this occurs, so that fminsearch will recognize such k as a bad parameter choice.

In any case, I don't recommend artificially truncating MaxIter. You should probably abort the ODE solver when it starts to choke.

[Tung]

Thank you for your answer. That makes sense to me. To deal with the error, I imagine that we should use syntax similar to try/except in Python?