Script stuck running forever

34 views (last 30 days)
Kevin Burke
Kevin Burke on 3 May 2022
Answered: Steven Lord on 3 May 2022
I have been trying to develop a script that can solve for unknown parameters in an SIRD model. I followed Star Strider's parameter estimation setup and it has lead to an unending loading loop when I press run. Can anybiody identify an issue here?
Thanks in Advance,
Kevin
P0=rand(3,1)*100; %Parameter initial guesses
US_Population = 330000000; %US population
USData = xlsread('USData.xlsx'); %Importing data
IRD=USData([40:70],[2:4]); %Isolating IRD March 2020 data
time=[1:31]'; %Creating a 31 day column vector
Infected = IRD(:,1); %Isolates cases column vector
S_data = US_Population - Infected; %S(t): Population - I(t)
US_March_2020=[S_data IRD];
[P,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Objective,P0,time,S_data);
function S = Objective(P,t)
%ODEs
%dS/dt = -beta*S*I
%dI/dt = beta*S*I - lambda*I -kd*I
%dR/dt = Lambda*I
%dD/dt = kd*I
%Variables: x(1)=S, x(2)=I, x(3)=R, x(4)=D
%Parameters: beta = P(1), lambda = P(2), kd = P(3)
%Adding this duplicate line of code so that X0 is definied in the function
US_Population = 330000000; %US population
USData = xlsread('USData.xlsx'); %Importing data
IRD=USData([40:70],[2:4]); %Isolating IRD March 2020 data
Infected = IRD(:,1); %Isolates cases column vector
S_data = US_Population - Infected; %S(t): Population - I(t)
US_March_2020=[S_data IRD];
x0= US_March_2020(1,:);
[T,Sv] = ode45(@SIRD, t, x0);
function dS = SIRD(t,x)
xdot = zeros(4,1);
xdot(1) = (-P(1).*x(1).*x(2));
xdot(2) = (P(1).*x(1).*x(2)-P(2).*x(2)-P(3).*x(2));
xdot(3) = (P(2).*x(2));
xdot(4) = (P(3).*x(2));
dS=xdot;
end
S = Sv(:,1);
end

Answers (2)

Walter Roberson
Walter Roberson on 3 May 2022
I doubt that code is creating an unending loop.
I think it much more likely that the ode is "stiff" for at least a range of inputs, and that it is simply taking a long long long time to calculate the answer.
Try a stiff solver such as ode15s.
  3 Comments
Walter Roberson
Walter Roberson on 3 May 2022
Yes, just replace ode45 by ode15s
Kevin Burke
Kevin Burke on 3 May 2022
This script runs now, except the only problem is that the output parameters are identical to the initial guesses that I put in. I have changed them several times to no avail. Do you know what is causing this and how I can fix it? I get the attached message as well. (included at bottom)
"Initial point is a local minimum. Optimization complete because the size of the gradient at the initial point is less than the value of the optimality tolerance."

Sign in to comment.


Steven Lord
Steven Lord on 3 May 2022
One obvious optimization relates to the section of code starting with this comment:
%Adding this duplicate line of code so that X0 is definied in the function
That code calls xlsread once per call that lsqcurvefit makes to your objective function. Use one of the techniques from this documentation page to pass the data that you read using xlsread into your objective function as an additional parameter so you don't keep accessing the disk to read the file.

Community Treasure Hunt

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

Start Hunting!