It took a few minutes to determine what the problems are with this (there were several).
I got it to work. Most of the problems were with respect to passing all the necessary arguments to the various funcitons. I will let you explore the changes I made, probably the most important of which was using ‘@SITRModel’ to pass the function as an argument.
Try this —
% Define the SITR model as a system of ODEs
function dydt = SITRModel(t, y, beta, gamma, delta, alpha, lambda, mu, N)
S = y(1); % Susceptible
I = y(2); % Infectious
Q = y(3); % Isolated
T = y(4); % Treated
R = y(5); % Recovered
dSdt = -beta * S * I / N;
dIdt = beta * S * I / N - gamma * I - delta * I - alpha * I;
dQdt = delta * I - lambda * Q;
dTdt = alpha * I - mu * T;
dRdt = gamma * I + lambda * Q + mu * T;
dydt = [dSdt; dIdt; dQdt; dTdt; dRdt];
end
% Observed data (replace with actual data)
% Format: [time, infected, isolated, treated, recovered]
observed_data = [0, 1, 0, 0, 0;
10, 50, 10, 5, 15;
20, 100, 25, 15, 50;
30, 150, 35, 30, 100;
40, 200, 50, 50, 200];
% Initial conditions
N = 1000000; % Total population
S0 = N - 1;
I0 = 1;
Q0 = 0;
T0 = 0;
R0 = 0;
y0 = [S0, I0, Q0, T0, R0];
% Time points for the solution (based on observed data)
tspan = observed_data(:, 1);
% Define the objective function for optimization
function error = objectiveFunction(params,tspan,SITRModel,y0,N,observed_data)
beta = params(1);
gamma = params(2);
delta = params(3);
alpha = params(4);
lambda = params(5);
mu = params(6);
% Solve the ODE with the current parameters
[t, y] = ode45(@(t, y) SITRModel(t, y, beta, gamma, delta, alpha, lambda, mu, N), tspan, y0);
% Interpolate the model's output to match the time points of observed data
model_values = interp1(t, y(:, 2:5), tspan);
% Calculate the sum of squared errors
error = sum((model_values - observed_data(:, 2:5)).^2, 'all');
end
% % Initial guess for parameters [beta, gamma, delta, alpha, lambda, mu]
initial_params = [0.3, 0.1, 0.05, 0.02, 0.1, 0.1];
initial_params = initial_params + rand(1,6)/100
% % Perform optimization using fminsearch
options = optimset('MaxFunEvals', 1000, 'MaxIter', 1000);
optimized_params = fminsearch(@(params)objectiveFunction(params,tspan,@SITRModel,y0,N,observed_data), initial_params, options);
% Display optimized parameters
disp('Optimized parameters:');
disp(optimized_params);
% optimized_params = initial_params;
% Solve the system with optimized parameters
[t, y] = ode45(@(t, y) SITRModel(t, y, optimized_params(1), optimized_params(2), optimized_params(3), optimized_params(4), optimized_params(5), optimized_params(6), N), linspace(0, 160, 100), y0);
% Plot the solution with optimized parameters
figure;
plot(t, y(:, 1), 'b-', t, y(:, 2), 'r-', t, y(:, 3), 'g-', t, y(:, 4), 'm-', t, y(:, 5), 'k-');
legend('Susceptible', 'Infectious', 'Isolated', 'Treated', 'Recovered');
title('Fitted SITR Model for COVID-19');
xlabel('Days');
ylabel('Population');
grid on;
.
