At one time, if I had a .mat file with a user-defined class variable, but lacked the classdef file, loading the .mat file would succeed, but convert the class property data to struct form. Now, in more recent Matlab, it fails to do that. Is there a supported way to load a user-defined object in a stripped down form, like a struct, when the class definition is absent?At one time, if I had a .mat file with a user-defined class variable, but lacked the classdef file, loading the .mat file would succeed, but convert the class property data to struct form. Now, in more recent Matlab, it fails to do that. Is there a supported way to load a user-defined object in a stripped down form, like a struct, when the class definition is absent? At one time, if I had a .mat file with a user-defined class variable, but lacked the classdef file, loading the .mat file would succeed, but convert the class property data to struct form. Now, in more recent Matlab, it fails to do that. Is there a supported way to load a user-defined object in a stripped down form, like a struct, when the class definition is absent? oop, mat file, load, struct MATLAB Answers — New Questions

​

I have been trying to solve a *Non linear Eigenvalue problem* using *fsolve and Newton’s iteration method* and have not been successful.
The matrix which I am looking to solve:

A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2); 0 0 cos(w) -sin(w)]; (Found in a paper, using it as a practice case)

*My Newton method code:*

%Define intial values and tolerances for the variable
w0=0.1;
tol=2;
maxiter=1000;
w=w0;
wold=w0;
lambda=0.1;

%Start Iteration
for i=1:maxiter
%Define A and B
A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2);
0 0 cos(w) -sin(w)];

B=[-2 0 0 0;-0.5*cos(w/2) 0.5*sin(w/2) 0.5*cos(w/2) -0.5*sin(w/2); sin(w/2) cos(w/2) -0.5*sin(w/2) -0.5*cos(w/2);
0 0 sin(w) cos(w)];
C=inv(B);

%Find Eigen value for the intermediate step
beta=eig(C*A);
epsilon=min(abs(beta));

%Update the variable
w=w0+epsilon;
err=abs(epsilon);
wold=w;
if(err<tol)
break;
end
end

*Fsolve code*

function fval=fun4evp(w)
A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2);
0 0 cos(w) -sin(w)];
fval=det(A);
end

wsol=fsolve(@(w)fun4evp,0.1);

ThanksI have been trying to solve a *Non linear Eigenvalue problem* using *fsolve and Newton’s iteration method* and have not been successful.
The matrix which I am looking to solve:

A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2); 0 0 cos(w) -sin(w)]; (Found in a paper, using it as a practice case)

*My Newton method code:*

%Define intial values and tolerances for the variable
w0=0.1;
tol=2;
maxiter=1000;
w=w0;
wold=w0;
lambda=0.1;

%Start Iteration
for i=1:maxiter
%Define A and B
A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2);
0 0 cos(w) -sin(w)];

B=[-2 0 0 0;-0.5*cos(w/2) 0.5*sin(w/2) 0.5*cos(w/2) -0.5*sin(w/2); sin(w/2) cos(w/2) -0.5*sin(w/2) -0.5*cos(w/2);
0 0 sin(w) cos(w)];
C=inv(B);

%Find Eigen value for the intermediate step
beta=eig(C*A);
epsilon=min(abs(beta));

%Update the variable
w=w0+epsilon;
err=abs(epsilon);
wold=w;
if(err<tol)
break;
end
end

*Fsolve code*

function fval=fun4evp(w)
A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2);
0 0 cos(w) -sin(w)];
fval=det(A);
end

wsol=fsolve(@(w)fun4evp,0.1);

Thanks I have been trying to solve a *Non linear Eigenvalue problem* using *fsolve and Newton’s iteration method* and have not been successful.
The matrix which I am looking to solve:

A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2); 0 0 cos(w) -sin(w)]; (Found in a paper, using it as a practice case)

*My Newton method code:*

%Define intial values and tolerances for the variable
w0=0.1;
tol=2;
maxiter=1000;
w=w0;
wold=w0;
lambda=0.1;

