how to draw graphs for optimal control

31 views (last 30 days)
mallela ankamma rao
mallela ankamma rao on 25 Jul 2022
Answered: Pavl M. about 5 hours ago
good evening to one and all
sir i have optimal control matlab code for covid-19 model. but i dont know how to write graphs for using the code.i attached SIDERAV model pdf also please tell me anyone how to get graphs.
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 365; % Final time
Delta = 0.00001; % 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
gamma_i = 1/14; % Recovery rate undetected: 14 days
gamma_d = 1/14; % Recovery rate detected: 14 days
gamma_a = 1/12.4; % Recovery rate threatend: 12.4 days
beta = 0.251; % Infection rate
xi_i = 0.0053; % Rate infected undetected to acutely symtomatic
xi_d = 0.0053; % Rate infected detected to acutely symtomatic
mu = 0.0085; % Mortality rate
mu_hat = 5*mu; % Mortality rate when hospital capacity is exceeded
nu = 0.1; % Testing rate (no, slow, fast testing = 0, 0.05, 0.10)
h_bar = 0.00333; % Hospital capacity rate (0.00222, 0.00333, 0.00444)
psi = 0; % Rate of vaccination
% WEIGHT FACTORS
c1 = 0; % Weight on threatened population
c2 = 0; % Weigth on deceased population
c3 = 0;
b1 = 1;
b2 = 1;
b3 = 1;
% INITIAL CONDITIONS MODEL
x1=zeros(1,M+1); % Susceptible
x2=zeros(1,M+1); % Infected - undetected
x3=zeros(1,M+1); % Infected - detected
x4=zeros(1,M+1); % Acutely symptomatic - Threatened
x5=zeros(1,M+1); % Recovered
x6=zeros(1,M+1); % Deceased
x7=zeros(1,M+1); % Vaccinated
x1(1) = 1-0.00001;
x2(1) = 0.00001;
x3(1) = 0;
x4(1) = 0;
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) = c2;
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
% IMPACT HEALTHCARE CAPACITY ON MORTALITY RATE
if x4(i) < h_bar
mu_bar = mu*x4(i);
else
mu_bar = mu*h_bar + mu_hat*(x4(i) - h_bar);
end
m11 = -beta*x1(i)*x2(i)*(1-u1(i)) - psi*u3(i)*x1(i);
m12 = beta*x1(i)*x2(i)*(1-u1(i)) - gamma_i*x2(i) - xi_i*x2(i) - ...
nu*x2(i)*(u2(i));
m13 = nu*x2(i)*(u2(i))-gamma_d*x3(i)-xi_d*x3(i);
m14 = xi_i*x2(i)+xi_d*x3(i)-gamma_a*x4(i)-mu_bar;
m15 = gamma_i*x2(i) + gamma_d*x3(i) + gamma_a*x4(i);
m16 = mu_bar;
m17 = psi*u3(i)*x1(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*x2(j)*(1-u1(j)) + L7(j)*psi*u3(j);
n12 = (L1(j)-L2(j))*beta*x1(j)*(1-u1(j)) + L2(j)*(gamma_i + xi_i) + ...
L2(j)*(nu*(u2(j))) - L3(j)*(nu*(u2(j))) - L4(j)*xi_i;
n13 = L3(j)*(gamma_d + xi_d) - L4(j)*xi_d;
n14 = L4(j)*(gamma_a + mu_bar) - L5(j)*mu_bar - L6(j)*mu_bar - c1*x4(j);
n15 = 0;
n16 = 0;
n17 = 0;
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(0.8,max(0,(L2-L1).*beta.*x1.*x2./(b1)));
u1 = 0.01.*U1 +0.99.*oldu1;
U2 = min(1, max(0, (((L2-L3).*nu.*x2)./(b2))));
u2 = 0.01.*U2 +0.99.*oldu2;
U3 = min(1, max(0, (((L1-L7).*psi.*x1)./(b3))));
u3 = 0.01.*U3 +0.99.*oldu3;
% COST FUNCTION
J = c1./2*sum(x4.^2)*h + b1./2*sum(u1.^2)*h + b2./2*sum(u2.^2)*h + ...
b3./2*sum(u3.^2)*h+ c2*max(x6);
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]);
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;
% IMMUNITY REACHED
imm = x5+x7.*100; % Percentage immune
x4_per = x4*100; % percentage threatened
x6_per = x6*100; % percentage deceased
results
SIDAREV
number of loops: 2
Cost function: 0
Portion deceased: 0.016393
please help me how to draw graphs for u1,u2 ,u3,x(1),x(2),x(3),.with different values

Sign in to answer this question.

Answers (2)

