(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)
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);
As sugested I include the Jacobian for the model:
function j = jacobian(~, y, pars)
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 = blkdiag(j(1:2,:), j(3:4,:), j(5:6,:), j(7:8,:));
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);
'Jacobian', @(t,y) jacobian(t, y, pars),...
'NonNegative', ones(1, numel(y0)),...
[~, 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)
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);
and likewise made the Jacobian the same way:
function j = simple_jacobian(~, y, pars)
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)...
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:
pars2 = 0.8389 0.7218 0.9540 0.6352 0.7432 0.7835 0.8443 0.2612 0.9067 0.5347;
'Jacobian', @(t,y) simple_jacobian(t, y, pars2),...
'NonNegative', ones(1, numel(y02)),...
ode15s(@(t,y) simple_system(t, y, pars2), tspan, y02, opts2);
[~,sim] = ode15s(@(t,y) simple_system(t, y, pars2), tspan, y02, opts2);