%Start Iteration
for i=1:maxiter
%Define A and B
A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2);
0 0 cos(w) -sin(w)];

B=[-2 0 0 0;-0.5*cos(w/2) 0.5*sin(w/2) 0.5*cos(w/2) -0.5*sin(w/2); sin(w/2) cos(w/2) -0.5*sin(w/2) -0.5*cos(w/2);
0 0 sin(w) cos(w)];
C=inv(B);

%Find Eigen value for the intermediate step
beta=eig(C*A);
epsilon=min(abs(beta));

%Update the variable
w=w0+epsilon;
err=abs(epsilon);
wold=w;
if(err<tol)
break;
end
end

*Fsolve code*

function fval=fun4evp(w)
A=[2*w -300 0 0;sin(w/2) cos(w/2) -sin(w/2) cos(w/2); 2*cos(w/2) -2*sin(w/2) -cos(w/2) sin(w/2);
0 0 cos(w) -sin(w)];
fval=det(A);
end

wsol=fsolve(@(w)fun4evp,0.1);

Thanks fsolve, eigenvalue, nonlinear, nonlinear eigenvalue problem MATLAB Answers — New Questions

​

<</matlabcentral/answers/uploaded_files/1822926/Nnt1.jpg>>

<</matlabcentral/answers/uploaded_files/1822927/Nnt2.jpg>>

I have a simulink circuit which has two PI controllers. I wanted to replace them with neural network blocks. I trained neural network blocks from data for normal conditions of the circuit. The neural network blocks give good results. But when I changed the normal condition of the circuit (changed input voltage) the output was unexpected. I need to know can the Ann blocks work for changed conditions of circuit? How can I make it able to work at variable conditions? The training results are given above.<</matlabcentral/answers/uploaded_files/1822926/Nnt1.jpg>>

<</matlabcentral/answers/uploaded_files/1822927/Nnt2.jpg>>

I have a simulink circuit which has two PI controllers. I wanted to replace them with neural network blocks. I trained neural network blocks from data for normal conditions of the circuit. The neural network blocks give good results. But when I changed the normal condition of the circuit (changed input voltage) the output was unexpected. I need to know can the Ann blocks work for changed conditions of circuit? How can I make it able to work at variable conditions? The training results are given above. <</matlabcentral/answers/uploaded_files/1822926/Nnt1.jpg>>

<</matlabcentral/answers/uploaded_files/1822927/Nnt2.jpg>>

I have a simulink circuit which has two PI controllers. I wanted to replace them with neural network blocks. I trained neural network blocks from data for normal conditions of the circuit. The neural network blocks give good results. But when I changed the normal condition of the circuit (changed input voltage) the output was unexpected. I need to know can the Ann blocks work for changed conditions of circuit? How can I make it able to work at variable conditions? The training results are given above. variable conditions of circuit MATLAB Answers — New Questions

​

Hi there!
For 3D surface plots, do you like the command, ‘shading interp’? When I use it, I notice that the color variation on the surface is quite smooth, with the faint lines that typically show the cross sections of the surface not visible anymore. However, I’m not sure that this is a good thing.
Thanks in advance,Hi there!
For 3D surface plots, do you like the command, ‘shading interp’? When I use it, I notice that the color variation on the surface is quite smooth, with the faint lines that typically show the cross sections of the surface not visible anymore. However, I’m not sure that this is a good thing.
Thanks in advance, Hi there!
For 3D surface plots, do you like the command, ‘shading interp’? When I use it, I notice that the color variation on the surface is quite smooth, with the faint lines that typically show the cross sections of the surface not visible anymore. However, I’m not sure that this is a good thing.
Thanks in advance, matlab surface plot shading MATLAB Answers — New Questions

​

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

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 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

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 code plotting MATLAB Answers — New Questions

​

