how to combine to get single graph in optimal control problems

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

[mallela ankamma rao]

i request you experts please anyone answer

Answer 1

[Torsten]

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

[mallela ankamma rao]

Thank you very much for your reply Torsten sir

As as you said, i have tried and i got graph sir

sir i request you please give me suggestions to decrease the code length

the code is

% code for with controls

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 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;

% code for u1

% 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

z1=zeros(1,M+1); % Susceptible

z2=zeros(1,M+1); % exposed

z3=zeros(1,M+1); % Asymptomatic Infected

z4=zeros(1,M+1); % symptomatic Infected

z5=zeros(1,M+1); % quarantined

z6=zeros(1,M+1); % hospitlized

z7=zeros(1,M+1); % recovered

z1(1) = 1390000000;

z2(1) = 80000;

z3(1) = 2;

z4(1) = 1;

z5(1) = 0;

z6(1) = 0;

z7(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;

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 DYNAMICCS

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))))))))))))));

end

% 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;

% code for u2

% 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

w1=zeros(1,M+1); % Susceptible

w2=zeros(1,M+1); % exposed

w3=zeros(1,M+1); % Asymptomatic Infected

w4=zeros(1,M+1); % symptomatic Infected

w5=zeros(1,M+1); % quarantined

w6=zeros(1,M+1); % hospitlized

w7=zeros(1,M+1); % recovered

w1(1) = 1390000000;

w2(1) = 80000;

w3(1) = 2;

w4(1) = 1;

w5(1) = 0;

w6(1) = 0;

w7(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;

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 DYNAMICCS

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))))))))))))));

end

% 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;

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(4)

hold on

plot(t1,y(:,4),'b-','linewidth',2)

plot(t,x4,'r-','linewidth',2)

plot(t,z4,'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 on

hold on

[mallela ankamma rao]

Torsten sir i request you please tell me how to minimize code because in every code half of the content was same

[rachel adi]

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 [e-mail].

Best wishes,

rachel.

Answer 2

[Pavl M.]

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