I have a code for optimal control of a deterministic SEIR model. I am unable to code the stochastic version of the model.

test = -1; % Test variable; as long as variable is negative …the while loops keeps repeating
t0 = 0; % Initial time
tf = 100; % 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
alpha = 1.5;
mu = 0.05;
b1 = 0.05;
b2 = 0.03;
lambda = 1.2;
epsilon2 = 0.55;
epsilon1 = 0.01;
gamma = 0.2;
% WEIGHT FACTORS
k1 = 0.05; % Weight on threatened population
k2 = 0.05; % Weigth on deceased population
l1 = 0.05;
l2 = 0.05;
m1 = 0.025;
m2 = 0.025;
m3 = 0.025;
m4 = 0.025;
% INITIAL CONDITIONS MODEL
S = zeros(1, M+1); % Susceptible
E = zeros(1, M+1); % Exposed
I = zeros(1, M+1); % Infected
R = zeros(1, M+1); % Recovered
S(1) = 99999;
E(1) = 100;
I(1) = 100;
R(1) = 50;
% 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
% INITIAL CONDITIONS ADOINT SYSTEM
lambda1 = zeros(1, M+1);
lambda2 = zeros(1, M+1);
lambda3 = zeros(1, M+1);
lambda4 = zeros(1, M+1);
lambda1(M+1) = 0;
lambda2(M+1) = 0;
lambda3(M+1) = 0;
lambda4(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldu2 = u2;
oldS = S;
oldE = E;
oldI = I;
oldR = R;
oldlambda1 = lambda1;
oldlambda2 = lambda2;
oldlambda3 = lambda3;
oldlambda4 = lambda4;

% SYSTEM DYNAMICCS
for i = 1:M

if i > round(tau1 / h)
E_tau1 = E(i-round(tau1/h));
else
E_tau1 = E(1);
end

if i > round(tau2 / h)
I_tau2 = I(i-round(tau2/h));
else
I_tau2 = I(1);
end
m11 = alpha – b1 * (1-u1(i)) * S(i) * E(i) / (S(i)+E(i)+I(i)+R(i)) – b2 * (1-u1(i)) * S(i) * I(i) / (S(i)+E(i)+I(i)+R(i))- mu * S(i);
m12 = b1 * (1-u1(i)) * S(i) * E(i) / (S(i)+E(i)+I(i)+R(i)) + b2 * (1-u1(i)) * S(i) * E(i) / (S(i)+E(i)+I(i)+R(i)) – lambda * E_tau1 – (epsilon1 + mu) * E(i);
m13 = lambda * E_tau1 – epsilon2 * (1+u2(i)) * I_tau2 – (gamma + mu) * I(i);
m14 = epsilon1 * (1+u2(i)) * E(i) + epsilon2 * I_tau2 – mu * R(i);

m21 = alpha – b1 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (E(i)+h2*m12) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) – b2 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (I(i)+h2*m13) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) – mu * (S(i)+h2*m11);
m22 = b1 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (E(i)+h2*m12) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) + b2 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (I(i)+h2*m13) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) – lambda * E_tau1 – (epsilon1 + mu) * (E(i)+h2*m12);
m23 = lambda * E_tau1 – epsilon2 * ( 1 + 0.5*(u2(i)+u2(i+1))) * I_tau2 – (gamma + mu) * (I(i)+h2*m13);
m24 = epsilon1 * ( 1 + 0.5*(u2(i)+u2(i+1))) * (E(i)+h2*m12) + epsilon2 * I_tau2 – mu * (R(i)+h2*m14);
m31 = alpha – b1 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (E(i)+h2*m22) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) – b2 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (I(i)+h2*m23) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) – mu * (S(i)+h2*m21);
m32 = b1 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (E(i)+h2*m22) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) + b2 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (I(i)+h2*m23) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) – lambda * E_tau1 – (epsilon1 + mu) * (E(i)+h2*m22);
m33 = lambda * E_tau1 – epsilon2 * ( 1 + 0.5*(u2(i)+u2(i+1))) * I_tau2 – (gamma + mu) * (I(i)+h2*m23);
m34 = epsilon1 * ( 1 + 0.5*(u2(i)+u2(i+1))) * (E(i)+h2*m22) + epsilon2 * I_tau2 – mu * (R(i)+h2*m24);
m41 = alpha – b1 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (E(i)+h2*m32) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) – b2 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (I(i)+h2*m33) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) – mu * (S(i)+h2*m31);
m42 = b1 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (E(i)+h2*m32) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) + b2 * ( 1 – 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (I(i)+h2*m33) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) – lambda * E_tau1 – (epsilon1 + mu) * (E(i)+h2*m32);
m43 = lambda * E_tau1 – epsilon2 * ( 1 + 0.5*(u2(i)+u2(i+1))) * I_tau2 – (gamma + mu) * (I(i)+h2*m33);
m44 = epsilon1 * ( 1 + 0.5*(u2(i)+u2(i+1))) * (E(i)+h2*m32) + epsilon2 * I_tau2 – mu * (R(i)+h2*m34);

