27 views (last 30 days)

Show older comments

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

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.

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

Start Hunting!