pid for dynamic malaria

1 visualización (últimos 30 días)
af af
af af el 21 de En. de 2021
Respondida: Sam Chak el 20 de Ag. de 2022
Can the pid controller be set to malaria dynamics?
  7 comentarios
af af
af af el 24 de En. de 2021
i put [T,Y]=ode45(@malariaSEIRS,tspan,y0);
but it still gives an error
Walter Roberson
Walter Roberson el 24 de En. de 2021
You appear to have functions named u1 and u2 and u3 that each take t as a parameter, but you have not posted the code for those functions.

Iniciar sesión para comentar.

Respuestas (1)

Sam Chak
Sam Chak el 20 de Ag. de 2022
Not sure what are PID controllers for the Malaria dynamics. However, by the mathematical manipulation, you can probably propose as shown below to keep the disease spread under control.
tspan = [0 10];
y0 = zeros(1, 7);
[T, Y] = ode45(@malariaSEIRS, tspan, y0);
plot(T, Y), grid, xlabel('Days')
legend('S_h', 'E_h', 'I_h', 'R_h', 'S_v', 'E_v', 'I_v', 'location', 'East')
Y(end, :)
ans = 1×7
0.0001 1.5163 0.4067 0.0021 0.9990 125.2982 0.8853
function dydt = malariaSEIRS(t, y)
dydt = zeros(7, 1);
% parameters
phi = 0.502;
epsilon = 0.2;
beta = 0.8333;
landa = 0.09;
muh = 0.00004;
muv = 0.1429;
k = 0.7902;
a1 = 1/17;
a2 = 1/18;
lambdah = 0.2;
lambdav = 1000;
tau = 0.01 - 0.7;
psi = 0.05;
b = 0.005;
p = 0.25;
Sh = 1100;
Eh = 200;
Ih = 400;
Rh = 0;
Sv = 800;
Ev = 250;
Iv = 80;
Nh = Sh + Eh + Ih + Rh;
% Nv = Sv + Ev + Iv;
betam = beta*epsilon*phi*Iv/Nh;
landav = landa*epsilon*phi*Ih/Nh;
% states
Sh = y(1);
Eh = y(2);
Ih = y(3);
Rh = y(4);
Sv = y(5);
Ev = y(6);
Iv = y(7);
% u1, u2, u3
k3 = 1; % adjust this parameter
u3 = ((k3 - muv)*Iv + a2*Ev)/p;
u2 = 0;
k1 = 1; % adjust this parameter
u1 = (- (k1 - p*u3 - muv - 1*landav)*Sv - lambdav)/landav;
% Malaria dynamics
dydt(1) = lambdah + (k*Rh) - (1 - u1)*betam*Sh - muh*Sh;
dydt(2) = (1 - u1)*betam*Sh - (a1 + muh)*Eh;
dydt(3) = a1*Eh - (b + tau*u2)*Ih - (psi + muh)*Ih;
dydt(4) = (b + tau*u2)*Ih - (k + muh)*Rh;
dydt(5) = lambdav - (1 - u1)*landav*Sv - p*u3*Sv - muv*Sv;
dydt(6) = (1 - u1)*landav*Sv - p*u3*Ev - (a2 + muv)*Ev;
dydt(7) = a2*Ev - p*u3*Iv - muv*Iv;
end

Community Treasure Hunt

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

Start Hunting!

Translated by