Parameter estimation of SEIHRDBP model

Question

University on 2 Mar 2024

0
Link

Direct link to this question

https://nl.mathworks.com/matlabcentral/answers/2089451-parameter-estimation-of-seihrdbp-model

Commented: Star Strider on 6 Mar 2024

observed.csv

% 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)  

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

Answer 1

Star Strider on 2 Mar 2024

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/2089451-parameter-estimation-of-seihrdbp-model#answer_1420031

Edited: Star Strider on 2 Mar 2024

Open in MATLAB Online

observed.csv

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 commentsHide 8 older comments

Star Strider on 2 Mar 2024

Edited: Star Strider on 2 Mar 2024

Open in MATLAB Online

observed.csv

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

.

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 on 4 Mar 2024

Edited: University on 4 Mar 2024

Open in MATLAB Online

mydata.csv

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

University on 6 Mar 2024

Thank you so much. I was about to reply you. I have gotten it.

Star Strider on 6 Mar 2024

As always, my pleasure!

Sign in to comment.

Parameter estimation of SEIHRDBP model

0 Comments
Show -2 older commentsHide -2 older comments

Accepted Answer

10 Comments
Show 8 older commentsHide 8 older comments

More Answers (0)

See Also

Categories

Tags

Products

Release

Community Treasure Hunt

Parameter estimation of SEIHRDBP model

0 Comments Show -2 older commentsHide -2 older comments

Accepted Answer

10 Comments Show 8 older commentsHide 8 older comments

More Answers (0)

See Also

Categories

Tags

Products

Release

Community Treasure Hunt

0 Comments
Show -2 older commentsHide -2 older comments

10 Comments
Show 8 older commentsHide 8 older comments