how to draw graphs for optimal control

31 visualizaciones (últimos 30 días)
mallela ankamma rao
mallela ankamma rao el 25 de Jul. de 2022
Respondida: Pavl M. hace alrededor de 8 horas
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

Iniciar sesión para responder a esta pregunta.

Respuestas (2)

Sam Chak
Sam Chak el 25 de Jul. de 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 comentarios
Mostrar 1 comentario más antiguo
Oladipupo Johnson
Oladipupo Johnson el 10 de Abr. de 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
Oladipupo Johnson
Oladipupo Johnson el 10 de Abr. de 2023
Here Great People. Thank you

Iniciar sesión para comentar.

Pavl M.
Pavl M. hace alrededor de 8 horas
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

Categorías

Más información sobre Get Started with Optimization Toolbox en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by