Optimal control of SEIR with RK4 method problem on updating

Hello everyone,

I am trying to implement an optimal control problem using Runge-Kutta 4th order for a SEIR model with two different categories. My code is running and provides an optimal control but the state variables S,E,I and R remain as if no intervention occurs, which means that the update of the second part is somehow not implemented in it? I don't understand where is the problem. I ran it some times and the results were fine but then all of a sudden it just gives me S,E,I and R as if no control is imposed on the model. Can you please have a look? code follows

function y=odeNEW

test=-1;

T=400;

deltaError=0.001;

M=1000;

t=linspace(0,T,M+1);

h=T/M;

h2=h/2;

C=0.001; K=1000; B=1;

g=0.0625; bcc=0.25; bca=0.11; bac=0.11; baa=0.34;

Sc=zeros(1,M+1);

Sa=zeros(1,M+1);

Ic=zeros(1,M+1);

Ia=zeros(1,M+1);

Rc=zeros(1,M+1);

Ra=zeros(1,M+1);

%init

Sc(1)=0.199; Sa(1)=0.699; Ic(1)=0.001; Ia(1)=0.001; Ra(1)=0; Rc(1)=0;

u=zeros(1,M+1);

LSc=zeros(1,M+1); LSa=zeros(1,M+1); LIc=zeros(1,M+1); LIa=zeros(1,M+1);

%final time

LSc(M+1)=0; LSa(M+1)=0; LIc(M+1)=0; LIa(M+1)=0;

J=zeros(1,M+1);

while (test<0)

oldu=u;

oldSc=Sc;

oldSa=Sa;

oldIc=Ic;

oldIa=Ia;

oldRc=Rc;

oldRa=Ra;

oldLSa=LSa;

oldLSc=LSc;

oldLIc=LIc;

oldLIa=LIa;

for i=1:M

m11=-u(i)*bcc*Sc(i)*Ic(i)-bca*u(i)*Sc(i)*Ia(i);

m12=bcc*u(i)*Sc(i)*Ic(i)+bca*u(i)*Sc(i)*Ia(i)-g*Ic(i);

m13=g*Ic(i);

m14=-bac*u(i)*Sa(i)*Ic(i)-baa*u(i)*Sa(i)*Ia(i);

m15=bac*u(i)*Sa(i)*Ic(i)+baa*Sa(i)*u(i)*Ia(i)-g*Ia(i);

m16=g*Ia(i);

%

aux=0.5*(u(i)+u(i+1));

m21=-aux*bcc*(Sc(i)+h2*m11)*(Ic(i)+h2*m12)-bca*aux*(Sc(i)+h2*m11)*(Ia(i)+h2*m15);

m22=aux*bcc*(Sc(i)+h2*m11)*(Ic(i)+h2*m12)+bca*aux*(Sc(i)+h2*m11)*(Ia(i)+h2*m15)-g*(Ic(i)+h2*m12) ;

m23=g*(Ic(i)+h2*m12) ;

m24=-aux*bac*(Sa(i)+h2*m14)*(Ic(i)+h2*m12)-baa*aux*(Sa(i)+h2*m14)*(Ia(i)+h2*m15);

m25=bac*aux*(Sa(i)+h2*m14)*(Ic(i)+h2*m12) +baa*aux*(Sa(i)+h2*m14)*(Ia(i)+h2*m15)-g*(Ia(i)+h2*m15) ;

m26=g*(Ia(i)+h2*m15) ;

%

m31=-aux*bcc*(Sc(i)+h2*m21)*(Ic(i)+h2*m22)-bca*aux*(Sc(i)+h2*m21)*(Ia(i)+h2*m25);

m32=aux*bcc*(Sc(i)+h2*m21)*(Ic(i)+h2*m22)+bca*aux*(Sc(i)+h2*m21)*(Ia(i)+h2*m25)-g*(Ic(i)+h2*m22) ;

m33=g*(Ic(i)+h2*m22) ;

