Fittg covid 19 SEIRD model

2 次查看（过去 30 天）
Hothayfa
Hothayfa 2025-3-3
回答： Star Strider 2025-3-3
fitted code does not align the real data
can you check it for me ?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start code
%% Step 1: Load Real Data from Excel
clc; clear; close all;
% Load data
filepath = 'C:\Users\Dell\Desktop\SEIRD_output.xlsx';
data = readtable(filepath);
% Extract relevant columns (first 256 days)
t_real = data.days(1:256); % Time vector
Real_S = data.Susceptibles(1:256);
Real_E = data.Exposed(1:256);
Real_I = data.Infected(1:256);
Real_R = data.Recovered(1:256);
Real_D = data.Died(1:256);
tspan = linspace(1, max(t_real), 256);
N = data.population(1);% Define Population Size
% Initial Conditions (from dataset)
y0 = [Real_S(1), 10, 5, Real_R(1), Real_D(1)]; % [S0, E0, I0, R0, D0]
%% Step 2: Define SEIRD Model ODE System
mu = 0.016; % Mortality rate due to COVID (fixed)
f = 3.5e-6; % Natural mortality rate (fixed)
b = 3e-6; % Birth rate (fixed)
seird_model = @(t, y, params) [
-(params(1) * y(1) * y(3) / N + f * y(1) + params(6) * y(1)) + (params(5) * y(4) + b * N); % dS/dt
-(f * y(2) + params(2) * y(2)) + (params(1) * y(1) * y(3) / N + params(4) * y(4)); % dE/dt
-(f * y(3) + params(3) * y(3)) + params(2) * y(2); % dI/dt
-(f * y(4) + params(5) * y(4) + params(4) * y(4)) + params(3) * (1 - mu) * y(3) + params(6) * y(1); % dR/dt
params(3) * mu * y(3)]; % dD/dt
%% Step 3: Define Grid Search Ranges for Parameters
beta_range = linspace(0.01, 1, 10);
sigma_range = linspace(0.001, 1, 10);
alpha_range = linspace(0.001, 1, 10);
eta_range = linspace(0.001, 1, 5);
xi_range = linspace(0.0001, 1, 5);
nu_range = linspace(0.00001, 0.001, 5);
best_params = [];
best_error = Inf;
%% Step 4: Grid Search for Rough Parameter Estimation
for beta = beta_range
for sigma = sigma_range
for alpha = alpha_range
for eta = eta_range
for xi = xi_range
for nu = nu_range
params = [beta, sigma, alpha, eta, xi, nu];
tspan = linspace(0, max(t_real), 500); % Increased resolution
[t, y] = ode45(@(t, y) seird_model(t, y, params), tspan, y0);
% Interpolate to match real data times
S = interp1(t, y(:,1), t_real, 'pchip');
E = interp1(t, y(:,2), t_real, 'pchip');
I = interp1(t, y(:,3), t_real, 'pchip');
R = interp1(t, y(:,4), t_real, 'pchip');
D = interp1(t, y(:,5), t_real, 'pchip');
% Compute Weighted Error
error = abs(Real_I - I) ./ max(Real_I) * 10 + ...
abs(Real_S - S) ./ max(Real_S) + ...
abs(Real_E - E) ./ max(Real_E) * 5 + ...
abs(Real_R - R) ./ max(Real_R) + ...
abs(Real_D - D) ./ max(Real_D);
% ✅ Fix: Corrected the if condition (minimization)
if sum(error) < best_error
best_error = sum(error);
best_params = params;
end
end
end
end
end
end
end
% ✅ Fix: Ensure best_params is always set
if isempty(best_params)
best_params = [0.3, 1/14, 0.1, 0.01, 0.001, 0.0001]; % Default values
end
%% Step 5: Fine-Tune with Nonlinear Optimization (lsqcurvefit)
objective_function = @(params, ~) simulate_seird_error(params, t_real, y0, Real_S, Real_E, Real_I, Real_R, Real_D, N);
t_dummy = ones(size(t_real)); % Dummy time input for lsqcurvefit
lb = [0.01, 0.001, 0.0001, 0.001, 0.001, 0.000001]; % Lower bounds
ub = [10, 1, 10, 10, 10, 0.001]; % Upper bounds
optimal_params = lsqcurvefit(objective_function, best_params, t_dummy, zeros(size(t_dummy)), lb, ub);
%% Step 6: Run the SEIRD Model with Optimized Parameters and Plot
[t, y] = ode45(@(t, y) seird_model(t, y, optimal_params), tspan, y0);
S = interp1(t, y(:,1), t_real, 'pchip');
E = interp1(t, y(:,2), t_real, 'pchip');
I = interp1(t, y(:,3), t_real, 'pchip');
R = interp1(t, y(:,4), t_real, 'pchip');
D = interp1(t, y(:,5), t_real, 'pchip');
figure;
subplot(5,1,1);
plot(t_real, Real_S, 'g', t_real, S, 'b--', 'LineWidth', 1.5);
title('Susceptibles (S)'); legend('Real', 'Fitted'); grid on;
subplot(5,1,2);
plot(t_real, Real_E, 'g', t_real, E, 'b--', 'LineWidth', 1.5);
title('Exposed (E)'); legend('Real', 'Fitted'); grid on;
subplot(5,1,3);
plot(t_real, Real_I, 'g', t_real, I, 'b--', 'LineWidth', 1.5);
title('Infected (I)'); legend('Real', 'Fitted'); grid on;
subplot(5,1,4);
plot(t_real, Real_R, 'g', t_real, R, 'b--', 'LineWidth', 1.5);
title('Recovered (R)'); legend('Real', 'Fitted'); grid on;
subplot(5,1,5);
plot(t_real, Real_D, 'g', t_real, D, 'b--', 'LineWidth', 1.5);
title('Deceased (D)'); legend('Real', 'Fitted'); grid on;
sgtitle('SEIRD Model Fitting to Real Data');
%% Supporting Function for Optimization
function error_vec = simulate_seird_error(params, t_real, y0, Real_S, Real_E, Real_I, Real_R, Real_D, N)
seird_model = @(t, y) [
-(params(1) * y(1) * y(3) / N + f * y(1) + params(6) * y(1)) + (params(5) * y(4) + b * N); % dS/dt
-(f * y(2) + params(2) * y(2)) + (params(1) * y(1) * y(3) / N + params(4) * y(4)); % dE/dt
-(f * y(3) + params(3) * y(3)) + params(2) * y(2); % dI/dt
-(f * y(4) + params(5) * y(4) + params(4) * y(4)) + params(3) * (1 - mu) * y(3) + params(6) * y(1); % dR/dt
params(3) * mu * y(3)]; % dD/dt
[t, y] = ode45(@(t, y) seird_model(t, y), t_real, y0);
S = interp1(t, y(:,1), t_real, 'pchip');
E = interp1(t, y(:,2), t_real, 'pchip');
I = interp1(t, y(:,3), t_real, 'pchip');
R = interp1(t, y(:,4), t_real, 'pchip');
D = interp1(t, y(:,5), t_real, 'pchip');
error_vec = log(1 + abs(Real_I - I) * 10) + log(1 + abs(Real_S - S)) + ...
log(1 + abs(Real_E - E) * 10) + log(1 + abs(Real_R - R)) + log(1 + abs(Real_D - D));
end

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2025-3-3
See: Parameter Estimation for a System of Differential Equations among other examples on how to do this.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by