Avoid ode15s from freezing in parameter optimization

(Version: R2018b) I am running PSO to try to find the parameters that minimizes the distance between my data and simulation. But, I noticed that ode15s keeps getting stuck/freezes and halting the search. I have two equations with four sets of inital values (i.e., 8 input-outputs). I have tried changing the RelTol, AbsTol,and InitialStep without any impact in speed.

TLDR; for a simpler system with same behaviour see next title.

The ODE model is coded as follows:

function dydt = system(~, y, pars)
    alpha = pars(1:2);
    delta = pars(3:4);
    power = pars(5:10);
    y = reshape(y, 2, 4)';
    dydt = zeros(4,2);
    dydt(:,1) = alpha(1).*y(:,1).^power(1) - delta(1).*y(:,1).^power(2).*y(:,2).^power(3);
    dydt(:,2) = alpha(2).*y(:,1).^power(4).*y(:,2).^power(5) - delta(2).*y(:,2).^power(6);
    dydt = dydt';
    dydt = dydt(:);
end

As sugested I include the Jacobian for the model:

function j = jacobian(~, y, pars)
    alpha = pars(1:2);
    delta = pars(3:4);
    power = pars(5:10);
    y = reshape(y, 2, 4);
    j = zeros(2, 8);
    j(1,:) = [alpha(1).*power(1).*y(1,:).^(power(1)-1)-delta(1).*power(2).*y(1,:).^(power(2)-1).*y(2,:).^power(3)...
        -delta(1).*power(3).*y(1,:).^power(2).*y(2,:).^(power(3)-1)...
        ];
    j(2,:) = [alpha(2).*power(4).*y(1,:).^(power(4)-1).*y(2,:).^power(5)...
         alpha(2).*power(5).*y(1,:).^power(4).*y(2,:).^(power(5)-1)-delta(2).*power(6).*y(2,:).^(power(6)-1)...
         ];
     j = reshape(j, 8, 2);
     j = blkdiag(j(1:2,:), j(3:4,:), j(5:6,:), j(7:8,:));
end

And a very slow example I found is the following:

y0 = [0.9477    0.6914    0.8712    1.3908    0.9905    1.1709    1.1806    0.8595];
pars = [1.6608    0.9739    4.7934    2.8388    3.4574    0.5525    3.1752    0.1559    4.7070    1.1896];
tspan = linspace(0, 720, 721);
opts = odeset(...
'Jacobian', @(t,y) jacobian(t, y, pars),...
'NonNegative', ones(1, numel(y0)),...
'Vectorized', 'on'...
);
[~, sim] = ode15s(@(t,y) system(t, y, pars), tspan, y0, opts);

TLDR; starts here.

The next thing I did was to only use one set of initial values to make the problem easier.

function dydt = simple_system(~, y, pars)
    alpha = pars(1:2);
    delta = pars(3:4);
    power = pars(5:10);
    dydt = zeros(2,1);
    dydt(1) = alpha(1).*y(1).^power(1) - delta(1).*y(1).^power(2).*y(2).^power(3);
    dydt(2) = alpha(2).*y(1).^power(4).*y(2).^power(5) - delta(2).*y(2).^power(6);
end

and likewise made the Jacobian the same way:

function j = simple_jacobian(~, y, pars)
    alpha = pars(1:2);
    delta = pars(3:4);
    power = pars(5:10);
    j = [alpha(1).*power(1).*y(1).^(power(1)-1)-delta(1).*power(2).*y(1).^(power(2)-1).*y(2).^power(3)...
        -delta(1).*power(3).*y(1).^power(2).*y(2).^(power(3)-1);...
        alpha(2).*power(4).*y(1).^(power(4)-1).*y(2).^power(5)...
         alpha(2).*power(5).*y(1).^power(4).*y(2).^(power(5)-1)-delta(2).*power(6).*y(2).^(power(6)-1)...
         ];
end

When I run the simple model, it no longer freezes on the parameters and the inital values used before and correctly throws a warning.

So I checked that my Jacobian function returned the correct value, and from what I could tell it did. That is, it returned a block diagonal matrix with 4x4 groups of values for each set of intial conditons with values equal to running the simple_jacobian. . Still, I tried to run simulations in sequential order instead of vectorized over intial values.

But, the solver still gets stuck (not if I output the plot though), and if anyone could help me understand why that is, it would be tremoundesly helpful.