Sam Chak
Sam Chak on 25 Jul 2022
Ran in:
Hi @mallela ankamma rao
Do expect plots like these?
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 365; % Final time
Delta = 0.00001; % 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
gamma_i = 1/14; % Recovery rate undetected: 14 days
gamma_d = 1/14; % Recovery rate detected: 14 days
gamma_a = 1/12.4; % Recovery rate threatend: 12.4 days
beta = 0.251; % Infection rate
xi_i = 0.0053; % Rate infected undetected to acutely symtomatic
xi_d = 0.0053; % Rate infected detected to acutely symtomatic
mu = 0.0085; % Mortality rate
mu_hat = 5*mu; % Mortality rate when hospital capacity is exceeded
nu = 0.1; % Testing rate (no, slow, fast testing = 0, 0.05, 0.10)
h_bar = 0.00333; % Hospital capacity rate (0.00222, 0.00333, 0.00444)
psi = 0; % Rate of vaccination
% WEIGHT FACTORS
c1 = 0; % Weight on threatened population
c2 = 0; % Weigth on deceased population
c3 = 0;
b1 = 1;
b2 = 1;
b3 = 1;
% INITIAL CONDITIONS MODEL
x1 = zeros(1, M+1); % Susceptible
x2 = zeros(1, M+1); % Infected - undetected
x3 = zeros(1, M+1); % Infected - detected
x4 = zeros(1, M+1); % Acutely symptomatic - Threatened
x5 = zeros(1, M+1); % Recovered
x6 = zeros(1, M+1); % Deceased
x7 = zeros(1, M+1); % Vaccinated
x1(1) = 1 - 0.00001;
x2(1) = 0.00001;
x3(1) = 0;
x4(1) = 0;
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) = c2;
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
% IMPACT HEALTHCARE CAPACITY ON MORTALITY RATE
if x4(i) < h_bar
mu_bar = mu*x4(i);
else
mu_bar = mu*h_bar + mu_hat*(x4(i) - h_bar);
end
m11 = -beta*x1(i)*x2(i)*(1-u1(i)) - psi*u3(i)*x1(i);
m12 = beta*x1(i)*x2(i)*(1-u1(i)) - gamma_i*x2(i) - xi_i*x2(i) - nu*x2(i)*(u2(i));
m13 = nu*x2(i)*(u2(i))-gamma_d*x3(i)-xi_d*x3(i);
m14 = xi_i*x2(i)+xi_d*x3(i)-gamma_a*x4(i)-mu_bar;
m15 = gamma_i*x2(i) + gamma_d*x3(i) + gamma_a*x4(i);
m16 = mu_bar;
m17 = psi*u3(i)*x1(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*x2(j)*(1-u1(j)) + L7(j)*psi*u3(j);
n12 = (L1(j)-L2(j))*beta*x1(j)*(1-u1(j)) + L2(j)*(gamma_i + xi_i) + L2(j)*(nu*(u2(j))) - L3(j)*(nu*(u2(j))) - L4(j)*xi_i;
n13 = L3(j)*(gamma_d + xi_d) - L4(j)*xi_d;
n14 = L4(j)*(gamma_a + mu_bar) - L5(j)*mu_bar - L6(j)*mu_bar - c1*x4(j);
n15 = 0;
n16 = 0;
n17 = 0;
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(0.8, max(0, (L2 - L1).*beta.*x1.*x2./(b1)));
u1 = 0.01.*U1 + 0.99.*oldu1;
U2 = min(1.0, max(0, (((L2 - L3).*nu.*x2)./(b2))));
u2 = 0.01.*U2 + 0.99.*oldu2;
U3 = min(1.0, max(0, (((L1 - L7).*psi.*x1)./(b3))));
u3 = 0.01.*U3 + 0.99.*oldu3;
% COST FUNCTION
J = c1./2*sum(x4.^2)*h + b1./2*sum(u1.^2)*h + b2./2*sum(u2.^2)*h + b3./2*sum(u3.^2)*h + c2*max(x6);
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]);
end
disp(['number of loops: ' num2str(loopcnt)]);
number of loops: 2
disp(['Cost function: ' num2str(J)]);
Cost function: 0
disp(['Portion deceased: ' num2str(max(x6))]);
Portion deceased: 0.016393
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;
% IMMUNITY REACHED
imm = x5 + x7.*100; % Percentage immune
x4_per = x4*100; % percentage threatened
x6_per = x6*100; % percentage deceased
plot(y(1,:), y(2,:)), grid on, xlabel('t'), ylabel('x_{1}') % plotting x1
plot(y(1,:), y(3,:)), grid on, xlabel('t'), ylabel('x_{2}') % plotting x2
plot(y(1,:), y(4,:)), grid on, xlabel('t'), ylabel('x_{3}') % plotting x3
  3 Comments
Show 1 older comment
Oladipupo Johnson
Oladipupo Johnson on 10 Apr 2023
Hi Everyone and Thanks to Mallela Ankamma Rao, your answers was helpful to our work. I also request you help with graph plot code for u1, u2 and u3 of SIDAREV model or kindly inbox me at oladipuposamueljohnson@gmail.com. Thank you Sir

Sign in to comment.

