pid for dynamic malaria

Can the pid controller be set to malaria dynamics?

[Star Strider]

That would appear to require a state space realisation.

However, it would likely be easier to simply use one of the ODE solvers (most likely ode15s) to integrate it.

[af af]

How do I do this?

function = malaria(dShdt,dEhdt;dIhdt;dRhdt;dSvdt;dEvdt;dIvdt)

phi=0.502;epsilon=0.2;beta=0.8333;landa=0.09;muh=0.00004;muv=0.1429;
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;
Nh=Sh+Eh+Ih+Rh;   %total human population
Nv=Sv+Ev+Iv;      %total mosquito population
betam=beta*epsilon*phi*Iv/Nh;
landav=landa*epsilon*phi*Ih/Nh;

dShdt=lambdah+(k*Rh)-(1-u1(t))*betam*Sh-(muh*Sh);
dEhdt=(1-u1(t))*betam*Sh-(a1+muh)*Eh;
dIhdt=(a1*Eh)-(b+tau*u2(t))*Ih-(psi+muh)*Ih;
dRhdt=(b+tau*u2(t))*Ih-(k+muh)*Rh;
dSvdt=lambdav-(1-u1(t))*landav*Sv-(p*u3(t)*Sv)-(muv*Sv);
dEvdt=(1-u1(t))*landav*Sv-p*u3(t)*Ev-(a2+muv)*Ev;
dIvdt=(a2*Ev)-(p*u3(t)*Iv)-(muv*Iv);
end

[Star Strider]

It will be necessary to pass the arguments to your function as a vector, not specific separate arguments. See the documentation for ode45, ode15s, and the others to understand how to code your ODE function correctly. (Also include the time variable as the first argument in the argument list, even if you do not use it in your code.)

[af af]

can you fixed? (It has error)

tspan=[0 100];
y0=zeros(1,7);
[T,Y]=ode45(malariaSEIRS,tspan,y0);
figure(1)
plot(T,y)
grid
function dydt=malariaSEIRS(t,y)
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;

dydt(1)=lambdah+(k*Rh)-(1-u1(t))*betam*Sh-muh*Sh;
dydt(2)=(1-u1(t))*betam*Sh-(a1+muh)*Eh;
dydt(3)=a1*Eh-(b+tau*u2(t))*Ih-(psi+muh)*Ih;
dydt(4)=(b+tau*u2(t))*Ih-(k+muh)*Rh;
dydt(5)=lambdav-(1-u1(t))*landav*Sv-p*u3(t)*Sv-muv*Sv;
dydt(6)=(1-u1(t))*landav*Sv-p*u3(t)*Ev-(a2+muv)*Ev;
dydt(7)=a2*Ev-p*u3(t)*Iv-muv*Iv;

dydt=[dydt(1),dydt(2),dydt(3),dydt(4),dydt(5),dydt(6),dydt(7)];
end

[Walter Roberson]

[T,Y]=ode45(@malariaSEIRS,tspan,y0);

[af af]

i put [T,Y]=ode45(@malariaSEIRS,tspan,y0);

but it still gives an error

[Walter Roberson]

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.

Answer 1

[Sam Chak]

Hi @af af

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