Month: January 2025
How to solve this error ?
Post Content Post Content @ cellular communication systems MATLAB Answers — New Questions
For 3D surface plots, do you like the command, ‘shading interp’?
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
Importing a PDF Bank Statement into MATLAB and splitting transactions correctly
Hi, my bank statement comes as a PDF, but I cannot for the life of me import it into Excel correctly in order to do some analysis. I decided to try MATLAB and see if it imports the data as intended.
The PDF starts with a bunch of text that isn’t useful, such as name and address etc. I want to start the import at the point in the document where the transactions begin.
The columns split into [Trans Data], [Post Date], [Description], [Amount]
Then each line below that has data corresponding to that data point.
The columns are only split by a space, but a data entry may also have spaces inside it, so spaces cannot be used to determine where one cell finishes and the next one starts.
Any ideas?
Bonus if I can export the final result as an Excel file. If not I can keep it in MATLABHi, my bank statement comes as a PDF, but I cannot for the life of me import it into Excel correctly in order to do some analysis. I decided to try MATLAB and see if it imports the data as intended.
The PDF starts with a bunch of text that isn’t useful, such as name and address etc. I want to start the import at the point in the document where the transactions begin.
The columns split into [Trans Data], [Post Date], [Description], [Amount]
Then each line below that has data corresponding to that data point.
The columns are only split by a space, but a data entry may also have spaces inside it, so spaces cannot be used to determine where one cell finishes and the next one starts.
Any ideas?
Bonus if I can export the final result as an Excel file. If not I can keep it in MATLAB Hi, my bank statement comes as a PDF, but I cannot for the life of me import it into Excel correctly in order to do some analysis. I decided to try MATLAB and see if it imports the data as intended.
The PDF starts with a bunch of text that isn’t useful, such as name and address etc. I want to start the import at the point in the document where the transactions begin.
The columns split into [Trans Data], [Post Date], [Description], [Amount]
Then each line below that has data corresponding to that data point.
The columns are only split by a space, but a data entry may also have spaces inside it, so spaces cannot be used to determine where one cell finishes and the next one starts.
Any ideas?
Bonus if I can export the final result as an Excel file. If not I can keep it in MATLAB import pdf data MATLAB Answers — New Questions
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
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
How to log signals in SIL Testing?
I am running a simulation as SIL (Software-In-The-Loop).
I would like to log some signals inside the model I am testing.
How can I do that? I am running a simulation as SIL (Software-In-The-Loop).
I would like to log some signals inside the model I am testing.
How can I do that? I am running a simulation as SIL (Software-In-The-Loop).
I would like to log some signals inside the model I am testing.
How can I do that? sil, pil, signal, logging MATLAB Answers — New Questions
How can I provide input data to my Speedgoat application using Simulink Real-Time R2020b onwards?
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
How to plot a phase line plot for an one-dimensional equation with MATLAB?
How do you plot a phase line plot for an one-dimensional equation with MATLAB?
e.g. f(t,x) = x*((m*(1-x)/(1+a-x))-1), m = 0.99, a = 0.0101How do you plot a phase line plot for an one-dimensional equation with MATLAB?
e.g. f(t,x) = x*((m*(1-x)/(1+a-x))-1), m = 0.99, a = 0.0101 How do you plot a phase line plot for an one-dimensional equation with MATLAB?
e.g. f(t,x) = x*((m*(1-x)/(1+a-x))-1), m = 0.99, a = 0.0101 phase, line, plot, one-dimensional MATLAB Answers — New Questions
How to find the area between two lines of different size matrices and fill that area?
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
Convert .txt to .csv with prespecified format
I have a (.txt) with this format
2023 JAN 1 00 31 34.1 38.5625 23.6833 14 2.0
2023 JAN 1 00 38 24.7 38.3304 20.4172 19 1.0
2023 JAN 1 00 47 15.0 38.3940 21.9118 7 0.9
and i want to make a file with the following format:
DATETIME;LAT;LONG;DEPTH;MAG
01-01-2023T00:31:34.1;38.5625;23.6833;14;2.0
01-01-2023T00:38:24.7;38.3304;20.4172;19;1.0
01-01-2023T00:47:15.0;38.3940;21.9118;7;0.9
Any help?I have a (.txt) with this format
2023 JAN 1 00 31 34.1 38.5625 23.6833 14 2.0
2023 JAN 1 00 38 24.7 38.3304 20.4172 19 1.0
2023 JAN 1 00 47 15.0 38.3940 21.9118 7 0.9
and i want to make a file with the following format:
DATETIME;LAT;LONG;DEPTH;MAG
01-01-2023T00:31:34.1;38.5625;23.6833;14;2.0
01-01-2023T00:38:24.7;38.3304;20.4172;19;1.0
01-01-2023T00:47:15.0;38.3940;21.9118;7;0.9
Any help? I have a (.txt) with this format
2023 JAN 1 00 31 34.1 38.5625 23.6833 14 2.0
2023 JAN 1 00 38 24.7 38.3304 20.4172 19 1.0
2023 JAN 1 00 47 15.0 38.3940 21.9118 7 0.9
and i want to make a file with the following format:
DATETIME;LAT;LONG;DEPTH;MAG
01-01-2023T00:31:34.1;38.5625;23.6833;14;2.0
01-01-2023T00:38:24.7;38.3304;20.4172;19;1.0
01-01-2023T00:47:15.0;38.3940;21.9118;7;0.9
Any help? convert, file MATLAB Answers — New Questions
Matlab course
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
“Error using parpool no space left on device” — even when there is enough space.
When I attempt to run my code, it immidiately fails with the warning:
“`
Error using parpool (line 133)
No space left on device
Error in testing (line 39)
parpool(2)
“`
Even though I have more than enough space available in my tmp folder and storage (>1TB).
What might be the cause of this error?When I attempt to run my code, it immidiately fails with the warning:
“`
Error using parpool (line 133)
No space left on device
Error in testing (line 39)
parpool(2)
“`
Even though I have more than enough space available in my tmp folder and storage (>1TB).
What might be the cause of this error? When I attempt to run my code, it immidiately fails with the warning:
“`
Error using parpool (line 133)
No space left on device
Error in testing (line 39)
parpool(2)
“`
Even though I have more than enough space available in my tmp folder and storage (>1TB).
What might be the cause of this error? parallel computing toolbox, parallel computing MATLAB Answers — New Questions
Linear actuator current simulation model
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
How to prevent wraparound using geoplot
I am trying to accomplish something similar to the below post, namely prevent my data from wrapping at the meridians. So far, no luck:
https://www.mathworks.com/matlabcentral/answers/2172159-how-to-wrap-axes-using-geoplot
If I set the lonLim to [-360 360] and call geolimits(latlim,lonlim), the longitude limits do NOT change to [-360 360] but to [-180.3914 181.1743]. Why is this? Is there a bult-in max limit in geolimits? It’s possible that the function is confused by the [-360 360] values, which are effectively equivalent.
But back to the wraparound issue. I am plotting NetCDF-type satellite observation data.
figure
gx = geoaxes;
geobasemap(gx,’usgsimagery’);
hold on
latLim = [-72 85];
lonLim = [-180 180];
geolimits(gx, latLim, lonLim);
lat = ncread(filepath, ‘Latitude’); % read NetCDF-type data
lon = ncread(filepath, ‘Longitude’);
lat1 = double(lat(:,:,1)) ; % read just the 1st of 6 layers
lon1 = double(lon(:,:,1));
lon1 = wrapTo180(lon1); % limit to +/- 180 degrees
k = boundary(lat1(:),lon1(:)); % get just the edge data
pgon = geopolyshape(lat1(k),lon1(k));
geoplot(gx,pgon);
I am able to plot the satellite observation patches on the Mercator projection geoaxes display. However, when a patch crosses the International Date Line, it does not straddle the line but whangs all the way to the right, to the Prime Meridian.
If I use geoplot instead of geopolyshape I get the same effect, only with vectors instead of shapes.When the longitude data reaches -180, the vectors jump all the way to the right instead of straddling the dateline.
geoplot(gx,lat1(k),lon1(k));
How can I massage my NetCDF data to prevent this?
My goal is to collect all the patches and merge them into a single surface as described here:
https://www.mathworks.com/matlabcentral/answers/80524-combine-surfaces-into-one-big-surfaceI am trying to accomplish something similar to the below post, namely prevent my data from wrapping at the meridians. So far, no luck:
https://www.mathworks.com/matlabcentral/answers/2172159-how-to-wrap-axes-using-geoplot
If I set the lonLim to [-360 360] and call geolimits(latlim,lonlim), the longitude limits do NOT change to [-360 360] but to [-180.3914 181.1743]. Why is this? Is there a bult-in max limit in geolimits? It’s possible that the function is confused by the [-360 360] values, which are effectively equivalent.
But back to the wraparound issue. I am plotting NetCDF-type satellite observation data.
figure
gx = geoaxes;
geobasemap(gx,’usgsimagery’);
hold on
latLim = [-72 85];
lonLim = [-180 180];
geolimits(gx, latLim, lonLim);
lat = ncread(filepath, ‘Latitude’); % read NetCDF-type data
lon = ncread(filepath, ‘Longitude’);
lat1 = double(lat(:,:,1)) ; % read just the 1st of 6 layers
lon1 = double(lon(:,:,1));
lon1 = wrapTo180(lon1); % limit to +/- 180 degrees
k = boundary(lat1(:),lon1(:)); % get just the edge data
pgon = geopolyshape(lat1(k),lon1(k));
geoplot(gx,pgon);
I am able to plot the satellite observation patches on the Mercator projection geoaxes display. However, when a patch crosses the International Date Line, it does not straddle the line but whangs all the way to the right, to the Prime Meridian.
If I use geoplot instead of geopolyshape I get the same effect, only with vectors instead of shapes.When the longitude data reaches -180, the vectors jump all the way to the right instead of straddling the dateline.
geoplot(gx,lat1(k),lon1(k));
How can I massage my NetCDF data to prevent this?
My goal is to collect all the patches and merge them into a single surface as described here:
https://www.mathworks.com/matlabcentral/answers/80524-combine-surfaces-into-one-big-surface I am trying to accomplish something similar to the below post, namely prevent my data from wrapping at the meridians. So far, no luck:
https://www.mathworks.com/matlabcentral/answers/2172159-how-to-wrap-axes-using-geoplot
If I set the lonLim to [-360 360] and call geolimits(latlim,lonlim), the longitude limits do NOT change to [-360 360] but to [-180.3914 181.1743]. Why is this? Is there a bult-in max limit in geolimits? It’s possible that the function is confused by the [-360 360] values, which are effectively equivalent.
But back to the wraparound issue. I am plotting NetCDF-type satellite observation data.
figure
gx = geoaxes;
geobasemap(gx,’usgsimagery’);
hold on
latLim = [-72 85];
lonLim = [-180 180];
geolimits(gx, latLim, lonLim);
lat = ncread(filepath, ‘Latitude’); % read NetCDF-type data
lon = ncread(filepath, ‘Longitude’);
lat1 = double(lat(:,:,1)) ; % read just the 1st of 6 layers
lon1 = double(lon(:,:,1));
lon1 = wrapTo180(lon1); % limit to +/- 180 degrees
k = boundary(lat1(:),lon1(:)); % get just the edge data
pgon = geopolyshape(lat1(k),lon1(k));
geoplot(gx,pgon);
I am able to plot the satellite observation patches on the Mercator projection geoaxes display. However, when a patch crosses the International Date Line, it does not straddle the line but whangs all the way to the right, to the Prime Meridian.
If I use geoplot instead of geopolyshape I get the same effect, only with vectors instead of shapes.When the longitude data reaches -180, the vectors jump all the way to the right instead of straddling the dateline.
geoplot(gx,lat1(k),lon1(k));
How can I massage my NetCDF data to prevent this?
My goal is to collect all the patches and merge them into a single surface as described here:
https://www.mathworks.com/matlabcentral/answers/80524-combine-surfaces-into-one-big-surface geoplot, geopolyshape, netcdf, geoaxes, geolimits MATLAB Answers — New Questions
Error: following inductors are connected with the following voltage source error matlab
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
Writing to table in matlab does not appear in a shared workbook?
Trying to write to a excel workbook that has shared privileges. The shared workbook works when two seperate people open it and make changes. Where there is conflict a box appears that allows the users to accept the changes etc…
When i write on the same line using matlab, then use the write command to save it. The changes do not appear on the opened workbook when clicking save. It just saves whatever is in the current workbook at that time with 0 conflicts, even though matlab has already wrote to it.
You only see the changes when you close the workbook and open it again without hitting save.
Is there a way round this ?Trying to write to a excel workbook that has shared privileges. The shared workbook works when two seperate people open it and make changes. Where there is conflict a box appears that allows the users to accept the changes etc…
When i write on the same line using matlab, then use the write command to save it. The changes do not appear on the opened workbook when clicking save. It just saves whatever is in the current workbook at that time with 0 conflicts, even though matlab has already wrote to it.
You only see the changes when you close the workbook and open it again without hitting save.
Is there a way round this ? Trying to write to a excel workbook that has shared privileges. The shared workbook works when two seperate people open it and make changes. Where there is conflict a box appears that allows the users to accept the changes etc…
When i write on the same line using matlab, then use the write command to save it. The changes do not appear on the opened workbook when clicking save. It just saves whatever is in the current workbook at that time with 0 conflicts, even though matlab has already wrote to it.
You only see the changes when you close the workbook and open it again without hitting save.
Is there a way round this ? matlab, excel, writetable MATLAB Answers — New Questions
For a parent variant + referenced subsystem with a child variant subsystem, can I have the variant control choice for the children be unique for each referenced parent?
Hi,
I have a parent variant subsystem with 2 choices. This is also a referenced subsystem and is used in 3 instances, let’s name them: A, B and C.
Within each parent there is a child variant subsystem with also two choices, ie. Parent_A -> Parent_A_Choice_X -> Child -> Child_Choice_Y
, where X is the variant choice for the parent subsystem and Y is the variant choice for the child subsystem.
My isssue is, that when I define the variant choice conditions (I do this in expression mode) for one of the children, the choice propagates to all the referenced parent subsystems A, B and C. The functionality I would like is that I can choose unique expressions for the children in each of the parent subsystems.
Is there any way to do this?
Best regards,
Emil MunkHi,
I have a parent variant subsystem with 2 choices. This is also a referenced subsystem and is used in 3 instances, let’s name them: A, B and C.
Within each parent there is a child variant subsystem with also two choices, ie. Parent_A -> Parent_A_Choice_X -> Child -> Child_Choice_Y
, where X is the variant choice for the parent subsystem and Y is the variant choice for the child subsystem.
My isssue is, that when I define the variant choice conditions (I do this in expression mode) for one of the children, the choice propagates to all the referenced parent subsystems A, B and C. The functionality I would like is that I can choose unique expressions for the children in each of the parent subsystems.
Is there any way to do this?
Best regards,
Emil Munk Hi,
I have a parent variant subsystem with 2 choices. This is also a referenced subsystem and is used in 3 instances, let’s name them: A, B and C.
Within each parent there is a child variant subsystem with also two choices, ie. Parent_A -> Parent_A_Choice_X -> Child -> Child_Choice_Y
, where X is the variant choice for the parent subsystem and Y is the variant choice for the child subsystem.
My isssue is, that when I define the variant choice conditions (I do this in expression mode) for one of the children, the choice propagates to all the referenced parent subsystems A, B and C. The functionality I would like is that I can choose unique expressions for the children in each of the parent subsystems.
Is there any way to do this?
Best regards,
Emil Munk simulink, variant subsystem, referenced subsystem, variant choice MATLAB Answers — New Questions
Band pass filtering of time series data
Hello, I am trying to replicate this study that filters the time series of the H-component of the magnetic field in the period range of 10-45 seconds to see ultralow frequency variations around the time of an earthquake or volcanic event, so may I ask anyone how to use the bandpass function in MATLAB right in this case?
Here’s my silly attempt
fs = 1/60; %since the data is per minute
bandpass(T.DAVH,[0.022 0.1],fs); % period range of 10-45 seconds to frequency in Hz
but it returns
Warning: Specified passband frequency is beyond the Nyquist range so signal has been
filtered with an allstop filter.
> In highpass>designFilter (line 129)
In highpass (line 97)
In bandpass>performHighpassFiltering (line 149)
In bandpass (line 116)
Attached the data below where T.DAVH is the H-component of the magnetic field measured in Davao station.
The description of the data is as follows:
Format IAGA-2002 |
Source of Data Kyushu University (KU) |
Station Name Davao |
IAGA CODE DAV (KU code) |
Geodetic Latitude 007.000 |
Geodetic Longitude 125.400 |
Elevation 8888.88 |
Reported HDZF |
Sensor Orientation HDZ |
Digital Sampling 1 seconds |
Data Interval Type Averaged 1-minute (00:30 – 01:29) |
Data Type Provisional
Thank you to anyone who can help!Hello, I am trying to replicate this study that filters the time series of the H-component of the magnetic field in the period range of 10-45 seconds to see ultralow frequency variations around the time of an earthquake or volcanic event, so may I ask anyone how to use the bandpass function in MATLAB right in this case?
Here’s my silly attempt
fs = 1/60; %since the data is per minute
bandpass(T.DAVH,[0.022 0.1],fs); % period range of 10-45 seconds to frequency in Hz
but it returns
Warning: Specified passband frequency is beyond the Nyquist range so signal has been
filtered with an allstop filter.
> In highpass>designFilter (line 129)
In highpass (line 97)
In bandpass>performHighpassFiltering (line 149)
In bandpass (line 116)
Attached the data below where T.DAVH is the H-component of the magnetic field measured in Davao station.
The description of the data is as follows:
Format IAGA-2002 |
Source of Data Kyushu University (KU) |
Station Name Davao |
IAGA CODE DAV (KU code) |
Geodetic Latitude 007.000 |
Geodetic Longitude 125.400 |
Elevation 8888.88 |
Reported HDZF |
Sensor Orientation HDZ |
Digital Sampling 1 seconds |
Data Interval Type Averaged 1-minute (00:30 – 01:29) |
Data Type Provisional
Thank you to anyone who can help! Hello, I am trying to replicate this study that filters the time series of the H-component of the magnetic field in the period range of 10-45 seconds to see ultralow frequency variations around the time of an earthquake or volcanic event, so may I ask anyone how to use the bandpass function in MATLAB right in this case?
Here’s my silly attempt
fs = 1/60; %since the data is per minute
bandpass(T.DAVH,[0.022 0.1],fs); % period range of 10-45 seconds to frequency in Hz
but it returns
Warning: Specified passband frequency is beyond the Nyquist range so signal has been
filtered with an allstop filter.
> In highpass>designFilter (line 129)
In highpass (line 97)
In bandpass>performHighpassFiltering (line 149)
In bandpass (line 116)
Attached the data below where T.DAVH is the H-component of the magnetic field measured in Davao station.
The description of the data is as follows:
Format IAGA-2002 |
Source of Data Kyushu University (KU) |
Station Name Davao |
IAGA CODE DAV (KU code) |
Geodetic Latitude 007.000 |
Geodetic Longitude 125.400 |
Elevation 8888.88 |
Reported HDZF |
Sensor Orientation HDZ |
Digital Sampling 1 seconds |
Data Interval Type Averaged 1-minute (00:30 – 01:29) |
Data Type Provisional
Thank you to anyone who can help! signal processing, bandpass filtering, bandpass MATLAB Answers — New Questions
Keeping the scale in subplots when on is changed
Hi, I have the following plot. When I usex the option "axis aquare" for figure (b), Figure (c) does not scale properly with figure (b). There is away how to keep the scale similar in the other figures despite I changed the properties of one of them. Thank you for your help.Hi, I have the following plot. When I usex the option "axis aquare" for figure (b), Figure (c) does not scale properly with figure (b). There is away how to keep the scale similar in the other figures despite I changed the properties of one of them. Thank you for your help. Hi, I have the following plot. When I usex the option "axis aquare" for figure (b), Figure (c) does not scale properly with figure (b). There is away how to keep the scale similar in the other figures despite I changed the properties of one of them. Thank you for your help. keeping scale properties of a subplots MATLAB Answers — New Questions
How to find the RMSE?
I have the following code:
clear;clc
ula = phased.ULA(‘NumElements’,10,’ElementSpacing’,0.5);
angs = [40 -20 20; 0 0 0];% angs=[azimuth; elevation]; My desired angles are 40, -20 and 20.
NumSignals = size(angs,2);
c = physconst(‘LightSpeed’);
fc = 300e6; % Operating frequency
lambda = c/fc;
pos = getElementPosition(ula)/lambda;
Nsamp = 1000;
nPower = 0.01;
rs = rng(2007);
signal = sensorsig(pos,Nsamp,angs,nPower);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MUSIC Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
musicspatialspect = phased.MUSICEstimator(‘SensorArray’,ula,…
‘OperatingFrequency’,fc,’ScanAngles’,-90:90,…
‘DOAOutputPort’,true,’NumSignalsSource’,’Property’,’NumSignals’,NumSignals);
[~,ang] = musicspatialspect(signal)
ymusic = musicspatialspect(signal);
helperPlotDOASpectra(musicspatialspect.ScanAngles,ymusic)% Changed by Me
function helperPlotDOASpectra(x2,y2)% Changed by Me
% Plot spectra in dB normalized to 0 dB at the peak
y2_dB = 20*log10(y2) – max(20*log10(y2));
plot(x2,y2_dB)
xlabel(‘Broadside Angle (degrees)’);
ylabel(‘Power (dB)’);
title(‘DOA Spatial Spectra’)
end
This estimates the angles and displays them in the command window. It also gives me the plot. Now I want to find the RMSE plot i.e., the root mean square error vs number of runs plot , where number of runs=100. The RMSE is between the actual angles i.e., 40 -20 20 assigend to "angs" and the estimated angles in variable "ang". But I don’t know how to do it?I have the following code:
clear;clc
ula = phased.ULA(‘NumElements’,10,’ElementSpacing’,0.5);
angs = [40 -20 20; 0 0 0];% angs=[azimuth; elevation]; My desired angles are 40, -20 and 20.
NumSignals = size(angs,2);
c = physconst(‘LightSpeed’);
fc = 300e6; % Operating frequency
lambda = c/fc;
pos = getElementPosition(ula)/lambda;
Nsamp = 1000;
nPower = 0.01;
rs = rng(2007);
signal = sensorsig(pos,Nsamp,angs,nPower);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MUSIC Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
musicspatialspect = phased.MUSICEstimator(‘SensorArray’,ula,…
‘OperatingFrequency’,fc,’ScanAngles’,-90:90,…
‘DOAOutputPort’,true,’NumSignalsSource’,’Property’,’NumSignals’,NumSignals);
[~,ang] = musicspatialspect(signal)
ymusic = musicspatialspect(signal);
helperPlotDOASpectra(musicspatialspect.ScanAngles,ymusic)% Changed by Me
function helperPlotDOASpectra(x2,y2)% Changed by Me
% Plot spectra in dB normalized to 0 dB at the peak
y2_dB = 20*log10(y2) – max(20*log10(y2));
plot(x2,y2_dB)
xlabel(‘Broadside Angle (degrees)’);
ylabel(‘Power (dB)’);
title(‘DOA Spatial Spectra’)
end
This estimates the angles and displays them in the command window. It also gives me the plot. Now I want to find the RMSE plot i.e., the root mean square error vs number of runs plot , where number of runs=100. The RMSE is between the actual angles i.e., 40 -20 20 assigend to "angs" and the estimated angles in variable "ang". But I don’t know how to do it? I have the following code:
clear;clc
ula = phased.ULA(‘NumElements’,10,’ElementSpacing’,0.5);
angs = [40 -20 20; 0 0 0];% angs=[azimuth; elevation]; My desired angles are 40, -20 and 20.
NumSignals = size(angs,2);
c = physconst(‘LightSpeed’);
fc = 300e6; % Operating frequency
lambda = c/fc;
pos = getElementPosition(ula)/lambda;
Nsamp = 1000;
nPower = 0.01;
rs = rng(2007);
signal = sensorsig(pos,Nsamp,angs,nPower);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MUSIC Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
musicspatialspect = phased.MUSICEstimator(‘SensorArray’,ula,…
‘OperatingFrequency’,fc,’ScanAngles’,-90:90,…
‘DOAOutputPort’,true,’NumSignalsSource’,’Property’,’NumSignals’,NumSignals);
[~,ang] = musicspatialspect(signal)
ymusic = musicspatialspect(signal);
helperPlotDOASpectra(musicspatialspect.ScanAngles,ymusic)% Changed by Me
function helperPlotDOASpectra(x2,y2)% Changed by Me
% Plot spectra in dB normalized to 0 dB at the peak
y2_dB = 20*log10(y2) – max(20*log10(y2));
plot(x2,y2_dB)
xlabel(‘Broadside Angle (degrees)’);
ylabel(‘Power (dB)’);
title(‘DOA Spatial Spectra’)
end
This estimates the angles and displays them in the command window. It also gives me the plot. Now I want to find the RMSE plot i.e., the root mean square error vs number of runs plot , where number of runs=100. The RMSE is between the actual angles i.e., 40 -20 20 assigend to "angs" and the estimated angles in variable "ang". But I don’t know how to do it? rmse, root mean sqaure error, music, algorithm MATLAB Answers — New Questions
Minimum of bivariable function
Hi everyone,
I just started using MATLAB and I’m trying to find the minmum of a bivariable function (CT(T1, T2)) and the correspoding 2D plot. Below the code used:
clear all ; clc ;
syms x;
syms T1;
syms T2;
H1=50;
%T2=10;
%T1=11;
Mp1=400;
Mc1=1500;
Mp2= 442.0684;
Mc2=1700;
H=100;
H2=H-H1;
w=0.1;
F1=@(T1)1-exp(-((T1-x)./100)).^2;
F2=@(T2) 1-exp(-((T2-x)./100)).^1.8;
f1=@(T1) 1-exp(-((H1-floor(H1./T1).*T1)-x)./100).^2;
f2=@(T2) 1-exp(-((H2-floor(H2./T2).*T2)-x)./100).^1.8;
A=@(T1)((H1./T1).*((-log(1-int(((30).*exp(-30).*x)*((1-exp(-((T1-x)./100)).^2)),0,T1)).*Mc1+(Mp1))))./H1;
B=@(T1)(-log(1-int((F1*((30).*exp(-30).*x)),0,(H1-(floor(H1./T1)).*T1)).*Mc1)./H1);
C=Mp2./H1;
D=@(T2)((H2./T2).*((-log(1-int(((1-exp(-((T2-x)./100)).^1.8))*((25).*exp(-25).*x),0,T2))).*(Mp1.*exp(w))))./H2;
E=@(T2)((-log(1-int((1-exp(-((H2-((H2./T2)).*T2)-x)./100))*((25).*exp(-25).*x),0,(H2-((H2./T2))).*T2))))*(Mc1.*exp(w))./H2;
CT=@(T1,T2)(A(T1)+B(T1)+C+D(T2)+E(T2))
Thanks.Hi everyone,
I just started using MATLAB and I’m trying to find the minmum of a bivariable function (CT(T1, T2)) and the correspoding 2D plot. Below the code used:
clear all ; clc ;
syms x;
syms T1;
syms T2;
H1=50;
%T2=10;
%T1=11;
Mp1=400;
Mc1=1500;
Mp2= 442.0684;
Mc2=1700;
H=100;
H2=H-H1;
w=0.1;
F1=@(T1)1-exp(-((T1-x)./100)).^2;
F2=@(T2) 1-exp(-((T2-x)./100)).^1.8;
f1=@(T1) 1-exp(-((H1-floor(H1./T1).*T1)-x)./100).^2;
f2=@(T2) 1-exp(-((H2-floor(H2./T2).*T2)-x)./100).^1.8;
A=@(T1)((H1./T1).*((-log(1-int(((30).*exp(-30).*x)*((1-exp(-((T1-x)./100)).^2)),0,T1)).*Mc1+(Mp1))))./H1;
B=@(T1)(-log(1-int((F1*((30).*exp(-30).*x)),0,(H1-(floor(H1./T1)).*T1)).*Mc1)./H1);
C=Mp2./H1;
D=@(T2)((H2./T2).*((-log(1-int(((1-exp(-((T2-x)./100)).^1.8))*((25).*exp(-25).*x),0,T2))).*(Mp1.*exp(w))))./H2;
E=@(T2)((-log(1-int((1-exp(-((H2-((H2./T2)).*T2)-x)./100))*((25).*exp(-25).*x),0,(H2-((H2./T2))).*T2))))*(Mc1.*exp(w))./H2;
CT=@(T1,T2)(A(T1)+B(T1)+C+D(T2)+E(T2))
Thanks. Hi everyone,
I just started using MATLAB and I’m trying to find the minmum of a bivariable function (CT(T1, T2)) and the correspoding 2D plot. Below the code used:
clear all ; clc ;
syms x;
syms T1;
syms T2;
H1=50;
%T2=10;
%T1=11;
Mp1=400;
Mc1=1500;
Mp2= 442.0684;
Mc2=1700;
H=100;
H2=H-H1;
w=0.1;
F1=@(T1)1-exp(-((T1-x)./100)).^2;
F2=@(T2) 1-exp(-((T2-x)./100)).^1.8;
f1=@(T1) 1-exp(-((H1-floor(H1./T1).*T1)-x)./100).^2;
f2=@(T2) 1-exp(-((H2-floor(H2./T2).*T2)-x)./100).^1.8;
A=@(T1)((H1./T1).*((-log(1-int(((30).*exp(-30).*x)*((1-exp(-((T1-x)./100)).^2)),0,T1)).*Mc1+(Mp1))))./H1;
B=@(T1)(-log(1-int((F1*((30).*exp(-30).*x)),0,(H1-(floor(H1./T1)).*T1)).*Mc1)./H1);
C=Mp2./H1;
D=@(T2)((H2./T2).*((-log(1-int(((1-exp(-((T2-x)./100)).^1.8))*((25).*exp(-25).*x),0,T2))).*(Mp1.*exp(w))))./H2;
E=@(T2)((-log(1-int((1-exp(-((H2-((H2./T2)).*T2)-x)./100))*((25).*exp(-25).*x),0,(H2-((H2./T2))).*T2))))*(Mc1.*exp(w))./H2;
CT=@(T1,T2)(A(T1)+B(T1)+C+D(T2)+E(T2))
Thanks. bivariable, optimization MATLAB Answers — New Questions