Here is an example for the simple system:

y02 = [1.4648    0.8389];
pars2 = 0.8389    0.7218    0.9540    0.6352    0.7432    0.7835    0.8443    0.2612    0.9067    0.5347;
opts2 = odeset(...
'Jacobian', @(t,y) simple_jacobian(t, y, pars2),...
'NonNegative', ones(1, numel(y02)),...
'Vectorized', 'on'...
);
% Throws an error
ode15s(@(t,y) simple_system(t, y, pars2), tspan, y02, opts2);
% Just freezes
[~,sim] = ode15s(@(t,y) simple_system(t, y, pars2), tspan, y02, opts2);

Answer 1

[Philip Berg]

The best solution to fix these types of problems when 'NonNegative' assumption is made is to add an event that cancels the integration.

Add the following event:

function [position,isterminal,direction] = simple_EventsFcn(t,y)
    %y1 event
    position(1) = y(1) - .0001; % The value that we want to be zero; can be similar to e.g., error tolerance
    isterminal(1) = 1;  % Halt integration 
    direction(1) = 0; % The zero can be approached from either direction
    %y2 event
    position(2) = y(2) - .0001; % The value that we want to be zero
    isterminal(2) = 1;  % Halt integration 
    direction(2) = 0; % The zero can be approached from either direction
end

The integration will now halt early and will not run very slow.

Redefining the ode options and runnig the same parameters:

y02 = [1.4648    0.8389];
pars2 = [0.8389    0.7218    0.9540    0.6352    0.7432    0.7835    0.8443    0.2612    0.9067    0.5347];
tspan = [0, 15.00001];
opts = odeset(...
'Jacobian', @(t,y) simple_jacobian(t, y, pars2),...
'Events', @simple_EventsFcn,...
'NonNegative', ones(1, numel(y02)),...
'Vectorized', 'on'...
);
tic;
[~, sim] = ode15s(@(t,y) simple_system(t, y, pars2), tspan, y02, opts);
toc;
%Elapsed time is 0.017697 seconds.

If one is doing optimisation like I am, checking if the number of rows of sim is the same as the length of tspan determines if the parameters could be integrated over the entire time-course or not. Else one can return NaN to the optimisation algorithm.

Answer 2

[Jan]

I've limite the time to [0, 15]. You see that one component explodes between t=15 and t=16. This let the step size of the integration shrink to tiny values. The integration does not freeze, but it runs extremely slow only.

This is a mathematical limitation, not a problem of Matlab.

y02 = [1.4648    0.8389];
pars2 = [0.8389    0.7218    0.9540    0.6352    0.7432    0.7835    0.8443    0.2612    0.9067    0.5347];
tspan = [0, 15.00001];
opts = odeset(...
'Jacobian', @(t,y) simple_jacobian(t, y, pars2),...
'NonNegative', ones(1, numel(y02)),...
'Vectorized', 'on'...
);

a = ode15s(@(t,y) simple_system(t, y, pars2), tspan, y02, opts);
plot(a.x, a.y)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

function dydt = simple_system(~, y, pars)
    alpha = pars(1:2);
    delta = pars(3:4);
    power = pars(5:10);
    dydt = zeros(2,1);
    dydt(1) = alpha(1).*y(1).^power(1) - delta(1).*y(1).^power(2).*y(2).^power(3);
    dydt(2) = alpha(2).*y(1).^power(4).*y(2).^power(5) - delta(2).*y(2).^power(6);
end

function j = simple_jacobian(~, y, pars)
    alpha = pars(1:2);
    delta = pars(3:4);
    p     = pars(5:10);
    j = [alpha(1).*p(1).*y(1).^(p(1)-1)-delta(1).*p(2).*y(1).^(p(2)-1).*y(2).^p(3)...
        -delta(1).*p(3).*y(1).^p(2).*y(2).^(p(3)-1);...
         alpha(2).*p(4).*y(1).^(p(4)-1).*y(2).^p(5)...
         alpha(2).*p(5).*y(1).^p(4).*y(2).^(p(5)-1)-delta(2).*p(6).*y(2).^(p(6)-1)...
         ];
end

[Philip Berg]

Yeah I see, but is there a way to have MATLAB throw a warning at these points instead (e.g., by some odeset options)? Why does it sometimes give a warning that the error tol is too small and sometimes just runs very slow? Would it help to reformulate the probelm as described in Fit ODE, Problem-Based?