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.