S(i+1) = S(i) + (h/6)*(m11 + 2*m21 + 2*m31 + m41);
E(i+1) = E(i) + (h/6)*(m12 + 2*m22 + 2*m32 + m42);
I(i+1) = I(i) + (h/6)*(m13 + 2*m23 + 2*m33 + m43);
R(i+1) = R(i) + (h/6)*(m14 + 2*m24 + 2*m34 + m44);
end

% ADJOINT SYSTEM
for i = 1:M % From initial to final value
j = M + 2 – i; % From final value to initial value
n11 = -m1 * S(j) + lambda1(j) * b1 * (1-u1(j)) * E(j) / (S(j)+E(j)+I(j)+R(j)) + lambda1(j) * b2 * (1-u1(j)) * I(j) / (S(j)+E(j)+I(j)+R(j)) + mu * lambda1(j) – b1 * lambda2(j) * (1-u1(j)) * E(j) / (S(j)+E(j)+I(j)+R(j)) – b2 * lambda2(j) * (1-u1(j)) * I(j) / (S(j)+E(j)+I(j)+R(j));
n12 = -k1 – m2 * E(j) + b1 * lambda1(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) – b1 * lambda2(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) + lambda * lambda2(j) + epsilon1 * lambda2(j) + mu * lambda2(j) – lambda * lambda3(j) – epsilon1 * lambda4(j) * (1+u2(j));
n13 = -k2 – m3 * I(j) + b2 * lambda1(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) – b2 * lambda2(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) + epsilon2 * lambda3(j) * (1+u2(j)) + ( gamma + mu) – epsilon2 * lambda4(j);
n14 = -m4 * R(j) + mu * lambda4(j);