m34=-aux*bac*(Sa(i)+h2*m24)*(Ic(i)+h2*m22)-baa*aux*(Sa(i)+h2*m24)*(Ia(i)+h2*m25);

m35=bac*aux*(Sa(i)+h2*m24)*(Ic(i)+h2*m22) +baa*aux*(Sa(i)+h2*m24)*(Ia(i)+h2*m25)-g*(Ia(i)+h2*m25);

m36=g*(Ia(i)+h2*m25);

%

aux=u(i+1);

m41=-aux*bcc*(Sc(i)+h*m31)*(Ic(i)+h*m32)-bca*aux*(Sc(i)+h*m31)*(Ia(i)+h*m35);

m42=aux*bcc*(Sc(i)+h*m31)*(Ic(i)+h*m32)+bca*aux*(Sc(i)+h*m31)*(Ia(i)+h*m35)-g*(Ic(i)+h*m32);

m43=g*(Ic(i)+h*m32);

m44=-aux*bac*(Sa(i)+h*m34)*(Ic(i)+h*m32)-baa*aux*(Sa(i)+h*m34)*(Ia(i)+h*m35);

m45=bac*aux*(Sa(i)+h*m34)*(Ic(i)+h*m32) +baa*aux*(Sa(i)+h*m34)*(Ia(i)+h*m35)-g*(Ia(i)+h*m35) ;

m46=g*(Ia(i)+h*m35) ;

%

Sc(i+1)=Sc(i)+(h/6)*(m11 + 2*m21 + 2*m31 + m41);

Ic(i+1)=Ic(i)+(h/6)*(m12 + 2*m22 + 2*m32 + m42);

Rc(i+1)=Rc(i)+(h/6)*(m13 + 2*m23 + 2*m33 + m43);

Sa(i+1)=Sa(i)+(h/6)*(m14 + 2*m24 + 2*m34 + m44);

Ia(i+1)=Ia(i)+(h/6)*(m15 + 2*m25 + 2*m35 + m45);

Ra(i+1)=Ra(i)+(h/6)*(m16 + 2*m26 + 2*m36 + m46);

end

for i=1:M %backward

j=M+2-i;

n11=LSc(j)*(bcc*u(j)*Ic(j)+bca*u(j)*Ia(j))-LIc(j)*(bcc*u(j)*Ic(j)+bca*u(j)*Ia(j));

auxx=B*K*exp(K*(C-(Ic(j)+(Ia(j)))));

n12=auxx + LSc(j)*bcc*u(j)*Sc(j) - LIc(j)*(bcc*u(j)*Sc(j)+bca*u(j)*Sc(j)-g)+ LSa(j)*bac*Sa(j)*u(j)+ LIa(j)*(-bac*Sa(j)*u(j)) ;

n13=LSa(j)*(bac*u(j)*Ic(j)+baa*u(j)*Ia(j))-LIa(j)*(bac*u(j)*Ic(j)+baa*u(j)*Ia(j)-g);

n14=auxx + LSc(j)*bca*u(j)*Sc(j) -LIc(j)*bca*u(j)*Sc(j)+LSa(j)*baa*u(j)*Sa(j)- LIa(j)*(baa*u(j)*Sa(j)-g);

%

n21=(LSc(j)-h2*n11)*(bcc*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1)))+(bca*0.5*(u(j)+u(j-1))*0.5*(Ia(j)+Ia(j-1)))-(LIc(j)-h2*n12)*bcc*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1))+bca*0.5*(u(j)+u(j-1))*0.5*(Ia(j)+Ia(j-1));

auxx=B*K*exp(K*(C-(Ic(j)+Ic(j-1)+(Ia(j)+Ia(j-1)))));

