Avoid ode15s from freezing in parameter optimization

1 次查看（过去 30 天）
Philip Berg
Philip Berg 2021-9-14
回答： Philip Berg 2021-9-14
(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);

请先登录，再回答此问题。

采纳的回答

Philip Berg
Philip Berg 2021-9-14
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.

更多回答（1 个）

Jan
Jan 2021-9-14
编辑：Jan 2021-9-14
Ran in:
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
  1 个评论
Philip Berg
Philip Berg 2021-9-14
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by