how to draw graphs for optimal control

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

Answer 1

[Sam Chak]

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

[mallela ankamma rao]

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

[Oladipupo Johnson]

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 [e-mail]. Thank you Sir

[Oladipupo Johnson]

Here Great People. Thank you

Answer 2

[Pavl M.]

  
  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