MATLAB Answers

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

11 views (last 30 days)
Gijs D
Gijs D on 8 Jan 2021
Answered: Uday Pradhan on 12 Jan 2021
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

Answers (1)

Uday Pradhan
Uday Pradhan on 12 Jan 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.

Community Treasure Hunt

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

Start Hunting!