how to draw graphs for optimal control

조회 수: 3(최근 30일)
mallela ankamma rao
mallela ankamma rao 2022년 7월 25일
댓글: mallela ankamma rao 2022년 7월 26일
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

답변(1개)

Sam Chak
Sam Chak 2022년 7월 25일
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
  댓글 수: 1
mallela ankamma rao
mallela ankamma rao 2022년 7월 26일
Thanks for your replying sir. my doubt how did they write the following graphs.
Figure Optimization control input u1 - step size 1,000
only u1(t)
only u2(t)
only u3(t)
here I = x2
i requesting you please help me in these graphs

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by