n22= auxx+(LSc(j)-h2*n11)*bcc*0.5*(u(j)+u(j-1))*0.5*(Sc(j)+Sc(j-1)) - (LIc(j)-h2*n12)*bcc*0.5*(u(j)+u(j-1))*0.5*(Sc(j)+Sc(j-1))+bca*0.5*(u(j)+u(j-1))*(0.5*(Sc(j)+Sc(j-1))-g)+ (LSa(j)-h2*n13)*bac*0.5*(Sa(j)+Sa(j-1))*0.5*(u(j)+u(j-1))+ (LIa(j)-h2*n14)*(-bac*0.5*(Sc(j)+Sc(j-1)));

n23=(LSa(j)-h2*n13)*(bac*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1))+baa*0.5*(u(j)+u(j-1))*0.5*(Ia(j)+Ia(j-1))) -(LIa(j)-h2*n14)*(bac*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1))-g);

n24= auxx+ (LSc(j)-h2*n11)*bca*0.5*(u(j)+u(j-1))*0.5*(Sc(j)-Sc(j-1)) -(LIc(j)-h2*n12)*bca*0.5*(u(j)+u(j-1))*0.5*(Sc(j)-Sc(j-1))+(LSa(j)-h2*n13)*baa*0.5*(u(j)+u(j-1))*0.5*(Sa(j)+Sa(j-1))- (LIa(j)-h2*n14)*baa*0.5*(u(j)+u(j-1))*(0.5*(Sa(j)+Sa(j-1))-g);

%

n31=(LSc(j)-h2*n21)*bcc*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1))+bca*0.5*(u(j)+u(j-1))*0.5*(Ia(j)+Ia(j-1))-(LIc(j)-h2*n22)*bcc*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1))+bca*0.5*(u(j)-u(j-1))*0.5*(Ia(j)+Ia(j-1));

auxx=B*K*exp(K*(C-(0.5*(Ic(j)+Ic(j-1))+0.5*((Ia(j)+Ia(j-1))))));

n32= auxx+(LSc(j)-h2*n21)*bcc*0.5*(u(j)+u(j-1))*0.5*(Sc(j)+Sc(j-1)) - (LIc(j)-h2*n22)*(bcc*0.5*(u(j)+u(j-1))*0.5*(Sc(j)+Sc(j-1))+bca*0.5*(u(j)+u(j-1))*0.5*(Sc(j)+Sc(j-1))-g)+ (LSa(j)-h2*n23)*bac*0.5*(Sa(j)+Sa(j-1))+ (LIa(j)-h2*n24)*(-bac*0.5*(Sc(j)+Sc(j-1)));

n33=(LSa(j)-h2*n23)*(bac*0.5*(u(j)+u(j-1))*(0.5*(Ic(j)+Ic(j-1))+baa*0.5*(u(j)+u(j-1))*0.5*(Ia(j)+Ia(j-1)))) -(LIa(j)-h2*n24)*(bac*0.5*(u(j)+u(j-1))*0.5*(Ic(j)+Ic(j-1))-g);

n34=auxx+ LSc(j)*bca*0.5*(u(j)+u(j-1))*0.5*(Sc(j)+Sc(j-1))-(LIc(j)-h2*n22)*bca*0.5*(u(j)+u(j-1))*0.5*(Sc(j)-Sc(j-1))+(LSa(j)-h2*n23)*baa*0.5*(u(j)+u(j-1))*0.5*(Sa(j)+Sa(j-1))-(LIa(j)-h2*n24)*baa*0.5*(u(j)-u(j-1))*(0.5*(Sa(j)+Sa(j-1))-g);

%

n41=(LSc(j)-h*n31)*(bcc*u(j-1)*Ic(j-1)+bca*u(j-1)*Ia(j-1))-(LIc(j)-h*n32)*bcc*u(j-1)*Ic(j-1)+bca*u(j-1)*Ia(j-1);

auxx=B*K*exp(K*(C-(Ic(j-1)+Ia(j-1))));

n42= auxx+(LSc(j)-h*n31)*bcc*u(j-1)*Sc(j-1)-(LIc(j)-h*n32)*(bcc*u(j-1)*Sc(j-1)+bca*u(j-1)*Sc(j-1))+(LSa(j)-h*n33)*bac*Sa(j-1)+(LIa(j)-h*n34)*(-bac*Sc(j-1));

