Hi all,
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
