how to combine to get single graph in optimal control problems

13 views (last 30 days)
mallela ankamma rao
mallela ankamma rao on 21 Nov 2022
Answered: Pavl M. on 21 Nov 2024 at 8:45
good afternoon to one and all
Sir i have done research on infectious diseases
in my paper i am using optimal technique with three control variables u1,u2 and u3
i knew how to get graphs separately for with controls, u1 and u2 but i dont know to combine three graphs in one graph
now i request you please help me anyone how to plot single graph
code1
code for all controls (u1,u2 and u3)
clear all; close all;
% SETUP FOR FORWARD-BACKWARD SWEEP METHOD
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 400; % Final time
delta = 0.0001; % Accepted tollerance
M = 999; % Number of nodes
t = linspace(t0,tf,M+1); % Time variable where linspace creates M+1 ...equally spaced nodes between t0 and tf, including t0 and tf.
h = tf/M; % Spacing between nodes
h2 = h/2; % Spacing equal to 2 for Runge-Kutta method
% MODEL PARAMETERS
N = 1390000000;
pi = 70 ;
zetaa = 0.5;
zetas = 0.3;
zetaq = 0.3;
omega = 0.2;
theta = 0.5;
mu = 0.00391;
mua = 0.001945;
muh = 0.001945;
beta= 0.5888;
lamdaa= 0.3679;
lamdas= 0.2393;
etas= 0.1049;
etaq= 0.1917;
gammaa=0.2524;
gammaq= 0.1050;
gammah=0.0930;
% WEIGHT FACTORS
c1 = 1;
c2 = 2;
c3 = 2;
c4 = 2;
c5 = 30;
c6 = 30;
c7 = 30;
% INITIAL CONDITIONS MODEL
x1=zeros(1,M+1);
x2=zeros(1,M+1);
x3=zeros(1,M+1);
x4=zeros(1,M+1);
x5=zeros(1,M+1);
x6=zeros(1,M+1);
x7=zeros(1,M+1);
x1(1) = 13900000;
x2(1) = 80000;
x3(1) = 2;
x4(1) = 1;
x5(1) = 0;
x6(1) = 0;
x7(1) = 0;
% INITIAL GUESS FOR OPTIMAL CONTROL INPUT
u1 = zeros(1,M+1); % Control input for government intervention
u2 = zeros(1,M+1); % Control input for testing
u3 = zeros(1,M+1); % Control input for vaccinating
% INITIAL CONDITIONS ADOINT SYSTEM
L1 = zeros(1,M+1);
L2 = zeros(1,M+1);
L3 = zeros(1,M+1);
L4 = zeros(1,M+1);
L5 = zeros(1,M+1);
L6 = zeros(1,M+1);
L7 = zeros(1,M+1);
L1(M+1) = 0;
L2(M+1) = 0;
L3(M+1) = 0;
L4(M+1) = 0;
L5(M+1) = 0;
L6(M+1) = 0;
L7(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldu2 = u2;
oldu3 = u3;
oldx1 = x1;
oldx2 = x2;
oldx3 = x3;
oldx4 = x4;
oldx5 = x5;
oldx6 = x6;
oldx7 = x7;
oldL1 = L1;
oldL2 = L2;
oldL3 = L3;
oldL4 = L4;
oldL5 = L5;
oldL6 = L6;
oldL7 = L7;
% SYSTEM DYNAMICCS
for i = 1:M
m11 = pi -beta*(1-u1(i))*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -mu*x1(i);
m12 = beta*(1-u1(i))*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -(omega+mu)*x2(i);
m13 = theta*omega*x2(i)-(lamdaa+u2(i)+gammaa+mua+mu)*x3(i);
m14 = (1-theta)*omega*x2(i)-(lamdas+u3(i)+etas+u3(i)+mu)*x4(i);
m15 = (lamdaa+u2(i))*x3(i)+(lamdas+u3(i))*x4(i)-(etaq+gammaq+mu)*x5(i);
m16 = (etas+u3(i))*x4(i)+etaq*x5(i)- (gammah+muh+mu)*x6(i);
m17 = gammaa*x3(i) + gammaq*x4(i) + gammah*x6(i)-mu*x7(i);
x1(i+1) = x1(i) + h*m11;
x2(i+1) = x2(i) + h*m12;
x3(i+1) = x3(i) + h*m13;
x4(i+1) = x4(i) + h*m14;
x5(i+1) = x5(i) + h*m15;
x6(i+1) = x6(i) + h*m16;
x7(i+1) = x7(i) + h*m17;
end
% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 - i; % From final value to initial value
n11 = (L1(j)-L2(j))*beta*(1-u1(j))*(zetaa*x3(i)+zetas*x4(i)+zetaq*x5(i))*(1/N) + L1(j)*u1(j);
n12 = L2(j)*(omega+mu) - L3(j)*theta*omega - L4(j)*(1-theta)*omega;
n13 = -c1+(L1(j)-L2(j))*beta*(1-u1(j))*zetaa*(x1(i)/N) + (L3(j)-L5(j))*(lamdaa+u2(j)) + gammaa*(L3(j)-L7(j)) + L3(j)*(mua+mu);
n14 = -c2+(L1(j)-L2(j))*beta*(1-u1(j))*zetas*(x1(i)/N)+ (L4(j)-L5(j))*(lamdas+u3(j))+ (L4(j)-L6(j))*(etas+u3(j))+ L4(j)*mu;
n15 = -c3+(L1(j)-L2(j))*beta*(1-u1(j))*zetaq*(x1(i)/N)+(L5(j)-L6(j))*etaq + (L5(j)-L7(j))*gammaq +L5(j)*mu;
n16 = -c4 + (L6(j)-L7(j))*gammah + L6(j)*(muh+mu);
n17 = L7(j)*mu;
L1(j-1) = L1(j) - h*n11;
L2(j-1) = L2(j) - h*n12;
L3(j-1) = L3(j) - h*n13;
L4(j-1) = L4(j) - h*n14;
L5(j-1) = L5(j) - h*n15;
L6(j-1) = L6(j) - h*n16;
L7(j-1) = L7(j) - h*n17;
end
% OPTIMALITY CONDITIONS
U1 = min(1,max(0,(L2-L1).*beta.*((zetaa.*x3(i)+zetas.*x4(i)+zetaq.*x5(i)).*x1(i))./(c5.*N)));
u1 = 0.01.*U1 +0.99.*oldu1;
U2 = min(1, max(0, (((L5-L1).*x3)./(c6))));
u2 = 0.01.*U2 +0.99.*oldu2;
U3 = min(1, max(0, ((((L4-L5)+(L4-L6)).*x4)./(c7))));
%u3 = 0.01.*U3 +0.99.*oldu3;
u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1*x3+c2*x4+c3*x5+c4*x6+c5./2*sum(u1.^2)*h + c6./2*sum(u2.^2)*h + c7./2*sum(u3.^2)*h;
% Cost1 = c1./2.*cumsum(x4.^2)*h; % Total cost of threatened population
% Cost2 = b1./2.*cumsum(u1.^2)*h; % Total cost of control input u1
% Cost3 = b2./2.*cumsum(u2.^2)*h; % Total cost of control input u2
% Cost4 = b2./2.*cumsum(u3.^2)*h; % Total cost of control input u3
% Cost5 = c2.*x6; % Total cost of deceased population
% J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5; % Cost at each time for ...plotting graphs
% CHECK CONVERGENCE TO STOP SWEEP METHOD
temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
temp4 = delta*sum(abs(x1)) - sum(abs(oldx1 - x1));
temp5 = delta*sum(abs(x2)) - sum(abs(oldx2 - x2));
temp6 = delta*sum(abs(x3)) - sum(abs(oldx3 - x3));
temp7 = delta*sum(abs(x4)) - sum(abs(oldx4 - x4));
temp8 = delta*sum(abs(x5)) - sum(abs(oldx5 - x5));
temp9 = delta*sum(abs(x6)) - sum(abs(oldx6 - x6));
temp10 = delta*sum(abs(x7)) - sum(abs(oldx7 - x7));
temp11 = delta*sum(abs(L1)) - sum(abs(oldL1 - L1));
temp12 = delta*sum(abs(L2)) - sum(abs(oldL2 - L2));
temp13 = delta*sum(abs(L3)) - sum(abs(oldL3 - L3));
temp14 = delta*sum(abs(L4)) - sum(abs(oldL4 - L4));
temp15 = delta*sum(abs(L5)) - sum(abs(oldL5 - L5));
temp16 = delta*sum(abs(L6)) - sum(abs(oldL6 - L6));
temp17 = delta*sum(abs(L7)) - sum(abs(oldL7 - L7));
%test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11 temp12 temp13 temp14 temp15 temp16 temp17]);
test = min(temp1,min(temp2,min(temp3,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,temp17))))))))))))))));
end
% disp(['number of loops: ' num2str(loopcnt)]);
% disp(['Cost function: ' num2str(J)]);
% disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:) = t;
y(2,:) = x1;
y(3,:) = x2;
y(4,:) = x3;
y(5,:) = x4;
y(6,:) = x5;
y(7,:) = x6;
y(8,:) = x7;
y(9,:) = L1;
y(10,:) = L2;
y(11,:) = L3;
y(12,:) = L4;
y(13,:) = L5;
y(14,:) = L6;
y(15,:) = L7;
y(16,:) = u1;
y(17,:) = u2;
y(18,:) = u3;
y(19,:) = J;
s0 = 1280000000;
e0 = 10000;
a0 = 2;
i0 = 1;
q0 = 0;
h0 = 0;
r0 = 0;
%d0 =0;
y0 = [s0;e0;a0;i0;q0;h0;r0];
tspan = [0,400];
p = [beta; lamdaa; lamdas; etas; etaq; gammaa; gammaq; gammah ];
[t1,y] = ode45(@immanuel, tspan, y0, [], p);
figure(1)
hold on
plot(t,x4,'r-','linewidth',2)
plot(t1,y(:,4),'b-','linewidth',2)
xlabel('Time(days)');
ylabel(' population');
legend('with control','Without control')
box on
hold on
code for control u1 (u2=0 and u3=0 in code1)
clear all; close all;
% SETUP FOR FORWARD-BACKWARD SWEEP METHOD
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 400; % Final time
delta = 0.0001; % Accepted tollerance
M = 999; % Number of nodes
t = linspace(t0,tf,M+1); % Time variable where linspace creates M+1 ...equally spaced nodes between t0 and tf, including t0 and tf.
h = tf/M; % Spacing between nodes
h2 = h/2; % Spacing equal to 2 for Runge-Kutta method
% MODEL PARAMETERS
N = 1390000000;
pi = 70 ;
zetaa = 0.5;
zetas = 0.3;
zetaq = 0.3;
omega = 0.2;
theta = 0.5;
mu = 0.00391;
mua = 0.001945;
muh = 0.001945;
beta= 0.5888;
lamdaa= 0.3679;
lamdas= 0.2393;
etas= 0.1049;
etaq= 0.1917;
gammaa=0.2524;
gammaq= 0.1050;
gammah=0.0930;
% WEIGHT FACTORS
c1 = 1;
c2 = 2;
c3 = 2;
c4 = 2;
c5 = 30;
c6 = 30;
c7 = 30;
% INITIAL CONDITIONS MODEL
x1=zeros(1,M+1);
x2=zeros(1,M+1);
x3=zeros(1,M+1);
x4=zeros(1,M+1);
x5=zeros(1,M+1);
x6=zeros(1,M+1);
x7=zeros(1,M+1);
x1(1) = 13900000;
x2(1) = 80000;
x3(1) = 2;
x4(1) = 1;
x5(1) = 0;
x6(1) = 0;
x7(1) = 0;
% INITIAL GUESS FOR OPTIMAL CONTROL INPUT
u1 = zeros(1,M+1); % Control input for government intervention
% Control input for vaccinating
% INITIAL CONDITIONS ADOINT SYSTEM
L1 = zeros(1,M+1);
L2 = zeros(1,M+1);
L3 = zeros(1,M+1);
L4 = zeros(1,M+1);
L5 = zeros(1,M+1);
L6 = zeros(1,M+1);
L7 = zeros(1,M+1);
L1(M+1) = 0;
L2(M+1) = 0;
L3(M+1) = 0;
L4(M+1) = 0;
L5(M+1) = 0;
L6(M+1) = 0;
L7(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldx1 = x1;
oldx2 = x2;
oldx3 = x3;
oldx4 = x4;
oldx5 = x5;
oldx6 = x6;
oldx7 = x7;
oldL1 = L1;
oldL2 = L2;
oldL3 = L3;
oldL4 = L4;
oldL5 = L5;
oldL6 = L6;
oldL7 = L7;
% SYSTEM DYNAMICCS
for i = 1:M
m11 = pi -beta*(1-u1(i))*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -mu*x1(i);
m12 = beta*(1-u1(i))*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -(omega+mu)*x2(i);
m13 = theta*omega*x2(i)-(lamdaa+gammaa+mua+mu)*x3(i);
m14 = (1-theta)*omega*x2(i)-(lamdas+etas+mu)*x4(i);
m15 = (lamdaa)*x3(i)+(lamdas)*x4(i)-(etaq+gammaq+mu)*x5(i);
m16 = (etas)*x4(i)+etaq*x5(i)- (gammah+muh+mu)*x6(i);
m17 = gammaa*x3(i) + gammaq*x4(i) + gammah*x6(i)-mu*x7(i);
x1(i+1) = x1(i) + h*m11;
x2(i+1) = x2(i) + h*m12;
x3(i+1) = x3(i) + h*m13;
x4(i+1) = x4(i) + h*m14;
x5(i+1) = x5(i) + h*m15;
x6(i+1) = x6(i) + h*m16;
x7(i+1) = x7(i) + h*m17;
end
% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 - i; % From final value to initial value
n11 = (L1(j)-L2(j))*beta*(1-u1(j))*(zetaa*x3(i)+zetas*x4(i)+zetaq*x5(i))*(1/N) + L1(j)*u1(j);
n12 = L2(j)*(omega+mu) - L3(j)*theta*omega - L4(j)*(1-theta)*omega;
n13 = -c1+(L1(j)-L2(j))*beta*(1-u1(j))*zetaa*(x1(i)/N) + (L3(j)-L5(j))*(lamdaa) + gammaa*(L3(j)-L7(j)) + L3(j)*(mua+mu);
n14 = -c2+(L1(j)-L2(j))*beta*(1-u1(j))*zetas*(x1(i)/N)+ (L4(j)-L5(j))*(lamdas)+ (L4(j)-L6(j))*(etas)+ L4(j)*mu;
n15 = -c3+(L1(j)-L2(j))*beta*(1-u1(j))*zetaq*(x1(i)/N)+(L5(j)-L6(j))*etaq + (L5(j)-L7(j))*gammaq +L5(j)*mu;
n16 = -c4 + (L6(j)-L7(j))*gammah + L6(j)*(muh+mu);
n17 = L7(j)*mu;
L1(j-1) = L1(j) - h*n11;
L2(j-1) = L2(j) - h*n12;
L3(j-1) = L3(j) - h*n13;
L4(j-1) = L4(j) - h*n14;
L5(j-1) = L5(j) - h*n15;
L6(j-1) = L6(j) - h*n16;
L7(j-1) = L7(j) - h*n17;
end
% OPTIMALITY CONDITIONS
U1 = min(1,max(0,(L2-L1).*beta.*((zetaa.*x3(i)+zetas.*x4(i)+zetaq.*x5(i)).*x1(i))./(c5.*N)));
u1 = 0.01.*U1 +0.99.*oldu1;
% U2 = min(1, max(0, (((L5-L1).*x3)./(c6))));
% u2 = 0.01.*U2 +0.99.*oldu2;
%
% U3 = min(1, max(0, ((((L4-L5)+(L4-L6)).*x4)./(c7))));
% %u3 = 0.01.*U3 +0.99.*oldu3;
% u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1*x3+c2*x4+c3*x5+c4*x6+c5./2*sum(u1.^2)*h ;
% Cost1 = c1./2.*cumsum(x4.^2)*h; % Total cost of threatened population
% Cost2 = b1./2.*cumsum(u1.^2)*h; % Total cost of control input u1
% Cost3 = b2./2.*cumsum(u2.^2)*h; % Total cost of control input u2
% Cost4 = b2./2.*cumsum(u3.^2)*h; % Total cost of control input u3
% Cost5 = c2.*x6; % Total cost of deceased population
% J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5; % Cost at each time for ...plotting graphs
% CHECK CONVERGENCE TO STOP SWEEP METHOD
temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
%temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
%temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
temp4 = delta*sum(abs(x1)) - sum(abs(oldx1 - x1));
temp5 = delta*sum(abs(x2)) - sum(abs(oldx2 - x2));
temp6 = delta*sum(abs(x3)) - sum(abs(oldx3 - x3));
temp7 = delta*sum(abs(x4)) - sum(abs(oldx4 - x4));
temp8 = delta*sum(abs(x5)) - sum(abs(oldx5 - x5));
temp9 = delta*sum(abs(x6)) - sum(abs(oldx6 - x6));
temp10 = delta*sum(abs(x7)) - sum(abs(oldx7 - x7));
temp11 = delta*sum(abs(L1)) - sum(abs(oldL1 - L1));
temp12 = delta*sum(abs(L2)) - sum(abs(oldL2 - L2));
temp13 = delta*sum(abs(L3)) - sum(abs(oldL3 - L3));
temp14 = delta*sum(abs(L4)) - sum(abs(oldL4 - L4));
temp15 = delta*sum(abs(L5)) - sum(abs(oldL5 - L5));
temp16 = delta*sum(abs(L6)) - sum(abs(oldL6 - L6));
temp17 = delta*sum(abs(L7)) - sum(abs(oldL7 - L7));
%test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11 temp12 temp13 temp14 temp15 temp16 temp17]);
test = min(temp1,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,temp17))))))))))))));
end
% disp(['number of loops: ' num2str(loopcnt)]);
% disp(['Cost function: ' num2str(J)]);
% disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:) = t;
y(2,:) = x1;
y(3,:) = x2;
y(4,:) = x3;
y(5,:) = x4;
y(6,:) = x5;
y(7,:) = x6;
y(8,:) = x7;
y(9,:) = L1;
y(10,:) = L2;
y(11,:) = L3;
y(12,:) = L4;
y(13,:) = L5;
y(14,:) = L6;
y(15,:) = L7;
y(16,:) = u1;
%y(17,:) = u2;
%y(18,:) = u3;
y(19,:) = J;
s0 = 1280000000;
e0 = 10000;
a0 = 2;
i0 = 1;
q0 = 0;
h0 = 0;
r0 = 0;
%d0 =0;
y0 = [s0;e0;a0;i0;q0;h0;r0];
tspan = [0,400];
p = [beta; lamdaa; lamdas; etas; etaq; gammaa; gammaq; gammah ];
[t1,y] = ode45(@immanuel, tspan, y0, [], p);
figure(2)
hold on
plot(t,x4,'r-','linewidth',2)
%plot(t1,y(:,4),'b-','linewidth',2)
xlabel('Time(days)');
ylabel(' population');
legend('u1 only')
box on
hold on
code for control u2 (u1=0 and u3=0 in code1)
clear all; close all;
% SETUP FOR FORWARD-BACKWARD SWEEP METHOD
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 400; % Final time
delta = 0.0001; % Accepted tollerance
M = 999; % Number of nodes
t = linspace(t0,tf,M+1); % Time variable where linspace creates M+1 ...equally spaced nodes between t0 and tf, including t0 and tf.
h = tf/M; % Spacing between nodes
h2 = h/2; % Spacing equal to 2 for Runge-Kutta method
% MODEL PARAMETERS
N = 1390000000;
pi = 70 ;
zetaa = 0.5;
zetas = 0.3;
zetaq = 0.3;
omega = 0.2;
theta = 0.5;
mu = 0.00391;
mua = 0.001945;
muh = 0.001945;
beta= 0.5888;
lamdaa= 0.3679;
lamdas= 0.2393;
etas= 0.1049;
etaq= 0.1917;
gammaa=0.2524;
gammaq= 0.1050;
gammah=0.0930;
% WEIGHT FACTORS
c1 = 1;
c2 = 2;
c3 = 2;
c4 = 2;
c5 = 30;
c6 = 30;
c7 = 30;
% INITIAL CONDITIONS MODEL
x1=zeros(1,M+1);
x2=zeros(1,M+1);
x3=zeros(1,M+1);
x4=zeros(1,M+1);
x5=zeros(1,M+1);
x6=zeros(1,M+1);
x7=zeros(1,M+1);
x1(1) = 13900000;
x2(1) = 80000;
x3(1) = 2;
x4(1) = 1;
x5(1) = 0;
x6(1) = 0;
x7(1) = 0;
% INITIAL GUESS FOR OPTIMAL CONTROL INPUT
u2 = zeros(1,M+1); % Control input for testing
% INITIAL CONDITIONS ADOINT SYSTEM
L1 = zeros(1,M+1);
L2 = zeros(1,M+1);
L3 = zeros(1,M+1);
L4 = zeros(1,M+1);
L5 = zeros(1,M+1);
L6 = zeros(1,M+1);
L7 = zeros(1,M+1);
L1(M+1) = 0;
L2(M+1) = 0;
L3(M+1) = 0;
L4(M+1) = 0;
L5(M+1) = 0;
L6(M+1) = 0;
L7(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu2 = u2;
% oldu = u3;
oldx1 = x1;
oldx2 = x2;
oldx3 = x3;
oldx4 = x4;
oldx5 = x5;
oldx6 = x6;
oldx7 = x7;
oldL1 = L1;
oldL2 = L2;
oldL3 = L3;
oldL4 = L4;
oldL5 = L5;
oldL6 = L6;
oldL7 = L7;
% SYSTEM DYNAMICCS
for i = 1:M
m11 = pi -beta*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -mu*x1(i);
m12 = beta*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -(omega+mu)*x2(i);
m13 = theta*omega*x2(i)-(lamdaa+u2(i)+gammaa+mua+mu)*x3(i);
m14 = (1-theta)*omega*x2(i)-(lamdas+etas+mu)*x4(i);
m15 = (lamdaa+u2(i))*x3(i)+(lamdas)*x4(i)-(etaq+gammaq+mu)*x5(i);
m16 = (etas)*x4(i)+etaq*x5(i)- (gammah+muh+mu)*x6(i);
m17 = gammaa*x3(i) + gammaq*x4(i) + gammah*x6(i)-mu*x7(i);
x1(i+1) = x1(i) + h*m11;
x2(i+1) = x2(i) + h*m12;
x3(i+1) = x3(i) + h*m13;
x4(i+1) = x4(i) + h*m14;
x5(i+1) = x5(i) + h*m15;
x6(i+1) = x6(i) + h*m16;
x7(i+1) = x7(i) + h*m17;
end
% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 - i; % From final value to initial value
n11 = (L1(j)-L2(j))*beta*(zetaa*x3(i)+zetas*x4(i)+zetaq*x5(i))*(1/N) ;
n12 = L2(j)*(omega+mu) - L3(j)*theta*omega - L4(j)*(1-theta)*omega;
n13 = -c1+(L1(j)-L2(j))*beta*zetaa*(x1(i)/N) + (L3(j)-L5(j))*(lamdaa+u2(j)) + gammaa*(L3(j)-L7(j)) + L3(j)*(mua+mu);
n14 = -c2+(L1(j)-L2(j))*beta*zetas*(x1(i)/N)+ (L4(j)-L5(j))*(lamdas)+ (L4(j)-L6(j))*(etas)+ L4(j)*mu;
n15 = -c3+(L1(j)-L2(j))*beta*zetaq*(x1(i)/N)+(L5(j)-L6(j))*etaq + (L5(j)-L7(j))*gammaq +L5(j)*mu;
n16 = -c4 + (L6(j)-L7(j))*gammah + L6(j)*(muh+mu);
n17 = L7(j)*mu;
L1(j-1) = L1(j) - h*n11;
L2(j-1) = L2(j) - h*n12;
L3(j-1) = L3(j) - h*n13;
L4(j-1) = L4(j) - h*n14;
L5(j-1) = L5(j) - h*n15;
L6(j-1) = L6(j) - h*n16;
L7(j-1) = L7(j) - h*n17;
end
% OPTIMALITY CONDITIONS
% U1 = min(1,max(0,(L2-L1).*beta.*((zetaa.*x3(i)+zetas.*x4(i)+zetaq.*x5(i)).*x1(i))./(c5.*N)));
% u1 = 0.01.*U1 +0.99.*oldu1;
U2 = min(1, max(0, (((L5-L1).*x3)./(c6))));
u2 = 0.01.*U2 +0.99.*oldu2;
% U3 = min(1, max(0, ((((L4-L5)+(L4-L6)).*x4)./(c7))));
% %u3 = 0.01.*U3 +0.99.*oldu3;
% u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1*x3+c2*x4+c3*x5+c4*x6 +c6./2*sum(u2.^2)*h;
% Cost1 = c1./2.*cumsum(x4.^2)*h; % Total cost of threatened population
% Cost2 = b1./2.*cumsum(u1.^2)*h; % Total cost of control input u1
% Cost3 = b2./2.*cumsum(u2.^2)*h; % Total cost of control input u2
% Cost4 = b2./2.*cumsum(u3.^2)*h; % Total cost of control input u3
% Cost5 = c2.*x6; % Total cost of deceased population
% J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5; % Cost at each time for ...plotting graphs
% CHECK CONVERGENCE TO STOP SWEEP METHOD
%temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
%temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
temp4 = delta*sum(abs(x1)) - sum(abs(oldx1 - x1));
temp5 = delta*sum(abs(x2)) - sum(abs(oldx2 - x2));
temp6 = delta*sum(abs(x3)) - sum(abs(oldx3 - x3));
temp7 = delta*sum(abs(x4)) - sum(abs(oldx4 - x4));
temp8 = delta*sum(abs(x5)) - sum(abs(oldx5 - x5));
temp9 = delta*sum(abs(x6)) - sum(abs(oldx6 - x6));
temp10 = delta*sum(abs(x7)) - sum(abs(oldx7 - x7));
temp11 = delta*sum(abs(L1)) - sum(abs(oldL1 - L1));
temp12 = delta*sum(abs(L2)) - sum(abs(oldL2 - L2));
temp13 = delta*sum(abs(L3)) - sum(abs(oldL3 - L3));
temp14 = delta*sum(abs(L4)) - sum(abs(oldL4 - L4));
temp15 = delta*sum(abs(L5)) - sum(abs(oldL5 - L5));
temp16 = delta*sum(abs(L6)) - sum(abs(oldL6 - L6));
temp17 = delta*sum(abs(L7)) - sum(abs(oldL7 - L7));
%test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11 temp12 temp13 temp14 temp15 temp16 temp17]);
test = min(temp2,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,temp17))))))))))))));
end
% disp(['number of loops: ' num2str(loopcnt)]);
% disp(['Cost function: ' num2str(J)]);
% disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:) = t;
y(2,:) = x1;
y(3,:) = x2;
y(4,:) = x3;
y(5,:) = x4;
y(6,:) = x5;
y(7,:) = x6;
y(8,:) = x7;
y(9,:) = L1;
y(10,:) = L2;
y(11,:) = L3;
y(12,:) = L4;
y(13,:) = L5;
y(14,:) = L6;
y(15,:) = L7;
%y(16,:) = u1;
y(17,:) = u2;
%y(18,:) = u3;
y(19,:) = J;
a0 = 1280000000;
b0 = 10000;
c0 = 2;
d0 = 1;
e0 = 0;
f0 = 0;
h0 = 0;
%d0 =0;
y0 = [a0;b0;c0;d0;e0;f0;g0];
tspan = [0,400];
p = [beta; lamdaa; lamdas; etas; etaq; gammaa; gammaq; gammah ];
[t1,y] = ode45(@immanuel, tspan, y0, [], p);
figure(3)
hold on
plot(t,x4,'r-','linewidth',2)
% plot(t1,y(:,4),'b-','linewidth',2)
xlabel('Time(days)');
ylabel(' population');
legend('u2 only')
box on
hold on
figures are
I need like this figure

