how to combine to get single graph in optimal control problems

17 views (last 30 days)
mallela ankamma rao on 21 Nov 2022
Edited: rachel adi on 8 May 2023
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
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
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
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 on 22 Nov 2022
i request you experts please anyone answer

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

Categories

Find more on Get Started with Model-Based Calibration Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!