Parameter estimation of SEIHRDBP model

2 Mar 2024
1 Answer
Answer Accepted
Updated 6 Mar 2024
7 Views (30 days)

0 votes

% I want simulate "parameter estimation of SEIHRDBP Model", where
% S=susceptible, E=exposed, I=infected, H=hospitalised, R=recovered,
% D=death, B=buried, P=pathogen. I tried Matlab codes from the Matlab
% community but couldn't get right. The parameters here are assumed values
% Parameter value
% Delta = params(1); mu = params(2); delta = params(3); gamma = params(4); gamma1 = params(5);
% eta = params(6); xi = params(7); d = params(8); b = params(9); phip = params(10);
% d1 = params(11); b1 = params(12); alpha = params(13); sigma = params(14);
% eta1 = params(15); beta = params(16); beta1 = params(17); beta2 = params(18); beta3 = params(19);
% betap = params(20);
% params= [10, 0.5, 0.5, 0.005, 0.12, 0.12, 0.6, 0.8, 0.5, 0, 0.04, 0.04, 0.03, 0.04, 0.04, 0.006, 0.002, 0.012, 0.0012; 0.01];
% dxdt = zeros(8, 1); % initializing the state variables
% params = zeros(20, 1); % initializing the parameters
%
% dxdt(1) = params(1) - (params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1)-params(2)*x(1);
% dxdt(2) =(params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1) - (params(2) + params(3))*x(2);
% dxdt(3) = params(3)*x(2) - (params(2) + params(4) + params(5))*x(1);
% dxdt(4) = params(5)*x(3) - (params(6) + params(7) + params(2))*x(4);
% dxdt(5) = params(7)*x(4) - params(2)*x(5);
% dxdt(6) = params(4)*x(3) + params(6)*x(4) -params(8)*x(6);
% dxdt(7) = params(8)*x(6) - params(9)*x(7);
% dxdt(8) = params(10) + params(14)*x(3) + params(15)*x(4) + params(11)*x(6) + params(12)*x(7) -params(13)*x(8);
% This is code I tried. I basically used this help site: https://uk.mathworks.com/matlabcentral/answers/1710260-need-help-with-ode-parameter-estimation
% The c1---c8 are the state variables, while the theta1---theta20 are the parameter values.
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218 0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192 0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289 0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233 0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756 0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255 0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526 0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401 0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277 0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323 0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702 0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835 0.002595 0.08301 0.08224 0.835];
theta0 = rand(20,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@SEIHRDBP,theta0,t,c,zeros(size(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = SEIHRDBP(theta, tv);
% figure
% hd = plot(t, c, 'p');
% for k1 = 1:size(c,2)
% CV(k1,:) = hd(k1).Color;
% hd(k1).MarkerFaceColor = CV(k1,:);
% end
% hold on
% hlp = plot(tv, Cfit);
% for k1 = 1:size(c,2)
% hlp(k1).Color = CV(k1,:);
% end
% hold off
% grid
% xlabel('Time')
% ylabel('Population')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
%
function C=SEIHRDBP(theta,t)
c0=theta(end-3:end);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(8,1);
dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1)-theta(2)*c(1);
dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
dcdt(8) = theta(10) + theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
dC=dcdt;
end
C=Cv;
end
% Please can some help check it for me. I have attached the observed caese of all the state variables(SEIHRDBP)

Sign in to answer this question.

 Accepted Answer

Star Strider
Star Strider on 2 Mar 2024
Edited: Star Strider on 2 Mar 2024
Ran in:

0 votes

Ran in:
Thank you for quoting my code!
The problem was here:
c0=theta(end-3:end);
use this instead:
c0=theta(end-7:end);
Since my code uses the last elements of the paramter vector for the initial conditions of the differential equations (it fits them as parameters, as well as the rate constants), that needs to be accounted for. Here there are 8 differential equations, so there must be 8 initial conditions.
Try this —
% I want simulate "parameter estimation of SEIHRDBP Model", where
% S=susceptible, E=exposed, I=infected, H=hospitalised, R=recovered,
% D=death, B=buried, P=pathogen. I tried Matlab codes from the Matlab
% community but couldn't get right. The parameters here are assumed values
% Parameter value
% Delta = params(1); mu = params(2); delta = params(3); gamma = params(4); gamma1 = params(5);
% eta = params(6); xi = params(7); d = params(8); b = params(9); phip = params(10);
% d1 = params(11); b1 = params(12); alpha = params(13); sigma = params(14);
% eta1 = params(15); beta = params(16); beta1 = params(17); beta2 = params(18); beta3 = params(19);
% betap = params(20);
% params= [10, 0.5, 0.5, 0.005, 0.12, 0.12, 0.6, 0.8, 0.5, 0, 0.04, 0.04, 0.03, 0.04, 0.04, 0.006, 0.002, 0.012, 0.0012; 0.01];
% dxdt = zeros(8, 1); % initializing the state variables
% params = zeros(20, 1); % initializing the parameters
%
% dxdt(1) = params(1) - (params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1)-params(2)*x(1);
% dxdt(2) =(params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) + params(20)*x(8))*x(1) - (params(2) + params(3))*x(2);
% dxdt(3) = params(3)*x(2) - (params(2) + params(4) + params(5))*x(1);
% dxdt(4) = params(5)*x(3) - (params(6) + params(7) + params(2))*x(4);
% dxdt(5) = params(7)*x(4) - params(2)*x(5);
% dxdt(6) = params(4)*x(3) + params(6)*x(4) -params(8)*x(6);
% dxdt(7) = params(8)*x(6) - params(9)*x(7);
% dxdt(8) = params(10) + params(14)*x(3) + params(15)*x(4) + params(11)*x(6) + params(12)*x(7) -params(13)*x(8);
T1 = readtable('observed.csv')
T1 = 35×8 table
susceptible exposed infected Hospitalized Recovered deaths Buried Pathogens ___________ _______ ________ ____________ _________ ______ ______ _________ 3444 323 164 164 0 55 55 0 3444 426 1 1 0 1 1 0 3444 300 5 5 0 5 5 0 3444 335 11 11 0 9 9 0 3444 386 12 12 0 6 6 0 3444 350 23 23 0 12 12 0 3444 404 130 130 0 55 55 0 3444 416 3470 3470 0 2287 2287 0 3444 449 54 54 0 33 33 0 3444 446 8 8 0 4 4 0 3444 374 66 66 0 49 49 0 34440 38646 28646 28646 0 11323 11323 0 3444 377 6 6 0 3 3 0 3444 491 36 36 0 13 13 0 3444 491 11 11 0 4 4 0 3444 437 1 1 0 1 1 0
VN = T1.Properties.VariableNames;
L = size(T1,1);
t = linspace(0, L-1, L).'; % Time Vector (NECESSARY) Units = ?
% % This is code I tried. I basically used this help site: https://uk.mathworks.com/matlabcentral/answers/1710260-need-help-with-ode-parameter-estimation
% % The c1---c8 are the state variables, while the theta1---theta20 are the parameter values.
% t=[0.1
% 0.2
% 0.4
% 0.6
% 0.8
% 1
% 1.5
% 2
% 3
% 4
% 5
% 6];
% c=[0.902 0.06997 0.02463 0.00218 0.902 0.06997 0.02463 0.00218
% 0.8072 0.1353 0.0482 0.008192 0.8072 0.1353 0.0482 0.008192
% 0.6757 0.2123 0.0864 0.0289 0.6757 0.2123 0.0864 0.0289
% 0.5569 0.2789 0.1063 0.06233 0.5569 0.2789 0.1063 0.06233
% 0.4297 0.3292 0.1476 0.09756 0.4297 0.3292 0.1476 0.09756
% 0.3774 0.3457 0.1485 0.1255 0.3774 0.3457 0.1485 0.1255
% 0.2149 0.3486 0.1821 0.2526 0.2149 0.3486 0.1821 0.2526
% 0.141 0.3254 0.194 0.3401 0.141 0.3254 0.194 0.3401
% 0.04921 0.2445 0.1742 0.5277 0.04921 0.2445 0.1742 0.5277
% 0.0178 0.1728 0.1732 0.6323 0.0178 0.1728 0.1732 0.6323
% 0.006431 0.1091 0.1137 0.7702 0.006431 0.1091 0.1137 0.7702
% 0.002595 0.08301 0.08224 0.835 0.002595 0.08301 0.08224 0.835];
c = table2array(T1);
theta0 = rand(28,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@SEIHRDBP,theta0,t,c,zeros(size(theta0)));
Local minimum possible. lsqcurvefit stopped because the size of the current step is less than the value of the step size tolerance.
fprintf(1,'\tRate Constants:\n')
Rate Constants:
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
Theta(1) = 2.40444 Theta(2) = 0.03125 Theta(3) = 1.41359 Theta(4) = 1.11936 Theta(5) = 1.11908 Theta(6) = 0.30996 Theta(7) = 0.06624 Theta(8) = 0.60836 Theta(9) = 0.41146 Theta(10) = 0.29405 Theta(11) = 0.29017 Theta(12) = 0.13046 Theta(13) = 1.72754 Theta(14) = 0.69723 Theta(15) = 0.38384 Theta(16) = 1.29961 Theta(17) = 0.99986 Theta(18) = 0.78319 Theta(19) = 0.54491 Theta(20) = 0.44378 Theta(21) = 0.41763 Theta(22) = 1.15486 Theta(23) = 1.41118 Theta(24) = 0.89700 Theta(25) = 0.14331 Theta(26) = 0.62479 Theta(27) = 0.24113 Theta(28) = 0.66698
tv = linspace(min(t), max(t));
Cfit = SEIHRDBP(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Population')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
legend(hlp, VN, 'Location','best')
%
function C=SEIHRDBP(theta,t)
c0=theta(end-7:end);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(8,1);
dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1)-theta(2)*c(1);
dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
dcdt(8) = theta(10) + theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
dC=dcdt;
end
C=Cv;
end
% Please can some help check it for me. I have attached the observed caese of all the state variables(SEIHRDBP)
That works, however it is not an escellent fit. I would keep funning it with random initial parameter vectyors until you get a good fit. Using the ga (genetic algorithm) function might produce better results.
.

10 Comments
Show 8 older comments Hide 8 older comments

University
University on 2 Mar 2024
Edited: University on 2 Mar 2024
Thank you so very much for you help. Please how does the ga (genetic algorithm) works? And why is it better?
I tried to run it my laptop, it was sucessful but subsequently, it returns this error:
Warning: Failure at t=5.304862e+00. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (1.421085e-14) at time t.
> In ode45 (line 350)
In sample>SEIHRDBP (line 58)
In lsqcurvefit (line 256)
In sample (line 29)
Error using lsqcurvefit
Function value and YDATA sizes are not equal.
Error using lsqcurvefit
Function value and YDATA sizes are not equal.
Error in sample (line 56)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@SEIHRDBP,theta0,t,c,zeros(size(theta0)));
Error using helpview
optim not found.
Star Strider
Star Strider on 2 Mar 2024
Edited: Star Strider on 2 Mar 2024
Here is the genetic algorithm approach to this problem. It takes longer than the allotted 55 seconds to run, so I cannot run it here. However if you have the Global Optimization Toolbox, it should give a better result than the gradient descent approach. Also, change the solver to ode15s, since the system is likely ‘stiff’.
Also, if you run it on your own computer on in MATLAB Online, use the ‘opts’ options structure instead of ‘optsAns’, and change that in the ga call as well (the last argument).
The error results in the ODE integrator encountering a singularity (±Inf) that then stops it because it can go no further. Start the code again and keep running the code, since it will eventually choose a set of parameters that will allow it to run completely and produce a useful result.
Try this —
T1 = readtable('observed.csv')
VN = T1.Properties.VariableNames;
L = size(T1,1);
t = linspace(0, L-1, L).'; % Time Vector (NECESSARY) Units = ?
c = table2array(T1);
ftns = @(theta) norm(c-SEIHRDBP(theta,t));
PopSz = 500;
Parms = 28;
optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-2, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems
% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
% optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1); % With 'HybridFcn'
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
legend(hlp, VN, 'Location','best')
function C=SEIHRDBP(theta,t)
c0=theta(end-7:end);
[T,Cv]=ode15s(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(8,1);
dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1)-theta(2)*c(1);
dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) + theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
dcdt(8) = theta(10) + theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
dC=dcdt;
end
C=Cv;
end
.
University
University on 2 Mar 2024
Thank you so very much for your help. Why do you choose theta0 = (1, 28)? Is it because we have 8 state variables and 20 parameters. Since there 20 parameters, why don't theta0 = (1, 20)?
Star Strider
Star Strider on 2 Mar 2024
That is because I have my code set up to estimate the initial conditions of the differential equations as separate parameters, rather than simply setting them as the initial concentrations in the data (although it is perfectly acceptable to use the first data entries as the initial estimates of the initial conditions).
The parameter vector ‘theta’ is set up so that the first 20 elements (in this instance) are the kinetic parameters, and the last 8 (in this instance) are the initial conditions. Since everything is measured and therefore subject to the same sort of measurement uncertainty, it actually makes sense to estimate the initial conditions as parameters.
The genetic algorithm will probably return the best estimates of the parameters, however it may take it a while. In situations such as yours, when you are estimating 28 parameters, I set ga to execute in a loop for several iterations (usually between 10 and 25, although occasionally as many as 100), saving the parameters and final fitness value (and execution times and the number of generations and other statistics that interest me) in each iteration in a cell array or numerical matrix, and then go off and do something else for a few minutes.
University
University on 2 Mar 2024
Oh, I see. Thank you
Star Strider
Star Strider on 2 Mar 2024
Edited: Star Strider on 3 Mar 2024
As always, my pleasure!
EDIT — (3 Mar 2024 at 17:18)
The data you are currently using are not appropriate for estimating the parameters of your system, because the data does not demonstrate any change. These data are essentially steady-state. You need a different data set.
.
University
University on 4 Mar 2024
Edited: University on 4 Mar 2024
Ran in:
Ran in:
Hi Strider,
I am back please. I have gotten the data correct data with date. The whole idea is that I want the parameter estimation for parameters for SEI. I.e., run the parameter estimation these classes (SEI). These are only the active human class. The parameters associated with classes are theta16, theta17, theta18, theta19, theta20, theta3, theta7 and theta8.
I tried to run the new code use the "ga" algorithm but I got the error:
Arrays have incompatible sizes for this operation.
Error in sample>@(theta)norm(c-SEIHRDBP(theta,t)) (line 9)
ftns = @(theta) norm(c-SEIHRDBP(theta,t));
Error in createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)
fcn_handle = @(x) fcn(x,FcnArgs{:});
In fact, this is my new corrected codes and new data is attached.
SEIHRDBP_MODEL_PARAMETER
Arrays have incompatible sizes for this operation.

