how to draw graphs for optimal control
19 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
0 Kommentare
Antworten (1)
Sam Chak
am 25 Jul. 2022
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)]);
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
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 Kommentare
Oladipupo Johnson
am 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
Siehe auch
Kategorien
Mehr zu Model Predictive Control Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!