n43=(LSa(j)-h*n33)*(bac*u(j-1)*Ic(j-1)+baa*u(j-1)*(Ia(j-1)-g));

n44=auxx+ (LSc(j)-h*n31)*bca*u(j-1)*Sc(j-1) -(LIc(j)-h*n32)*bca*u(j-1)*Sc(j-1)+(LSa(j)-h*n33)*baa*u(j-1)*Sa(j-1)- (LIa(j)-h*n34)*(baa*u(j-1)*Sa(j-1)-g);

%

LSc(j-1) = LSc(j)-(h/6)*( n11 + 2*n21 + 2*n31 + n41 ) ;

LIc(j-1) = LIc(j)-(h/6)*( n12 + 2*n22 + 2*n32 + n42 ) ;

LSa(j-1) = LSa(j)-(h/6)*( n13 + 2*n23 + 2*n33 + n43 ) ;

LIa(j-1) = LIa(j)-(h/6)*( n14 + 2*n24 + 2*n34 + n44 ) ;

end

%new control vector

for i=1:M+1

vAux(i)=0.5*(bcc*LSc(i)*Sc(i)*Ic(i)+bca*Sc(i)*Ia(i)-LIc(i)*(bcc*Sc(i)*Ic(i)+bca*Sc(i)*Ia(i))+LSa(i)*(bac*Sa(i)*Ic(i)+baa*Sc(i)*Ia(i))-LIa(i)*(bac*Sa(i)*Ic(i)+baa*Sa(i)*Ia(i)));

auxU = min([max([0 vAux(i)]) 0.9]);

u(i) = 0.5* (auxU + oldu(i));

end

b=10^2;

J= Ic(M+1)+Rc(M+1)+Ia(M+1)+Ra(M+1)-trapz( t,b*(u .^2) );

%absolute error for convergence

temp1=deltaError*sum(abs(Sc))-sum(abs(oldSc-Sc));

temp2=deltaError*sum(abs(Sa))-sum(abs(oldSa-Sa));

temp3=deltaError*sum(abs(Ic)-sum(abs(oldIc-Ic)));

temp4=deltaError*sum(abs(Ia))-sum(abs(oldIa-Ia));

temp5=deltaError*sum(abs(u))-sum(abs(oldu-u));

temp6=deltaError*sum(abs(LSc))-sum(abs(oldLSc-LSc));

temp7=deltaError*sum(abs(LSa))-sum(abs(oldLSa-LSa));

temp8=deltaError*sum(abs(LIc))-sum(abs(oldLIc-LIc));

% temp11=deltaError*sum(abs(LRc))-sum(abs(oldLRc-LRc));

%temp12=deltaError*sum(abs(LRa))-sum(abs(oldLRa-LRa));

temp9=deltaError*sum(abs(Rc))-sum(abs(oldRc-Rc));

temp10=deltaError*sum(abs(Ra))-sum(abs(oldRa-Ra));

test = min(temp1,min(temp2,min(temp3,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10))))))))));

plot(t,u)

end

[Sonu Lamba]

Hi @nota siachouli,

I am also facing the similar problem. Have solved the issue? If so, can you please help me too?

[nota siachouli]

Hi! I have solved my problem. What can I help you with?

[Mahmood Dad]

Dear nota siachouli,

Hope you feel good.

I want to have the model you are reviewing as well as the modified MATLAB code if possible.

I would be very grateful if you could send it to me.

[e-mail]

With best wishes

Mahmood

[mallela ankamma rao]

good evening sir

i am also doing optimal control theory for covid-19 model SEIAQHR

i dont how to draw graphs for infected , susceptible with control and without control like below jpg

if you can give code relating these graphs ,I would be very grateful to you.

if possible please send any reference codes for graphs to mail id [e-mail]

Thanking you