Sign in to answer this question.

Answers (2)

Torsten
Torsten on 22 Nov 2022
Ran in:
Use "hold on" and "hold off" to plot several graphs in one figure:
f = @(x) x.^2;
g = @(x)sqrt(x);
x = 0:0.01:1;
hold on
plot(x,f(x))
plot(x,g(x))
hold off
grid on
  3 Comments
Show 1 older comment
mallela ankamma rao
mallela ankamma rao on 23 Nov 2022
Torsten sir i request you please tell me how to minimize code because in every code half of the content was same
rachel adi
rachel adi on 8 May 2023
Edited: rachel adi on 8 May 2023
Hi @mallela ankamma rao,
I am currently learning about optimal control for epidemic models with more than one control or intervention involved. If I may, I would like to have the model and Matlab code you implemented here.
I would be very grateful if you could send it to me at adi.rachel17@gmail.com.
Best wishes,
rachel.

Sign in to comment.

Pavl M.
Pavl M. on 21 Nov 2024 at 8:45
Ran in:
Code adjustment (descreasing, length, beautifying, refromatting, refactoring):
I detected in the programme code script that
the setup
SETUP FOR FORWARD-BACKWARD SWEEP METHOD
repeated 3 times, so refactored it
clear all; close all;
% SETUP FOR FORWARD-BACKWARD SWEEP METHOD
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 400; % Final time
delta = 0.0001; % Accepted tollerance
M = 999; % Number of nodes
t = linspace(t0,tf,M+1); % Time variable where linspace creates M+1 ...equally spaced nodes between t0 and tf, including t0 and tf.
h = tf/M; % Spacing between nodes
h2 = h/2; % Spacing equal to 2 for Runge-Kutta method
% MODEL PARAMETERS
N = 1390000000;
pi = 70 ;
zetaa = 0.5;
zetas = 0.3;
zetaq = 0.3;
omega = 0.2;
theta = 0.5;
mu = 0.00391;
mua = 0.001945;
muh = 0.001945;
beta= 0.5888;
lamdaa= 0.3679;
lamdas= 0.2393;
etas= 0.1049;
etaq= 0.1917;
gammaa=0.2524;
gammaq= 0.1050;
gammah=0.0930;
% WEIGHT FACTORS
c1 = 1;
c2 = 2;
c3 = 2;
c4 = 2;
c5 = 30;
c6 = 30;
c7 = 30;
% INITIAL CONDITIONS MODEL
x1=zeros(1,M+1); % Susceptible
x2=zeros(1,M+1); % exposed
x3=zeros(1,M+1); % Asymptomatic Infected
x4=zeros(1,M+1); % symptomatic Infected
x5=zeros(1,M+1); % quarantined
x6=zeros(1,M+1); % hospitlized
x7=zeros(1,M+1); % recovered
x1(1) = 1390000000;
x2(1) = 80000;
x3(1) = 2;
x4(1) = 1;
x5(1) = 0;
x6(1) = 0;
x7(1) = 0;
% INITIAL GUESS FOR OPTIMAL CONTROL INPUT
u1 = zeros(1,M+1); % Control input for government intervention
u2 = zeros(1,M+1); % Control input for testing
u3 = zeros(1,M+1); % Control input for vaccinating
% INITIAL CONDITIONS ADJOINT SYSTEM
L1 = zeros(1,M+1);
L2 = zeros(1,M+1);
L3 = zeros(1,M+1);
L4 = zeros(1,M+1);
L5 = zeros(1,M+1);
L6 = zeros(1,M+1);
L7 = zeros(1,M+1);
L1(M+1) = 0;
L2(M+1) = 0;
L3(M+1) = 0;
L4(M+1) = 0;
L5(M+1) = 0;
L6(M+1) = 0;
L7(M+1) = 0;
%where repeated three times so far
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldu2 = u2;
oldu3 = u3;
oldx1 = x1;
oldx2 = x2;
oldx3 = x3;
oldx4 = x4;
oldx5 = x5;
oldx6 = x6;
oldx7 = x7;
oldL1 = L1;
oldL2 = L2;
oldL3 = L3;
oldL4 = L4;
oldL5 = L5;
oldL6 = L6;
oldL7 = L7;
% SYSTEM DYNAMICS
for i = 1:M
m11 = pi -beta*(1-u1(i))*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -mu*x1(i);
m12 = beta*(1-u1(i))*(zetaa*x3(i) +zetas*x4(i)+zetaq*x5(i))*(x1(i)/N) -(omega+mu)*x2(i);
m13 = theta*omega*x2(i)-(lamdaa+u2(i)+gammaa+mua+mu)*x3(i);
m14 = (1-theta)*omega*x2(i)-(lamdas+u3(i)+etas+u3(i)+mu)*x4(i);
m15 = (lamdaa+u2(i))*x3(i)+(lamdas+u3(i))*x4(i)-(etaq+gammaq+mu)*x5(i);
m16 = (etas+u3(i))*x4(i)+etaq*x5(i)- (gammah+muh+mu)*x6(i);
m17 = gammaa*x3(i) + gammaq*x4(i) + gammah*x6(i)-mu*x7(i);
x1(i+1) = x1(i) + h*m11;
x2(i+1) = x2(i) + h*m12;
x3(i+1) = x3(i) + h*m13;
x4(i+1) = x4(i) + h*m14;
x5(i+1) = x5(i) + h*m15;
x6(i+1) = x6(i) + h*m16;
x7(i+1) = x7(i) + h*m17;
end
% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 - i; % From final value to initial value
n11 = (L1(j)-L2(j))*beta*(1-u1(j))*(zetaa*x3(i)+zetas*x4(i)+zetaq*x5(i))*(1/N) + L1(j)*u1(j);
n12 = L2(j)*(omega+mu) - L3(j)*theta*omega - L4(j)*(1-theta)*omega;
n13 = -c1+(L1(j)-L2(j))*beta*(1-u1(j))*zetaa*(x1(i)/N) + (L3(j)-L5(j))*(lamdaa+u2(j)) + gammaa*(L3(j)-L7(j)) + L3(j)*(mua+mu);
n14 = -c2+(L1(j)-L2(j))*beta*(1-u1(j))*zetas*(x1(i)/N)+ (L4(j)-L5(j))*(lamdas+u3(j))+ (L4(j)-L6(j))*(etas+u3(j))+ L4(j)*mu;
n15 = -c3+(L1(j)-L2(j))*beta*(1-u1(j))*zetaq*(x1(i)/N)+(L5(j)-L6(j))*etaq + (L5(j)-L7(j))*gammaq +L5(j)*mu;
n16 = -c4 + (L6(j)-L7(j))*gammah + L6(j)*(muh+mu);
n17 = L7(j)*mu;
L1(j-1) = L1(j) - h*n11;
L2(j-1) = L2(j) - h*n12;
L3(j-1) = L3(j) - h*n13;
L4(j-1) = L4(j) - h*n14;
L5(j-1) = L5(j) - h*n15;
L6(j-1) = L6(j) - h*n16;
L7(j-1) = L7(j) - h*n17;
end
% OPTIMALITY CONDITIONS
U1 = min(1,max(0,(L2-L1).*beta.*((zetaa.*x3(i)+zetas.*x4(i)+zetaq.*x5(i)).*x1(i))./(c5.*N)));
u1 = 0.01.*U1 +0.99.*oldu1;
U2 = min(1, max(0, (((L5-L1).*x3)./(c6))));
u2 = 0.01.*U2 +0.99.*oldu2;
U3 = min(1, max(0, ((((L4-L5)+(L4-L6)).*x4)./(c7))));
%u3 = 0.01.*U3 +0.99.*oldu3;
u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1*x3+c2*x4+c3*x5+c4*x6+c5./2*sum(u1.^2)*h + c6./2*sum(u2.^2)*h + c7./2*sum(u3.^2)*h;
% Cost1 = c1./2.*cumsum(x4.^2)*h; % Total cost of threatened population
% Cost2 = b1./2.*cumsum(u1.^2)*h; % Total cost of control input u1
% Cost3 = b2./2.*cumsum(u2.^2)*h; % Total cost of control input u2
% Cost4 = b2./2.*cumsum(u3.^2)*h; % Total cost of control input u3
% Cost5 = c2.*x6; % Total cost of deceased population
% J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5; % Cost at each time for ...plotting graphs
% CHECK CONVERGENCE TO STOP SWEEP METHOD
temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
temp4 = delta*sum(abs(x1)) - sum(abs(oldx1 - x1));
temp5 = delta*sum(abs(x2)) - sum(abs(oldx2 - x2));
temp6 = delta*sum(abs(x3)) - sum(abs(oldx3 - x3));
temp7 = delta*sum(abs(x4)) - sum(abs(oldx4 - x4));
temp8 = delta*sum(abs(x5)) - sum(abs(oldx5 - x5));
temp9 = delta*sum(abs(x6)) - sum(abs(oldx6 - x6));
temp10 = delta*sum(abs(x7)) - sum(abs(oldx7 - x7));
temp11 = delta*sum(abs(L1)) - sum(abs(oldL1 - L1));
temp12 = delta*sum(abs(L2)) - sum(abs(oldL2 - L2));
temp13 = delta*sum(abs(L3)) - sum(abs(oldL3 - L3));
temp14 = delta*sum(abs(L4)) - sum(abs(oldL4 - L4));
temp15 = delta*sum(abs(L5)) - sum(abs(oldL5 - L5));
temp16 = delta*sum(abs(L6)) - sum(abs(oldL6 - L6));
temp17 = delta*sum(abs(L7)) - sum(abs(oldL7 - L7));
%test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11 temp12 temp13 temp14 temp15 temp16 temp17]);
test = min(temp1,min(temp2,min(temp3,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,temp17))))))))))))))));
end
% disp(['number of loops: ' num2str(loopcnt)]);
% disp(['Cost function: ' num2str(J)]);
% disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:) = t;
y(2,:) = x1;
y(3,:) = x2;
y(4,:) = x3;
y(5,:) = x4;
y(6,:) = x5;
y(7,:) = x6;
y(8,:) = x7;
y(9,:) = L1;
y(10,:) = L2;
y(11,:) = L3;
y(12,:) = L4;
y(13,:) = L5;
y(14,:) = L6;
y(15,:) = L7;
y(16,:) = u1;
y(17,:) = u2;
y(18,:) = u3;
y(19,:) = J;
% code for u1
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldz1 = z1;
oldz2 = z2;
oldz3 = z3;
oldz4 = z4;
oldz5 = z5;
oldz6 = z6;
oldz7 = z7;
oldL1 = L1;
oldL2 = L2;
oldL3 = L3;
oldL4 = L4;
oldL5 = L5;
oldL6 = L6;
oldL7 = L7;
% SYSTEM DYNAMICS
for i = 1:M
m11 = pi -beta*(1-u1(i))*(zetaa*z3(i) +zetas*z4(i)+zetaq*z5(i))*(z1(i)/N) -mu*z1(i);
m12 = beta*(1-u1(i))*(zetaa*z3(i) +zetas*z4(i)+zetaq*z5(i))*(z1(i)/N) -(omega+mu)*z2(i);
m13 = theta*omega*z2(i)-(lamdaa+gammaa+mua+mu)*z3(i);
m14 = (1-theta)*omega*z2(i)-(lamdas+etas+mu)*z4(i);
m15 = (lamdaa)*z3(i)+(lamdas)*z4(i)-(etaq+gammaq+mu)*z5(i);
m16 = (etas)*z4(i)+etaq*z5(i)- (gammah+muh+mu)*z6(i);
m17 = gammaa*z3(i) + gammaq*z4(i) + gammah*z6(i)-mu*z7(i);
z1(i+1) = z1(i) + h*m11;
z2(i+1) = z2(i) + h*m12;
z3(i+1) = z3(i) + h*m13;
z4(i+1) = z4(i) + h*m14;
z5(i+1) = z5(i) + h*m15;
z6(i+1) = z6(i) + h*m16;
z7(i+1) = z7(i) + h*m17;
end
% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 - i; % From final value to initial value
n11 = (L1(j)-L2(j))*beta*(1-u1(j))*(zetaa*z3(i)+zetas*z4(i)+zetaq*z5(i))*(1/N) + L1(j)*u1(j);
n12 = L2(j)*(omega+mu) - L3(j)*theta*omega - L4(j)*(1-theta)*omega;
n13 = -c1+(L1(j)-L2(j))*beta*(1-u1(j))*zetaa*(z1(i)/N) + (L3(j)-L5(j))*(lamdaa) + gammaa*(L3(j)-L7(j)) + L3(j)*(mua+mu);
n14 = -c2+(L1(j)-L2(j))*beta*(1-u1(j))*zetas*(z1(i)/N)+ (L4(j)-L5(j))*(lamdas)+ (L4(j)-L6(j))*(etas)+ L4(j)*mu;
n15 = -c3+(L1(j)-L2(j))*beta*(1-u1(j))*zetaq*(z1(i)/N)+(L5(j)-L6(j))*etaq + (L5(j)-L7(j))*gammaq +L5(j)*mu;
n16 = -c4 + (L6(j)-L7(j))*gammah + L6(j)*(muh+mu);
n17 = L7(j)*mu;
L1(j-1) = L1(j) - h*n11;
L2(j-1) = L2(j) - h*n12;
L3(j-1) = L3(j) - h*n13;
L4(j-1) = L4(j) - h*n14;
L5(j-1) = L5(j) - h*n15;
L6(j-1) = L6(j) - h*n16;
L7(j-1) = L7(j) - h*n17;
end
% OPTIMALITY CONDITIONS
U1 = min(1,max(0,(L2-L1).*beta.*((zetaa.*z3(i)+zetas.*z4(i)+zetaq.*z5(i)).*z1(i))./(c5.*N)));
u1 = 0.01.*U1 +0.99.*oldu1;
% U2 = min(1, max(0, (((L5-L1).*x3)./(c6))));
% u2 = 0.01.*U2 +0.99.*oldu2;
%
% U3 = min(1, max(0, ((((L4-L5)+(L4-L6)).*x4)./(c7))));
% %u3 = 0.01.*U3 +0.99.*oldu3;
% u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1*z3+c2*z4+c3*z5+c4*z6+c5./2*sum(u1.^2)*h ;
% Cost1 = c1./2.*cumsum(x4.^2)*h; % Total cost of threatened population
% Cost2 = b1./2.*cumsum(u1.^2)*h; % Total cost of control input u1
% Cost3 = b2./2.*cumsum(u2.^2)*h; % Total cost of control input u2
% Cost4 = b2./2.*cumsum(u3.^2)*h; % Total cost of control input u3
% Cost5 = c2.*x6; % Total cost of deceased population
% J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5; % Cost at each time for ...plotting graphs
% CHECK CONVERGENCE TO STOP SWEEP METHOD
temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
%temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
%temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
temp4 = delta*sum(abs(z1)) - sum(abs(oldz1 - z1));
temp5 = delta*sum(abs(z2)) - sum(abs(oldz2 - z2));
temp6 = delta*sum(abs(z3)) - sum(abs(oldz3 - z3));
temp7 = delta*sum(abs(z4)) - sum(abs(oldz4 - z4));
temp8 = delta*sum(abs(z5)) - sum(abs(oldz5 - z5));
temp9 = delta*sum(abs(z6)) - sum(abs(oldz6 - z6));
temp10 = delta*sum(abs(z7)) - sum(abs(oldz7 - z7));
temp11 = delta*sum(abs(L1)) - sum(abs(oldL1 - L1));
temp12 = delta*sum(abs(L2)) - sum(abs(oldL2 - L2));
temp13 = delta*sum(abs(L3)) - sum(abs(oldL3 - L3));
temp14 = delta*sum(abs(L4)) - sum(abs(oldL4 - L4));
temp15 = delta*sum(abs(L5)) - sum(abs(oldL5 - L5));
temp16 = delta*sum(abs(L6)) - sum(abs(oldL6 - L6));
temp17 = delta*sum(abs(L7)) - sum(abs(oldL7 - L7));
%test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11 temp12 temp13 temp14 temp15 temp16 temp17]);
test = min(temp1,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,temp17))))))))))))));
% disp(['number of loops: ' num2str(loopcnt)]);
% disp(['Cost function: ' num2str(J)]);
% disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:) = t;
y(2,:) = z1;
y(3,:) = z2;
y(4,:) = z3;
y(5,:) = z4;
y(6,:) = z5;
y(7,:) = z6;
y(8,:) = z7;
y(9,:) = L1;
y(10,:) = L2;
y(11,:) = L3;
y(12,:) = L4;
y(13,:) = L5;
y(14,:) = L6;
y(15,:) = L7;
y(16,:) = u1;
%y(17,:) = u2;
%y(18,:) = u3;
y(19,:) = J;
end
% code for u2
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu2 = u2;
% oldu = u3;
oldw1 = w1;
oldw2 = w2;
oldw3 = w3;
oldw4 = w4;
oldw5 = w5;
oldw6 = w6;
oldw7 = w7;
oldL1 = L1;
oldL2 = L2;
oldL3 = L3;
oldL4 = L4;
oldL5 = L5;
oldL6 = L6;
oldL7 = L7;
% SYSTEM DYNAMICS
for i = 1:M
m11 = pi -beta*(zetaa*w3(i) +zetas*w4(i)+zetaq*w5(i))*(w1(i)/N) -mu*w1(i);
m12 = beta*(zetaa*w3(i) +zetas*w4(i)+zetaq*w5(i))*(w1(i)/N) -(omega+mu)*w2(i);
m13 = theta*omega*w2(i)-(lamdaa+u2(i)+gammaa+mua+mu)*w3(i);
m14 = (1-theta)*omega*w2(i)-(lamdas+etas+mu)*w4(i);
m15 = (lamdaa+u2(i))*w3(i)+(lamdas)*w4(i)-(etaq+gammaq+mu)*w5(i);
m16 = (etas)*w4(i)+etaq*w5(i)- (gammah+muh+mu)*w6(i);
m17 = gammaa*w3(i) + gammaq*w4(i) + gammah*w6(i)-mu*w7(i);
w1(i+1) = w1(i) + h*m11;
w2(i+1) = w2(i) + h*m12;
w3(i+1) = w3(i) + h*m13;
w4(i+1) = w4(i) + h*m14;
w5(i+1) = w5(i) + h*m15;
w6(i+1) = w6(i) + h*m16;
w7(i+1) = w7(i) + h*m17;
end
% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 - i; % From final value to initial value
n11 = (L1(j)-L2(j))*beta*(zetaa*w3(i)+zetas*w4(i)+zetaq*w5(i))*(1/N) ;
n12 = L2(j)*(omega+mu) - L3(j)*theta*omega - L4(j)*(1-theta)*omega;
n13 = -c1+(L1(j)-L2(j))*beta*zetaa*(w1(i)/N) + (L3(j)-L5(j))*(lamdaa+u2(j)) + gammaa*(L3(j)-L7(j)) + L3(j)*(mua+mu);
n14 = -c2+(L1(j)-L2(j))*beta*zetas*(w1(i)/N)+ (L4(j)-L5(j))*(lamdas)+ (L4(j)-L6(j))*(etas)+ L4(j)*mu;
n15 = -c3+(L1(j)-L2(j))*beta*zetaq*(w1(i)/N)+(L5(j)-L6(j))*etaq + (L5(j)-L7(j))*gammaq +L5(j)*mu;
n16 = -c4 + (L6(j)-L7(j))*gammah + L6(j)*(muh+mu);
n17 = L7(j)*mu;
L1(j-1) = L1(j) - h*n11;
L2(j-1) = L2(j) - h*n12;
L3(j-1) = L3(j) - h*n13;
L4(j-1) = L4(j) - h*n14;
L5(j-1) = L5(j) - h*n15;
L6(j-1) = L6(j) - h*n16;
L7(j-1) = L7(j) - h*n17;
end
% OPTIMALITY CONDITIONS
% U1 = min(1,max(0,(L2-L1).*beta.*((zetaa.*x3(i)+zetas.*x4(i)+zetaq.*x5(i)).*x1(i))./(c5.*N)));
% u1 = 0.01.*U1 +0.99.*oldu1;
U2 = min(1, max(0, (((L5-L1).*w3)./(c6))));
u2 = 0.01.*U2 +0.99.*oldu2;
% U3 = min(1, max(0, ((((L4-L5)+(L4-L6)).*x4)./(c7))));
% %u3 = 0.01.*U3 +0.99.*oldu3;
% u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1*w3+c2*w4+c3*w5+c4*w6 +c6./2*sum(u2.^2)*h;
% Cost1 = c1./2.*cumsum(x4.^2)*h; % Total cost of threatened population
% Cost2 = b1./2.*cumsum(u1.^2)*h; % Total cost of control input u1
% Cost3 = b2./2.*cumsum(u2.^2)*h; % Total cost of control input u2
% Cost4 = b2./2.*cumsum(u3.^2)*h; % Total cost of control input u3
% Cost5 = c2.*x6; % Total cost of deceased population
% J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5; % Cost at each time for ...plotting graphs
% CHECK CONVERGENCE TO STOP SWEEP METHOD
%temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
%temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
temp4 = delta*sum(abs(w1)) - sum(abs(oldw1 - w1));
temp5 = delta*sum(abs(w2)) - sum(abs(oldw2 - w2));
temp6 = delta*sum(abs(w3)) - sum(abs(oldw3 - w3));
temp7 = delta*sum(abs(w4)) - sum(abs(oldw4 - w4));
temp8 = delta*sum(abs(w5)) - sum(abs(oldw5 - w5));
temp9 = delta*sum(abs(w6)) - sum(abs(oldw6 - w6));
temp10 = delta*sum(abs(w7)) - sum(abs(oldw7 - w7));
temp11 = delta*sum(abs(L1)) - sum(abs(oldL1 - L1));
temp12 = delta*sum(abs(L2)) - sum(abs(oldL2 - L2));
temp13 = delta*sum(abs(L3)) - sum(abs(oldL3 - L3));
temp14 = delta*sum(abs(L4)) - sum(abs(oldL4 - L4));
temp15 = delta*sum(abs(L5)) - sum(abs(oldL5 - L5));
temp16 = delta*sum(abs(L6)) - sum(abs(oldL6 - L6));
temp17 = delta*sum(abs(L7)) - sum(abs(oldL7 - L7));
%test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11 temp12 temp13 temp14 temp15 temp16 temp17]);
test = min(temp2,min(temp4,min(temp5,min(temp6,min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,min(temp13,min(temp14,min(temp15,min(temp16,temp17))))))))))))));
% disp(['number of loops: ' num2str(loopcnt)]);
% disp(['Cost function: ' num2str(J)]);
% disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:) = t;
y(2,:) = w1;
y(3,:) = w2;
y(4,:) = w3;
y(5,:) = w4;
y(6,:) = w5;
y(7,:) = w6;
y(8,:) = w7;
y(9,:) = L1;
y(10,:) = L2;
y(11,:) = L3;
y(12,:) = L4;
y(13,:) = L5;
y(14,:) = L6;
y(15,:) = L7;
%y(16,:) = u1;
y(17,:) = u2;
%y(18,:) = u3;
y(19,:) = J;
end
a0 = 1280000000;
b0 = 10000;
c0 = 2;
d0 = 1;
e0 = 0;
f0 = 0;
h0 = 0;
g0 = 0;
%d0 =0;
y0 = [a0;b0;c0;d0;e0;f0;g0];
tspan = [0,400];
p = [beta; lamdaa; lamdas; etas; etaq; gammaa; gammaq; gammah ];
%[t1,y] = ode45(@immanuel, tspan, y0, [], p);
figure
%hold on
%plot(t1,y(:,4),'b-','linewidth',2)
plot(1:10:400,x4(1:10:400),'r-','linewidth',2)
figure
plot(1:10:400,y(1:10:400),'m-','linewidth',2)
%plot(t,w4/2,'g-','linewidth',2)
%xlabel('Time(days)');
ylabel(' population');
%legend('without controls','with controls','u1 only','u2 only')
%title('Symptomatic infected')
box off

Categories

Find more on Path Planning and Following in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!