Error in solution>SEIHRDBP_MODEL_PARAMETER (line 59)
opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,12)*1E-3 ones(1,7)*10],'MaxGenerations',2E3)%,'PlotFcn','gaplotbestf');
function SEIHRDBP_MODEL_PARAMETER
%7 ODEs to describe the model (S,E,I,H,R,D,B, P)
function M = SEIHRDBP(p,t)
%Initial conditions with population
M0 = p(21:28);
[~,m] = ode15s(@DifEq,t,M0);
function dM = DifEq(~,k)
dMdt = zeros(8,1);
dMdt(1) = p(1) - (p(16)*k(3) + p(17)*k(4) + p(18)*k(6)+ p(19)*k(7) + p(20)*k(8))*k(1)-p(2)*k(1);
dMdt(2) = (p(16)*k(3) + p(17)*k(4) + p(18)*k(6)+ p(19)*k(7) + p(20)*k(8))*k(1) - (p(2) + p(3))*k(2);
dMdt(3) = p(3)*k(2) - (p(2) + p(4) + p(5))*k(1);
dMdt(4) = p(5)*k(3) - (p(6) + p(7) + p(2))*k(4);
dMdt(5) = p(7)*k(4) - p(2)*k(5);
dMdt(6) = p(4)*k(3) + p(6)*k(4) -p(8)*k(6);
dMdt(7) = p(8)*k(6) - p(9)*k(7);
dMdt(8) = p(10) + p(14)*k(3) + p(15)*k(4) + p(11)*k(6) + p(12)*k(7) -p(13)*k(8);
dM = dMdt;
end
M = m;
end
%From day Aug to nov 20
t = [1
2
3
4
5
6
7
8
9
10];
%Official data for I, D for 10 days
k = [24 23
62 35
68 41
70 42
70 42
70 43
68 49
67 49
67 49
66 49];
% [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
function fv = ftns(p,t)
SEIHRDBPout = SEIHRDBP(p,t);
fv = norm(k-SEIHRDBPout(:,[3 6]));
end
% ftns = @(p) norm(k-EIDAGG2(p,t));
PopSz = 50;
% Parameter = 12;
Parameter = 28;
opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,12)*1E-3 ones(1,7)*10],'MaxGenerations',2E3)%,'PlotFcn','gaplotbestf');
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n',t0)
[p,fval,exitflag,output] = ga(@(p)ftns(p,t),Parameter,[],[],[],[],zeros(Parameter,1),[ones(20,1); Inf(8,1)],[],[],opts);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n',t1)
GA_Time = etime(t1,t0);
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf(1,'\n\tRate Constants:\n')
for k1 = 1:12
fprintf(1,'\t\tp(%2d) = %11.8f\n',k1,p(k1))
end
fprintf(['\n\tInitial Conditions:\n'])
for k1 = 1:8
fprintf('\t\tdMdt(%d) = %11.8f\n',k1,p(k1))
end
tv = linspace(min(t),max(t));
Mfit = SEIHRDBP(p,tv);
figure(1)
plot(t,k,'p')
hold on
hlp = plot(tv,Mfit);
hold off
grid
xlabel('time')
ylabel('Population')
legend(hlp,'S','E','I','H','R','D','B','P',Location','best')
end
Star Strider
Star Strider on 5 Mar 2024
I posted a reply to this in your other question Parameter estimation of SEIHRDBP model a few minutes ago.
University
University on 6 Mar 2024
Thank you so much. I was about to reply you. I have gotten it.
Star Strider
Star Strider on 6 Mar 2024
As always, my pleasure!

Sign in to comment.

More Answers (0)

Categories

Find more on Vibration Analysis in Help Center and File Exchange

Products

Release

R2023b

Tags

Community Treasure Hunt

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

Start Hunting!