Optimal control of SEIR with RK4 method problem on updating

18 次查看（过去 30 天）
nota siachouli
nota siachouli 2020-8-12
评论： Shivam 2024-9-28
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
  5 个评论
显示 3更早的评论
mallela ankamma rao
mallela ankamma rao 2022-7-26
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 ushanand.mallela@gmail.com
Thanking you
Shivam
Shivam 2024-9-28
Hey @nota siachouli can you please provide the corrected code(if you got). Actually I am working on the similar type of problem. If you provide the scheme to solve optimal control problem for SEIR model, it will be a great help. My email-id for your reference is d21021@students.iitmandi.ac.in

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Biological and Health Sciences 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by