Azzera filtri
Azzera filtri

Same code but different result with optimal control - forward-backward sweep method

19 visualizzazioni (ultimi 30 giorni)
Hi all,
I have got some issues with Matlab. I copied a code that I found in the appendix on this website: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7123754/#Equ28
However, I get a slightly different result than as shown in the figure in the article. I also reproduced another optimal control code from another article in which my results differ from those of the article too.
I have a feeling that there are problems with the while loop.
Could it be that the while loop only repeat once? Can I check how convergence happens over time?
I would be nice if someone can help me.
Thanks in advance!
The code is as follows:
function ocmodel1
% This function computes the optimal control
% and the corresponding solution using forward−backward sweep
clc;
clear all;
test = -1;
delta1 = 0.001; %set tolerance
N = 100; %number of subdivisions
h = 1/N; %step
t = 0:h:1; % t−variable mesh
u = zeros(1,length(t)); %initialization
x = zeros(1,length(t));
lam = zeros(1,length(t));
x(1) = 10; %initial value assigned to x(0)
beta = 0.05; %parameters
mu = 0.01;
gamma = 0.5;
P = 100;
w1 = 1;
while (test<0) % while the tolerance is reached, repeat
oldu = u;
oldx = x;
oldlam = lam;
for i=1:N %loop that solve the forward differential equation
k1 = beta*(P-x(i)) * x(i) -(mu + gamma) * x(i) - u(i) * x(i);
k2 = beta*(P-x(i)-0.5 * k1 * h)*(x(i)+0.5 * k1 * h) - (mu+gamma)*(x(i)+0.5 * k1 * h)-0.5*(u(i)+u(i+1))*(x(i)+0.5 * k1 * h);
k3 = beta*(P-x(i)-0.5 * k2 * h)*(x(i)+0.5 * k2 * h) - (mu+gamma)*(x(i)+0.5 * k2 * h)-0.5*(u(i)+u(i+1))*(x(i)+0.5 * k2 * h);
k4 = beta*(P-x(i)-k3 * h)*(x(i)+k3 * h) - (mu+gamma)*(x(i)+k3 * h)-u(i+1)*(x(i)+k3 * h);
x(i+1) = x(i) + (h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N %loop that solves the backward differential equation of the adjoint system
j = N + 2 -i;
k1 = -w1-lam(j)*(beta*(P-x(j))-beta * x(j)-(mu+gamma) - u(j));
k2 = -w1-(lam(j)-0.5 * k1 * h)*(beta*(P-x(j)+0.5 * k1 * h) -(mu+gamma) -0.5*(u(j)+u(j-1)));
k3 = -w1-(lam(j)-0.5 * k2 * h)*(beta*(P-x(j)+0.5 * k2 * h) -(mu+gamma) -0.5*(u(j)+u(j-1)));
k4 = -w1 -(lam(j)-k3 * h)*(beta*(P-x(j)+k3 * h) -(mu+gamma) - u(j-1));
lam(j-1) = lam(j) - (h/6)*(k1+2*k2+2*k3+k4);
end
u1 = min(100,max(0,lam.* x/2));
u = 0.5*(u1 + oldu);
temp1 = delta1 * sum(abs(u)) - sum(abs(oldu - u));
temp2 = delta1 * sum(abs(x)) - sum(abs(oldx - x));
temp3 = delta1 * sum(abs(lam)) - sum(abs(oldlam -lam));
test = min(temp1,min(temp2,temp3));
end
figure(1) %plotting
plot(t,u)
figure(2)
plot(t,x)
end

Risposte (1)

Uday Pradhan
Uday Pradhan il 12 Gen 2021
Hi,
You may want to debug your program with the in-built debugging tools MATLAB provides. You may create a loop - counter variable and increment it inside the loop to count the iterations, also consider printing out the required output to command line to see how the values change.

Categorie

Scopri di più su Loops and Conditional Statements in Help Center e File Exchange

Prodotti


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by