Parameter Estimation for SIR model using Particle Swarm Optimization.

Hullo Matlab community. I am stack with Parameter Estimation of SIR (Susceptible-Infected-Recovered) model using Particle Swarm Optimization(PSO) method. I need to estimate two parameters, beta and delta. I am trying to modify a program to suit my needs. The program has 3 file i.e. ofun (objective function), sir (ode model) and pso (script file with PSO code and data). I don’t know how to modify the objective function(ofun) so that the code can run past line 42 of pso script and estimate the parameters. Below are the original and my modified files for parameter estimation. Please help, I will be grateful. Thank you.

link for original code: https://www.researchgate.net/post/curve_fitting_using_PSOhttps://www.researchgate.net/post/curve_fitting_using_PSO

%Sir model

function dx = sir(t, x)

global beta delta;

dx = [0; 0; 0];

dx(1) = -beta * x(1) * x(2);

dx(2) = beta * x(1) * x(2) - delta * x(2);

dx(3) = delta * x(2);

end

%objective function ofun

function f=ofun(x,data)

% for i = 1:length(data(:,1))

% [t, y] = ode15s('sir', x(i), [ 500 2 0 ]);

% end

[t, y] = ode15s('sir', [1 77], [ 500 2 0 ]);

yvalue = interp1(t, y(:,2), data(:,1));

SSE = sum((yvalue - data(:,2)).^2);

f=SSE;

%pso script

tic

clc

clear all

%parameter estimation of SIR model using PSO

%% Input data

tdata= [1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;...

26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;...

49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;...

72;73;74;75;76;77];

ydata=[19;40;38;34;39;43;58;89;89;106;133;129;119;139;400;441;415;392;450;600;...

832;812;513;825;700;460;493;852;998;1117;1131;832;1043;800;1087;784;...

780;838;760;515;374;553;382;458;382;347;245;231;496;434;267;320;318;200;...

138;270;202;168;155;191;104;236;202;184;121;121;99;104;76;87;66;127;144;131;23;103;45];

data = [tdata;ydata]; %*****%

%*****%

%%

LB=[0 0 0 ]; %lower bounds of variables %*****%

UB=[1 10 10]; %upper bounds of variables %*****%

% pso parameters values

m=3; % number of variables(parameters???) %*****%

n=100; % population size

wmax=0.9; % inertia weight

wmin=0.4; % inertia weight

c1=2; % acceleration factor

c2=2; % acceleration factor

% pso main program----------------------------------------------------start

maxite=500; % set maximum number of iteration %*****%

maxrun=10; % set maximum number of runs need to be

for run=1:maxrun

run

% pso initialization----------------------------------------------start

for i=1:n

for j=1:m

x0(i,j)=round(LB(j)+rand()*(UB(j)-LB(j)));

end

end

x=x0; % initial population

v=0.1*x0; % initial velocity

for i=1:n

f0(i,1)=ofun(x0(i,:),data) %*****%

end

[fmin0,index0]=min(f0);

pbest=x0; % initial pbest

gbest=x0(index0,:); % initial gbest

% pso initialization------------------------------------------------end

% pso algorithm---------------------------------------------------start

ite=1;

tolerance=1;

while ite<=maxite && tolerance>10^-12

w=wmax-(wmax-wmin)*ite/maxite; % update inertial weight

% pso velocity updates

for i=1:n

for j=1:m

v(i,j)=w*v(i,j)+c1*rand()*(pbest(i,j)-x(i,j))...

+c2*rand()*(gbest(1,j)-x(i,j));

end

end

% pso position update

for i=1:n

for j=1:m

x(i,j)=x(i,j)+v(i,j);

end

end

% handling boundary violations

for i=1:n

for j=1:m

if x(i,j)<LB(j)

x(i,j)=LB(j);

elseif x(i,j)>UB(j)

x(i,j)=UB(j);

end

end

end

% evaluating fitness

for i=1:n

f(i,1)=ofun(x(i,:),data); %*****%

end

% updating pbest and fitness

for i=1:n

if f(i,1)<f0(i,1)

pbest(i,:)=x(i,:);

f0(i,1)=f(i,1);

end

end

[fmin,index]=min(f0); % finding out the best particle

ffmin(ite,run)=fmin; % storing best fitness

ffite(run)=ite; % storing iteration count

% updating gbest and best fitness

if fmin<fmin0

gbest=pbest(index,:);

fmin0=fmin;

end

% calculating tolerance

if ite>100;

tolerance=abs(ffmin(ite-100,run)-fmin0);

end

% displaying iterative results

if ite==1

disp(sprintf('Iteration Best particle Objective fun'));

end

disp(sprintf('%8g %8g %8.4f',ite,index,fmin0));

ite=ite+1;

end

% pso algorithm-----------------------------------------------------end

gbest;

fvalue=ofun(gbest,data); %*****%

%fvalue=10*(gbest(1)-1)^2+20*(gbest(2)-2)^2+30*(gbest(3)-3)^2;

fff(run)=fvalue;

rgbest(run,:)=gbest;

disp(sprintf('--------------------------------------'));

end

% pso main program------------------------------------------------------end

disp(sprintf('\n'));

disp(sprintf('*********************************************************'));

disp(sprintf('Final Results-----------------------------'));

[bestfun,bestrun]=min(fff)

best_variables=rgbest(bestrun,:)

disp(sprintf('*********************************************************'));

toc

% PSO convergence characteristic

plot(ffmin(1:ffite(bestrun),bestrun),'-k');

xlabel('Iteration');

ylabel('Fitness function value');

title('PSO convergence characteristic')

[crt56 ak]

Used star stride's example