There is no "From File" block in the Simulink Real-Time (SLRT) library in R2020b. When I ran Upgrade Advisor on my model, it just recommended to "Delete the block".
Is there an alternative workflow to read data from the target into my Speedgoat real-time application in R2020b onwards? Or other ways to specify different input data from MATLAB without the need to rebuild the application?There is no "From File" block in the Simulink Real-Time (SLRT) library in R2020b. When I ran Upgrade Advisor on my model, it just recommended to "Delete the block".
Is there an alternative workflow to read data from the target into my Speedgoat real-time application in R2020b onwards? Or other ways to specify different input data from MATLAB without the need to rebuild the application? There is no "From File" block in the Simulink Real-Time (SLRT) library in R2020b. When I ran Upgrade Advisor on my model, it just recommended to "Delete the block".
Is there an alternative workflow to read data from the target into my Speedgoat real-time application in R2020b onwards? Or other ways to specify different input data from MATLAB without the need to rebuild the application? slrt, speedgoat, input, read, data, file, import MATLAB Answers — New Questions

​

I have two lines that im plotting for a vehicles location. One is the intended path I wanted the vehicle to travel, its a perfectly straight line, the other is the path the vehicle actually took. Its close but not perfect. My end goal is to find how much the vehicle actual path deviated from the intended path as a proxy for navigational accuracy. My idea is that I could calculate the area between the two lines and use that as an indicator for how accuratly the vehicle followed its intended path. When I do this calculation I also want to highlight the difference between the two lines in the graph. This is all part of a much larger program intended for use by the people operating these vehicles. Below is how im plotting the position of both the vehicle and the planned path, I also included a jpg of what the path looks like compared to the intened path.
Ive looked into this and the error I keep getting is:
Specify the coordinates as matrices of the same size or as vectors with the same number of elements.
The problem is, there is only 71 points for the track path and 508486 points for vehicle position, I tried interpolating the track path to make the matrices the same size but I get NaN values for everything after the 71 points so im not sure why thats not working.
I tried patch and fill and got the same error with both. Below is my code:
P = uigetdir(‘C:’);
S1 = dir(fullfile(P,’AHR2.csv’));
S3 = dir(fullfile(P,’CMD.csv’));

F = fullfile(S3.folder,S3.name);
M = readmatrix(F);
F1 = fullfile(S1(k).folder,S1(k).name);
M1 = readmatrix(F1);

TrackLat1 = M(:,10);
TrackLong1 = M(:,11);

Lat2 = M1(:,7);
Long2 = M1(:,8);

TrackLatLong = [TrackLat1, TrackLong1];
[MM, ~, ~] = rmoutliers(TrackLatLong, 1);
TrackLat = MM(:,1);
TrackLong = MM(:,2);

LatLong = [Lat2, Long2];
[MM1, ~, ~] = rmoutliers(LatLong, 1);
Lat1 = MM1(:,1);
Long1 = MM1(:,2);

plot(Long1, Lat1)
hold on
plot(TrackLong, TrackLat)
hold on
X = [Long1;Lat1];
Y = [TrackLong; TrackLat];
patch([X fliplr(X)], [Y fliplr(Y)], ‘red’)

title ‘FAT Position’
legend(‘USV Position’, ‘Mission Plan’)