n21 = -m1 * 0.5*(S(j) + S(j-1)) + (lambda1(j)-h2*n11) * b1 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + (lambda1(j)-h2*n11) * b2 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + mu * (lambda1(j)-h2*n11) – b1 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1))) – b2 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1)));
n22 = -k1 – m2 * 0.5*(E(j) + E(j-1)) + b1 * (lambda1(j)-h2*n11) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) – b1 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + lambda * (lambda2(j)-h2*n12) + epsilon1 * (lambda2(j)-h2*n12) + mu * (lambda2(j)-h2*n12) – lambda * (lambda3(j)-h2*n13) – epsilon1 * (lambda4(j)-h2*n14) * (1+(0.5*(u2(j) + u2(j-1))));
n23 = -k2 – m3 * 0.5*(I(j) + I(j-1)) + b2 * (lambda1(j)-h2*n11) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) – b2 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + epsilon2 * (lambda3(j)-h2*n13) * (1+0.5*(u2(j)+u2(j-1))) + ( gamma + mu) – epsilon2 * (lambda4(j)-h2*n14);
n24 = -m4 * 0.5*(R(j) + R(j-1)) + mu * (lambda4(j)-h2*n14);
n31 = -m1 * 0.5*(S(j) + S(j-1)) + (lambda1(j)-h2*n21) * b1 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + (lambda1(j)-h2*n21) * b2 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + mu * (lambda1(j)-h2*n21) – b1 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1))) – b2 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1)));
n32 = -k1 – m2 * 0.5*(E(j) + E(j-1)) + b1 * (lambda1(j)-h2*n21) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) – b1 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + lambda * (lambda2(j)-h2*n22) + epsilon1 * (lambda2(j)-h2*n22) + mu * (lambda2(j)-h2*n22) – lambda * (lambda3(j)-h2*n23) – epsilon1 * (lambda4(j)-h2*n24) * (1+(0.5*(u2(j) + u2(j-1))));
n33 = -k2 – m3 * 0.5*(I(j) + I(j-1)) + b2 * (lambda1(j)-h2*n21) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) – b2 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + epsilon2 * (lambda3(j)-h2*n23) * (1+0.5*(u2(j)+u2(j-1))) + ( gamma + mu) – epsilon2 * (lambda4(j)-h2*n24);
n34 = -m4 * 0.5*(R(j) + R(j-1)) + mu * (lambda4(j)-h2*n24);
n41 = -m1 * 0.5*(S(j) + S(j-1)) + (lambda1(j)-h2*n31) * b1 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + (lambda1(j)-h2*n31) * b2 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + mu * (lambda1(j)-h2*n31) – b1 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1))) – b2 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1)));
n42 = -k1 – m2 * 0.5*(E(j) + E(j-1)) + b1 * (lambda1(j)-h2*n31) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) – b1 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + lambda * (lambda2(j)-h2*n32) + epsilon1 * (lambda2(j)-h2*n32) + mu * (lambda2(j)-h2*n32) – lambda * (lambda3(j)-h2*n33) – epsilon1 * (lambda4(j)-h2*n34) * (1+(0.5*(u2(j) + u2(j-1))));
n43 = -k2 – m3 * 0.5*(I(j) + I(j-1)) + b2 * (lambda1(j)-h2*n31) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) – b2 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + epsilon2 * (lambda3(j)-h2*n33) * (1+0.5*(u2(j)+u2(j-1))) + ( gamma + mu) – epsilon2 * (lambda4(j)-h2*n34);
n44 = -m4 * 0.5*(R(j) + R(j-1)) + mu * (lambda4(j)-h2*n34);

lambda1(j-1) = lambda1(j) – (h/6)*(n11 + 2*n21 + 2*n31 + n41);
lambda2(j-1) = lambda2(j) – (h/6)*(n12 + 2*n22 + 2*n32 + n42);
lambda3(j-1) = lambda3(j) – (h/6)*(n13 + 2*n23 + 2*n33 + n43);
lambda4(j-1) = lambda4(j) – (h/6)*(n14 + 2*n24 + 2*n34 + n44);
end

% OPTIMALITY CONDITIONS
U1 = max(0, min(1.0, ((lambda2 – lambda1).*b1.*S.*E + (lambda2 – lambda1).*b2.*S.*I ) ./ (l1 .* (S+E+I+R))));
u1 = 0.01.*U1 + 0.99.*oldu1;
U2 = max(0, min(1.0, ((lambda3 .* epsilon2 .* I – lambda4 .* epsilon1 .* E) ./(l2))));
u2 = 0.01.*U2 + 0.99.*oldu2;

