MATLAB Answers

Avoid ode15s from freezing in parameter optimization

조회 수: 4(최근 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일
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?

댓글을 달려면 로그인하십시오.

태그

Community Treasure Hunt

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

Start Hunting!

Translated by