how to combine to get single graph in optimal control problems

4 views (last 30 days)
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

Answers (1)

Torsten
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
  2 Comments
mallela ankamma rao
mallela ankamma rao on 23 Nov 2022
Torsten sir i request you please tell me how to minimize code because in every code half of the content was same

Sign in to comment.

Categories

Find more on Simulink Functions 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!