% U1 = max(0, min(1, (L2 – L1).*b1.*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 = k1*sum(E)*h + k2*sum(I)*h + l1*sum(u1.^2)*h + l2*sum(u2.^2)*h + m1./2*sum(S.^2)*h + m2./2*sum(E.^2)*h+ m3./2*sum(I.^2)*h+ m4/2*sum(R.^2)*h;
Cost1 = k1.*cumsum(E)*h; % Total cost of threatened population
Cost2 = k2.*cumsum(I)*h;
Cost3 = l1.*cumsum(u1.^2)*h;
Cost4 = l2.*cumsum(u2.^2)*h;
Cost5 = m1./2.*cumsum(S.^2)*h;
Cost6 = m2./2.*cumsum(E.^2)*h; % Total cost of control input u1
Cost7 = m3./2.*cumsum(I.^2)*h; % Total cost of control input u2
Cost8 = m4./2.*cumsum(R.^2)*h; % Total cost of control input u3
J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5 + Cost6 + Cost7+ Cost8; % Cost at each time for …plotting graphs
% % 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(S)) – sum(abs(oldS – S));
temp4 = Delta*sum(abs(E)) – sum(abs(oldE – E));
temp5 = Delta*sum(abs(I)) – sum(abs(oldI – I));
temp6 = Delta*sum(abs(R)) – sum(abs(oldR – R));
temp7 = Delta*sum(abs(lambda1)) – sum(abs(oldlambda1 – lambda1));
temp8 = Delta*sum(abs(lambda2)) – sum(abs(oldlambda2 – lambda2));
temp9 = Delta*sum(abs(lambda3)) – sum(abs(oldlambda3 – lambda3));
temp10 = Delta*sum(abs(lambda4)) – sum(abs(oldlambda4 – lambda4));
test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10]);
end
disp([‘number of loops: ‘ num2str(loopcnt)]);
disp([‘Cost function: ‘ num2str(J)]);
disp([‘Portion deceased: ‘ num2str(max(x6))]);
y(1,:) = t;
y(2,:) = S;
y(3,:) = E;
y(4,:) = I;
y(5,:) = R;
y(9,:) = lambda1;
y(10,:) = lambda2;
y(11,:) = lambda3;
y(12,:) = lambda4;
y(16,:) = u1;
y(17,:) = u2;
y(19,:) = J;
% IMMUNITY REACHED
% imm = x5 + x7.*100; % Percentage immune
% R_per = R*100; % percentage threatened
% x6_per = x6*100; % percentage deceased
plot(y(1,:), y(2,:)), grid on, xlabel(‘t’), ylabel(‘S’) % plotting S
plot(y(1,:), y(3,:)), grid on, xlabel(‘t’), ylabel(‘E’) % plotting E
plot(y(1,:), y(4,:)), grid on, xlabel(‘t’), ylabel(‘I’) % plotting I
plot(y(1,:), y(5,:)), grid on, xlabel(‘t’), ylabel(‘R’) % plotting R
I want the deterministic version of the codetest = -1; % Test variable; as long as variable is negative …the while loops keeps repeating
t0 = 0; % Initial time
tf = 100; % 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
alpha = 1.5;
mu = 0.05;
b1 = 0.05;
b2 = 0.03;
lambda = 1.2;
epsilon2 = 0.55;
epsilon1 = 0.01;
gamma = 0.2;
% WEIGHT FACTORS
k1 = 0.05; % Weight on threatened population
k2 = 0.05; % Weigth on deceased population
l1 = 0.05;
l2 = 0.05;
m1 = 0.025;
m2 = 0.025;
m3 = 0.025;
m4 = 0.025;
% INITIAL CONDITIONS MODEL
S = zeros(1, M+1); % Susceptible
E = zeros(1, M+1); % Exposed
I = zeros(1, M+1); % Infected
R = zeros(1, M+1); % Recovered
S(1) = 99999;
E(1) = 100;
I(1) = 100;
R(1) = 50;
% 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
% INITIAL CONDITIONS ADOINT SYSTEM
lambda1 = zeros(1, M+1);
lambda2 = zeros(1, M+1);
lambda3 = zeros(1, M+1);
lambda4 = zeros(1, M+1);
lambda1(M+1) = 0;
lambda2(M+1) = 0;
lambda3(M+1) = 0;
lambda4(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldu2 = u2;
oldS = S;
oldE = E;
oldI = I;
oldR = R;
oldlambda1 = lambda1;
oldlambda2 = lambda2;
oldlambda3 = lambda3;
oldlambda4 = lambda4;

% SYSTEM DYNAMICCS
for i = 1:M

if i > round(tau1 / h)
E_tau1 = E(i-round(tau1/h));
else
E_tau1 = E(1);
end

% COST FUNCTION
J = k1*sum(E)*h + k2*sum(I)*h + l1*sum(u1.^2)*h + l2*sum(u2.^2)*h + m1./2*sum(S.^2)*h + m2./2*sum(E.^2)*h+ m3./2*sum(I.^2)*h+ m4/2*sum(R.^2)*h;
Cost1 = k1.*cumsum(E)*h; % Total cost of threatened population
Cost2 = k2.*cumsum(I)*h;
Cost3 = l1.*cumsum(u1.^2)*h;
Cost4 = l2.*cumsum(u2.^2)*h;
Cost5 = m1./2.*cumsum(S.^2)*h;
Cost6 = m2./2.*cumsum(E.^2)*h; % Total cost of control input u1
Cost7 = m3./2.*cumsum(I.^2)*h; % Total cost of control input u2
Cost8 = m4./2.*cumsum(R.^2)*h; % Total cost of control input u3
J2 = Cost1 + Cost2 + Cost3 + Cost4 + Cost5 + Cost6 + Cost7+ Cost8; % Cost at each time for …plotting graphs
% % 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(S)) – sum(abs(oldS – S));
temp4 = Delta*sum(abs(E)) – sum(abs(oldE – E));
temp5 = Delta*sum(abs(I)) – sum(abs(oldI – I));
temp6 = Delta*sum(abs(R)) – sum(abs(oldR – R));
temp7 = Delta*sum(abs(lambda1)) – sum(abs(oldlambda1 – lambda1));
temp8 = Delta*sum(abs(lambda2)) – sum(abs(oldlambda2 – lambda2));
temp9 = Delta*sum(abs(lambda3)) – sum(abs(oldlambda3 – lambda3));
temp10 = Delta*sum(abs(lambda4)) – sum(abs(oldlambda4 – lambda4));
test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10]);
end
disp([‘number of loops: ‘ num2str(loopcnt)]);
disp([‘Cost function: ‘ num2str(J)]);
disp([‘Portion deceased: ‘ num2str(max(x6))]);
y(1,:) = t;
y(2,:) = S;
y(3,:) = E;
y(4,:) = I;
y(5,:) = R;
y(9,:) = lambda1;
y(10,:) = lambda2;
y(11,:) = lambda3;
y(12,:) = lambda4;
y(16,:) = u1;
y(17,:) = u2;
y(19,:) = J;
% IMMUNITY REACHED
% imm = x5 + x7.*100; % Percentage immune
% R_per = R*100; % percentage threatened
% x6_per = x6*100; % percentage deceased
plot(y(1,:), y(2,:)), grid on, xlabel(‘t’), ylabel(‘S’) % plotting S
plot(y(1,:), y(3,:)), grid on, xlabel(‘t’), ylabel(‘E’) % plotting E
plot(y(1,:), y(4,:)), grid on, xlabel(‘t’), ylabel(‘I’) % plotting I
plot(y(1,:), y(5,:)), grid on, xlabel(‘t’), ylabel(‘R’) % plotting R
I want the deterministic version of the code test = -1; % Test variable; as long as variable is negative …the while loops keeps repeating
t0 = 0; % Initial time
tf = 100; % 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
alpha = 1.5;
mu = 0.05;
b1 = 0.05;
b2 = 0.03;
lambda = 1.2;
epsilon2 = 0.55;
epsilon1 = 0.01;
gamma = 0.2;
% WEIGHT FACTORS
k1 = 0.05; % Weight on threatened population
k2 = 0.05; % Weigth on deceased population
l1 = 0.05;
l2 = 0.05;
m1 = 0.025;
m2 = 0.025;
m3 = 0.025;
m4 = 0.025;
% INITIAL CONDITIONS MODEL
S = zeros(1, M+1); % Susceptible
E = zeros(1, M+1); % Exposed
I = zeros(1, M+1); % Infected
R = zeros(1, M+1); % Recovered
S(1) = 99999;
E(1) = 100;
I(1) = 100;
R(1) = 50;
% 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
% INITIAL CONDITIONS ADOINT SYSTEM
lambda1 = zeros(1, M+1);
lambda2 = zeros(1, M+1);
lambda3 = zeros(1, M+1);
lambda4 = zeros(1, M+1);
lambda1(M+1) = 0;
lambda2(M+1) = 0;
lambda3(M+1) = 0;
lambda4(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
loopcnt = loopcnt + 1;
oldu1 = u1;
oldu2 = u2;
oldS = S;
oldE = E;
oldI = I;
oldR = R;
oldlambda1 = lambda1;
oldlambda2 = lambda2;
oldlambda3 = lambda3;
oldlambda4 = lambda4;

% SYSTEM DYNAMICCS
for i = 1:M

if i > round(tau1 / h)
E_tau1 = E(i-round(tau1/h));
else
E_tau1 = E(1);
end

Cart

Cart