COVID-19 SEIR model

I tries to estimate the parameters of SEIR model using ode15s and fmincon. But it shows the following issue:

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than

the value of the step size tolerance and constraints are

satisfied to within the value of the constraint tolerance.

the main code and functions are as following

Main Code

close all; clear all; clc

global t Confirmed Recovered Deaths Active N

%insert reported data

[Date, Days, Confirm, Recover, Death] = getDataCOVID_BD("C:\academic\MIST\COVID-19\SEIR-master\SEIR-master\Corona .xlsx", "Req_data_BD_new", [24,54]);

Confirmed = Confirm;

Recovered = Recover;

Deaths = Death;

Active = Confirmed-Recovered-Deaths; %quarantined cases

time = Date;

N = numel(time); %% number of time points

dt = 1;

tn = [0:N-1].*dt;

t = tn';

% Tdt =1/24;

% Ttime = datetime(time(1)):Tdt:datetime(Date(end));

% TN = numel(Ttime);

% Ttn = [0:TN-1].*dt;

% Tt = Ttn';

% Population = N

N = 1.61e8;

%% Parameter initial guess (WHO's prescribed minimum values are taken)

beta = 1;

alpha= 0.1;

gamma= 0.5;

delta= 0.07;

lamda= [0.1,0.05];

kappa= [0.1,0.05];

p0 = [beta,alpha,gamma,delta,lamda,kappa];

%% show initial objective

disp(['Initial SSE Objective: ' num2str(obj_1(p0))])

% optimize parameters

% no linear constraints

A = [];

b = [];

Aeq = [];

beq = [];

nlcon = [];

% optimize with fmincon

lb = [0.5,0.01,0.2,0.07,10,2,5,9]; % lower bound

ub = [1.2,0.12,1,0.5,45,15,9,10]; % upper bound

% options = optimoptions(@fmincon,'Algorithm','interior-point');

p = fmincon(@obj_1,p0,A,b,Aeq,beq,lb,ub,nlcon); %,options);

% show final objective

disp(['Final SSE Objective: ' num2str(obj_1(p))])

% optimized parameter values

beta = p(1);

alpha = p(2);

gamma = p(3);

delta =p(4);

lamda = [p(5),p(6)];

kappa = [p(7),p(8)];

disp(['beta: ' num2str(beta)])

disp(['alpha: ' num2str(alpha)])

disp(['gamma: ' num2str(gamma)])

disp(['delta: ' num2str(delta)])

disp(['lamda: ' num2str(lamda)])

disp(['kappa: ' num2str(kappa)])

% calculate model with updated parameters

Tp = sim_1(p);

%%

%%plot result

figure

semilogy(t,Tp(:,4),'r',t,Tp(:,5),'b',t,Tp(:,6),'k');

hold on

semilogy(t,Active,'ro',t,Recovered,'bo',t,Deaths,'ko');

% ylim([0,1.1*Npop])

ylabel('Number of cases')

xlabel('time (days)')

leg = {'Active(fitted)',...

'Recovered (fitted)','Death (fitted)',...

'Active(reported)',...

'Recovered (reported)','Death (reported)'};

legend(leg{:},'location','southoutside')

set(gcf,'color','w')

grid on

axis tight

set(gca,'yscale','lin')

sim_1 function

function T = sim_1(p)

global t Confirmed Recovered Deaths Active N

T = zeros(length(t),7);

T (1,1)= N -(3*Confirmed(1)); %S0

T (1,2)= Confirmed(1); %E0

T (1,3)= Confirmed(1); %I0

T (1,4)= Confirmed(1)-Recovered(1)-Deaths(1);%Q0

T (1,5)= Recovered(1); %R0

T (1,6)= Deaths (1); %D0

T (1,7)= 0; %P0

T0 = [ T(1,1),T(1,2),T(1,3),T(1,4),T(1,5),T(1,6),T(1,7)];

for i = 1:length(t)-1

ts = [t(i),t(i+1)];

sol = ode15s(@(t,x)SEIRQPD(t,x,N,p),ts,T0);

T0 = [sol.y(1,end),sol.y(2,end),sol.y(3,end),sol.y(4,end)...

sol.y(5,end),sol.y(6,end),sol.y(7,end)];

T(i+1,1) = T0(1);

T(i+1,2) = T0(2);

T(i+1,3) = T0(3);

T(i+1,4) = T0(4);

T(i+1,5) = T0(5);

T(i+1,6) = T0(6);

T(i+1,7) = T0(7);

end

end

obj_1 function

% define objective

function obj = obj_1(p)

global Confirmed Recovered Deaths Active

% simulate model

Tp = sim_1(p);

% calculate objective

obj = sum(((Tp(:,4)-Active)./Active).^2) ...

+ sum(((Tp(:,5)-Recovered)./Recovered).^2)...

+ sum(((Tp(:,6)-Deaths)./Deaths).^2);

end

SEIRQPD function

% define SEIRQDP model

function dTdt = SEIRQPD(t,x,N,p)

beta =p(1);

alpha=p(2);

gamma=p(3);

delta=p(4);

lamda=p(5);

kappa=p(7);

% SEIRQPD States

S = x(1);

E = x(2);

I = x(3);

Q = x(4);

R = x(5);

D = x(6);

P = x(7);

% SEIRQPD ODE

dSdt = -(beta*S*I/N)-(alpha*S);

dEdT = (beta*S*I/N)- (gamma*E);

dIdT = (gamma*E)-(delta*I);

dQdT = (delta*I)-(lamda*Q)-(kappa*Q);

dRdT = lamda*Q;

dDdT = kappa*Q;

dPdT = alpha*S;

dTdt = [dSdt,dEdT,dIdT,dQdT,dRdT,dDdT,dPdT]';

end

[Star Strider]

There is no problem.

The fmincon function found a solution, although it warns that it could be a local minimum. Use the estimated ‘p’ parameters in the function to see if it fits your data.

[Samia Sarothi]

Actually the model doesn't fit with data also. i think there is some issue with sim_1.