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)
endHello everyone,

function y=odeNEW
test=-1;
T=400;
deltaError=0.001;
M=1000;
t=linspace(0,T,M+1);
h=T/M;
h2=h/2;

for i=1:M %backward
j=M+2-i;

function y=odeNEW
test=-1;
T=400;
deltaError=0.001;
M=1000;
t=linspace(0,T,M+1);
h=T/M;
h2=h/2;

for i=1:M %backward
j=M+2-i;

Cart

Cart