xlim tight
ylim tight
ylabel ‘Latitude’
xlabel ‘Longitude’
attached is a screen shot of the plot plus a zoomed in section because you might not be able to tell. Thanks in advance for any help I might recieve, if any one has any ideas on a simpler way to do this please let me know!I have two lines that im plotting for a vehicles location. One is the intended path I wanted the vehicle to travel, its a perfectly straight line, the other is the path the vehicle actually took. Its close but not perfect. My end goal is to find how much the vehicle actual path deviated from the intended path as a proxy for navigational accuracy. My idea is that I could calculate the area between the two lines and use that as an indicator for how accuratly the vehicle followed its intended path. When I do this calculation I also want to highlight the difference between the two lines in the graph. This is all part of a much larger program intended for use by the people operating these vehicles. Below is how im plotting the position of both the vehicle and the planned path, I also included a jpg of what the path looks like compared to the intened path.
Ive looked into this and the error I keep getting is:
Specify the coordinates as matrices of the same size or as vectors with the same number of elements.
The problem is, there is only 71 points for the track path and 508486 points for vehicle position, I tried interpolating the track path to make the matrices the same size but I get NaN values for everything after the 71 points so im not sure why thats not working.
I tried patch and fill and got the same error with both. Below is my code:
P = uigetdir(‘C:’);
S1 = dir(fullfile(P,’AHR2.csv’));
S3 = dir(fullfile(P,’CMD.csv’));

F = fullfile(S3.folder,S3.name);
M = readmatrix(F);
F1 = fullfile(S1(k).folder,S1(k).name);
M1 = readmatrix(F1);

TrackLat1 = M(:,10);
TrackLong1 = M(:,11);

Lat2 = M1(:,7);
Long2 = M1(:,8);

TrackLatLong = [TrackLat1, TrackLong1];
[MM, ~, ~] = rmoutliers(TrackLatLong, 1);
TrackLat = MM(:,1);
TrackLong = MM(:,2);

LatLong = [Lat2, Long2];
[MM1, ~, ~] = rmoutliers(LatLong, 1);
Lat1 = MM1(:,1);
Long1 = MM1(:,2);

plot(Long1, Lat1)
hold on
plot(TrackLong, TrackLat)
hold on
X = [Long1;Lat1];
Y = [TrackLong; TrackLat];
patch([X fliplr(X)], [Y fliplr(Y)], ‘red’)

title ‘FAT Position’
legend(‘USV Position’, ‘Mission Plan’)

xlim tight
ylim tight
ylabel ‘Latitude’
xlabel ‘Longitude’
attached is a screen shot of the plot plus a zoomed in section because you might not be able to tell. Thanks in advance for any help I might recieve, if any one has any ideas on a simpler way to do this please let me know! I have two lines that im plotting for a vehicles location. One is the intended path I wanted the vehicle to travel, its a perfectly straight line, the other is the path the vehicle actually took. Its close but not perfect. My end goal is to find how much the vehicle actual path deviated from the intended path as a proxy for navigational accuracy. My idea is that I could calculate the area between the two lines and use that as an indicator for how accuratly the vehicle followed its intended path. When I do this calculation I also want to highlight the difference between the two lines in the graph. This is all part of a much larger program intended for use by the people operating these vehicles. Below is how im plotting the position of both the vehicle and the planned path, I also included a jpg of what the path looks like compared to the intened path.
Ive looked into this and the error I keep getting is:
Specify the coordinates as matrices of the same size or as vectors with the same number of elements.
The problem is, there is only 71 points for the track path and 508486 points for vehicle position, I tried interpolating the track path to make the matrices the same size but I get NaN values for everything after the 71 points so im not sure why thats not working.
I tried patch and fill and got the same error with both. Below is my code:
P = uigetdir(‘C:’);
S1 = dir(fullfile(P,’AHR2.csv’));
S3 = dir(fullfile(P,’CMD.csv’));

F = fullfile(S3.folder,S3.name);
M = readmatrix(F);
F1 = fullfile(S1(k).folder,S1(k).name);
M1 = readmatrix(F1);

TrackLat1 = M(:,10);
TrackLong1 = M(:,11);

Lat2 = M1(:,7);
Long2 = M1(:,8);

TrackLatLong = [TrackLat1, TrackLong1];
[MM, ~, ~] = rmoutliers(TrackLatLong, 1);
TrackLat = MM(:,1);
TrackLong = MM(:,2);

LatLong = [Lat2, Long2];
[MM1, ~, ~] = rmoutliers(LatLong, 1);
Lat1 = MM1(:,1);
Long1 = MM1(:,2);

