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

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

Answer 1

[Uday Pradhan]

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.