Pavl M.
Pavl M. 5 minutes ago
Ran in:
clc
clear all
close all
test = -1; % Test variable; as long as variable is negative ...the while loops keeps repeating
t0 = 0; % Initial time
tf = 365; % Final time
Delta = 0.00001; % 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
gamma_i = 1/14; % Recovery rate undetected: 14 days
gamma_d = 1/14; % Recovery rate detected: 14 days
gamma_a = 1/12.4; % Recovery rate threatend: 12.4 days
beta = 0.251; % Infection rate
xi_i = 0.0053; % Rate infected undetected to acutely symtomatic
xi_d = 0.0053; % Rate infected detected to acutely symtomatic
mu = 0.0085; % Mortality rate
mu_hat = 5*mu; % Mortality rate when hospital capacity is exceeded
nu = 0.1; % Testing rate (no, slow, fast testing = 0, 0.05, 0.10)
h_bar = 0.00333; % Hospital capacity rate (0.00222, 0.00333, 0.00444)
psi = 0; % Rate of vaccination
% WEIGHT FACTORS
c1 = 0; % Weight on threatened population
c2 = 0; % Weigth on deceased population
c3 = 0;
b1 = 1;
b2 = 1;
b3 = 1;
% INITIAL CONDITIONS MODEL
x1 = zeros(1, M+1); % Susceptible
x2 = zeros(1, M+1); % Infected - undetected
x3 = zeros(1, M+1); % Infected - detected
x4 = zeros(1, M+1); % Acutely symptomatic - Threatened
x5 = zeros(1, M+1); % Recovered
x6 = zeros(1, M+1); % Deceased
x7 = zeros(1, M+1); % Vaccinated
x1(1) = 1 - 0.00001;
x2(1) = 0.00001;
x3(1) = 0;
x4(1) = 0;
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) = c2;
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 DYNAMICS
for i = 1:M
% IMPACT HEALTHCARE CAPACITY ON MORTALITY RATE
if x4(i) < h_bar
mu_bar = mu*x4(i);
else
mu_bar = mu*h_bar + mu_hat*(x4(i) - h_bar);
end
m11 = -beta*x1(i)*x2(i)*(1-u1(i)) - psi*u3(i)*x1(i);
m12 = beta*x1(i)*x2(i)*(1-u1(i)) - gamma_i*x2(i) - xi_i*x2(i) - nu*x2(i)*(u2(i));
m13 = nu*x2(i)*(u2(i))-gamma_d*x3(i)-xi_d*x3(i);
m14 = xi_i*x2(i)+xi_d*x3(i)-gamma_a*x4(i)-mu_bar;
m15 = gamma_i*x2(i) + gamma_d*x3(i) + gamma_a*x4(i);
m16 = mu_bar;
m17 = psi*u3(i)*x1(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*x2(j)*(1-u1(j)) + L7(j)*psi*u3(j);
n12 = (L1(j)-L2(j))*beta*x1(j)*(1-u1(j)) + L2(j)*(gamma_i + xi_i) + L2(j)*(nu*(u2(j))) - L3(j)*(nu*(u2(j))) - L4(j)*xi_i;
n13 = L3(j)*(gamma_d + xi_d) - L4(j)*xi_d;
n14 = L4(j)*(gamma_a + mu_bar) - L5(j)*mu_bar - L6(j)*mu_bar - c1*x4(j);
n15 = 0;
n16 = 0;
n17 = 0;
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(0.8, max(0, (L2 - L1).*beta.*x1.*x2./(b1)));
u1 = 0.01.*U1 + 0.99.*oldu1;
U2 = min(1.0, max(0, (((L2 - L3).*nu.*x2)./(b2))));
u2 = 0.01.*U2 + 0.99.*oldu2;
U3 = min(1.0, max(0, (((L1 - L7).*psi.*x1)./(b3))));
u3 = 0.01.*U3 + 0.99.*oldu3;
% COST FUNCTION
J = c1./2*sum(x4.^2)*h + b1./2*sum(u1.^2)*h + b2./2*sum(u2.^2)*h + b3./2*sum(u3.^2)*h + c2*max(x6);
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]);
end
disp(['number of loops: ' num2str(loopcnt)]);
number of loops: 2
disp(['Cost function: ' num2str(J)]);
Cost function: 0
disp(['Portion deceased: ' num2str(max(x6))]);
Portion deceased: 0.016393
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;
% IMMUNITY REACHED
imm = x5 + x7.*100; % Percentage immune
x4_per = x4*100; % percentage threatened
x6_per = x6*100; % percentage deceased
MAP = colorcube (M)
MAP = 999×3
0.1250 0.1250 0 0.1250 0.2500 0 0.1250 0.3750 0 0.1250 0.5000 0 0.1250 0.6250 0 0.1250 0.7500 0 0.1250 0.8750 0 0.1250 1.0000 0 0.2500 0.1250 0 0.2500 0.2500 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
figure
colormap %spring, rainbow
ans = 256×3
0.2422 0.1504 0.6603 0.2444 0.1534 0.6728 0.2464 0.1569 0.6847 0.2484 0.1607 0.6961 0.2503 0.1648 0.7071 0.2522 0.1689 0.7179 0.2540 0.1732 0.7286 0.2558 0.1773 0.7393 0.2576 0.1814 0.7501 0.2594 0.1854 0.7610
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
colorbar
for i = 1:19
plot(y(1,:), y(i,:)), grid on, xlabel('t'), ylabel('x_{1}') % plotting x1
hold on
colorbar
end

Categories

Find more on Get Started with Optimization Toolbox in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!