plot(Long1, Lat1)
hold on
plot(TrackLong, TrackLat)
hold on
X = [Long1;Lat1];
Y = [TrackLong; TrackLat];
patch([X fliplr(X)], [Y fliplr(Y)], ‘red’)

title ‘FAT Position’
legend(‘USV Position’, ‘Mission Plan’)

xlim tight
ylim tight
ylabel ‘Latitude’
xlabel ‘Longitude’
attached is a screen shot of the plot plus a zoomed in section because you might not be able to tell. Thanks in advance for any help I might recieve, if any one has any ideas on a simpler way to do this please let me know! patch MATLAB Answers — New Questions

​

Hi, I went through a course that I found through the mathworks website and I am having trouble finding the link back to the section where it listed all the courses where the course material had been uploaded online. Does anyone know where is is located?Hi, I went through a course that I found through the mathworks website and I am having trouble finding the link back to the section where it listed all the courses where the course material had been uploaded online. Does anyone know where is is located? Hi, I went through a course that I found through the mathworks website and I am having trouble finding the link back to the section where it listed all the courses where the course material had been uploaded online. Does anyone know where is is located? course MATLAB Answers — New Questions

​

Hello,

I am trying to model a linear actuator in simscape matlab but have been quite unsuccessful. The way I want the model to work is to apply a load at the end of the lead screw and set a speed at which the motor runs and from that information measure a current that is drawn by the motor.

My main issue is that I am unable to make the BLDC (https://uk.mathworks.com/help/sps/ref/bldc.html) drive my gears. I am using a BLDC current controller with PWM generation to power the gates of a converter that should then power the BLDC. I have attached the part of the circuit I think is the problem but honestly I don’t really know if that is the problem. What I think the problem is that the converter is not getting the correct pulse to drive which is most likely cause by the way i have set up my BLDC controller. I will also attach the whole circuit.

Many thanks in advance.Hello,

I am trying to model a linear actuator in simscape matlab but have been quite unsuccessful. The way I want the model to work is to apply a load at the end of the lead screw and set a speed at which the motor runs and from that information measure a current that is drawn by the motor.

My main issue is that I am unable to make the BLDC (https://uk.mathworks.com/help/sps/ref/bldc.html) drive my gears. I am using a BLDC current controller with PWM generation to power the gates of a converter that should then power the BLDC. I have attached the part of the circuit I think is the problem but honestly I don’t really know if that is the problem. What I think the problem is that the converter is not getting the correct pulse to drive which is most likely cause by the way i have set up my BLDC controller. I will also attach the whole circuit.

Many thanks in advance. Hello,

I am trying to model a linear actuator in simscape matlab but have been quite unsuccessful. The way I want the model to work is to apply a load at the end of the lead screw and set a speed at which the motor runs and from that information measure a current that is drawn by the motor.

My main issue is that I am unable to make the BLDC (https://uk.mathworks.com/help/sps/ref/bldc.html) drive my gears. I am using a BLDC current controller with PWM generation to power the gates of a converter that should then power the BLDC. I have attached the part of the circuit I think is the problem but honestly I don’t really know if that is the problem. What I think the problem is that the converter is not getting the correct pulse to drive which is most likely cause by the way i have set up my BLDC controller. I will also attach the whole circuit.

Many thanks in advance. simscape, simulink, model, electric_motor_control MATLAB Answers — New Questions

​

Hi,
We are trying to simulate a circuit and we face attached error. Is there any particular reason why this comes and kindly let us know the solution for the same as well.Hi,
We are trying to simulate a circuit and we face attached error. Is there any particular reason why this comes and kindly let us know the solution for the same as well. Hi,
We are trying to simulate a circuit and we face attached error. Is there any particular reason why this comes and kindly let us know the solution for the same as well. simulink, error MATLAB Answers — New Questions

​