Adjoint function in an optimal control problem returns NaN values

I am performing optimal control simulations for the first time, but almost all entries of control functions (vectors) u1 and u2 go to the upper bound which is 1. I found out that this is because the adjoint functions lamda1, lambda2, ..., lambda11 return NaN values. I already carefully checked the adjoint equations to check for possible singularities. I suspect that it's because of the looping. What else could be the possible reason and how could I fix it? Thank you so much in advance.

close all
clf
clear
clc
beta1 = 0.2945;
beta2 = 0.3711;
rho = 0.3081;
Lambda = 3148798;
mu = 0.01444;
db = 0.01;
dh = 0.333;
dbh = 0.01;
du = 0.2;
epsilonH = 0.1;
epsilonB = 0.5;
epsilonU = 0.75;
xi = 1.1;
zeta1 = 1.5;
zeta2 = 1.3;
etaA = 1.35;
etaBH = 1.155;
pi=0.006;
tau1=0.33;
tau2=0.3079;
tau3=0.3;
tau4=0.015;
kappa = 0.002;
omega = 0.01;
phi = 0.2;
varphi=0.25;
p = 0.63;
r = 0.025;
T = 20;
N0 = 109465287;
IB0 = 72598;
IH0 = 9688;
AH0 = 1415;
IBH0 = 143;
P0 = 0.001*N0;
TB0=0.1*IB0;
TH0=9405;
TBH0=0.5*IBH0;
U0 = 200000;
S0 = N0-U0-IB0-IH0-AH0-IBH0-P0-TB0-TH0-TBH0;
x0 = [S0; P0; U0; IB0; IH0; AH0; IBH0; TB0; TH0; TBH0; N0];
test = -1;
delta= 0.001;
N = 1000;
t = linspace(2017,2037,N+1);
h = T/N;
h2 = h/2;
u1 = zeros(1,N+1);
u2 = zeros(1,N+1);
x = zeros(11,N+1);
x(:,1) = x0;
lambda = zeros(11,N+1);
lambda(T)=0;
while(test < 0) % test for convergence
    oldu1 = u1;
    oldu2 = u2;
    oldx = x;
    oldlambda = lambda;
    
    %New values of x and lambda
    x = RK4fwd(t,x,u1,u2,Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r,N,h,h2);
    lambda = RK4bwd(t,x,u1,u2,lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r,N,h,h2);
    c1 = 10^6;
    c2 = 10^5;
    N1 = x(1,:) + x(2,:) + x(3,:) + x(4,:) + x(5,:) + x(6,:) + x(7,:) + x(8,:) + x(9,:) + x(10,:);
    
    %Updating controls u1 and u2 using the new x and lambda
    opt1 = ((lambda(5,:)-lambda(1,:)).*epsilonH*beta2.*(x(5,:)+etaA.*x(6,:)+etaBH.*x(7,:)+zeta2.*x(3,:)).*x(1,:)./N1 + (lambda(4,:)-lambda(1,:)).*(1-epsilonH)*epsilonB*beta1.*(x(4,:)+xi.*x(7,:)+zeta1.*x(3,:)).*x(1,:)./N1)./c1;
    u1new = max(0,min(opt1,1));
    u1 = 0.5.*(oldu1 + u1new);
    opt2 = ((lambda(1,:)-lambda(2,:)).*(1-epsilonH)*(1-epsilonB)*(1-epsilonU).*x(1,:))./c2;
    u2new = max(0,min(opt2,1));
    u2 = 0.5*(oldu2 + u2new);
    
    temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
    temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
    temp3 = delta*sum(abs(x)) - sum(abs(oldx - x));
    temp4 = delta*sum(abs(lambda)) - sum(abs(oldlambda - lambda));
    test = min(temp1,min(temp2,min(temp3,temp4)));
end
%Plotting
X = x(3,:)+x(4,:)+x(5,:);
figure(1)
set(gcf, 'Units', 'Normalized', ...
    'OuterPosition', [0, 0, 0.45, 0.75]);
set(gcf,'Color','white')
subplot(3,1,1)
plot(t,X,'LineWidth',3)
%set(gca,'FontSize',24)
%ylabel({'$U+I_{B}+I_{H}$'},'Interpreter','latex')
xlim([2017 2037])
subplot(3,1,2)
plot(t,u1,'LineWidth',3)
%set(gca,'FontSize',24)
%ylabel({'$u_{1}(t)$'},'Interpreter','latex')
xlim([2017 2037])
subplot(3,1,3)
plot(t,u2,'LineWidth',3)
%set(gca,'FontSize',24)
xlabel('Time','Interpreter','latex')
%ylabel({'$u_{2}(t)$'},'Interpreter','latex')
xlim([2017 2037])

%optimal control problem
function dx = state_func(t,x,u1,u2,Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r)
S = x(1);
P = x(2);
U = x(3);
IB = x(4);
IH = x(5);
AH = x(6);
IBH = x(7);
TB = x(8);
TH = x(9);
TBH = x(10);
N1 = S + P + U + IB + AH + IH + IBH + TB + TH + TBH;
dS = Lambda*(1-pi) - (1-u1)*epsilonH*beta2*(IH+etaA*AH+etaBH*IBH+zeta2*U)*S/N1 - (1-u1)*(1-epsilonH)*epsilonB*beta1*(IB+xi*IBH+zeta1*U)*S/N1 - (1-epsilonH)*(1-epsilonB)*epsilonU*(beta1*(IB+xi*IBH+zeta1*U)/N1 + beta2*(IH+etaA*AH+etaBH*IBH+zeta2*U)/N1)*S - (1-epsilonH)*(1-epsilonB)*(1-epsilonU)*u2*S - mu*S;
dP = Lambda*pi + (1-epsilonH)*(1-epsilonB)*(1-epsilonU)*u2*S + r*TB - mu*P;
dU = (1-epsilonH)*(1-epsilonB)*epsilonU*(beta1*(IB+xi*IBH+zeta1*U)/N1 + beta2*(IH+etaA*AH+etaBH*IBH+zeta2*U)/N1)*S - (du+mu)*U;
dIB = (1-u1)*(1-epsilonH)*epsilonB*beta1*(IB+xi*IBH+zeta1*U)*S/N1 + phi*TB - beta2*(IH+etaA*AH+etaBH*IBH+zeta2*U)*IB/N1 - (tau1+db+mu)*IB; 
dIH = (1-u1)*epsilonH*beta2*(IH+etaA*AH+etaBH*IBH+zeta2*U)*S/N1 + omega*TH - beta1*(IB+xi*IBH+zeta1*U)*IH/N1 - (rho+tau2+mu)*IH;
dAH = rho*IH + kappa*rho*TH - beta1*(IB+xi*IBH+zeta1*U)*AH/N1 - (tau3+dh+mu)*AH;
dIBH = beta2*(IH+etaA*AH+etaBH*IBH+zeta2*U)*IB/N1 + beta1*(IB+xi*IBH+zeta1*U)*(IH+AH)/N1 + varphi*TBH - (tau4+dbh+mu)*IBH;
dTB = tau1*IB - (phi+r+mu)*TB;
dTH = tau2*IH + tau3*AH + p*TBH - (omega+kappa*rho+mu)*TH;
dTBH = tau4*IBH - (varphi+p+mu)*TBH;
dN1 = Lambda - du*U - db*IB - dh*IH - dbh*IBH - mu*N1;
dx = [dS; dP; dU; dIB; dIH; dAH; dIBH; dTB; dTH; dTBH; dN1];
end

%Solves states forward in time using initial values and controls
function x = RK4fwd(t,x,u1,u2,Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r,N,h,h2)
h2 = 0.5*h;
for i = 1:N
    k1 = state_func(t,x(:,i),u1(i),u2(i),Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    k2 = state_func(t+h2,x(:,i)+h2.*k1,0.5*(u1(i)+u1(i+1)),0.5*(u2(i)+u2(i+1)),Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    k3 = state_func(t+h2,x(:,i)+h2.*k2,0.5*(u1(i)+u1(i+1)),0.5*(u2(i)+u2(i+1)),Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    k4 = state_func(t+h,x(:,i)+h.*k3,u1(i+1),u2(i+1),Lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,pi,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    x(:,i+1) = x(:,i)+(h/6).*(k1 + 2.*k2 + 2.*k3 + k4);
end
end

% adjoint system
function dlambda = lambda_func(t,x,u1,u2,lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r)
N1 = x(1,:) + x(2,:) + x(3,:) + x(4,:) + x(5,:) + x(6,:) + x(7,:) + x(8,:) + x(9,:) + x(10,:);
lambda1 = lambda(1);
lambda2 = lambda(2);
lambda3 = lambda(3);
lambda4 = lambda(4);
lambda5 = lambda(5);
lambda6 = lambda(6);
lambda7 = lambda(7);
lambda8 = lambda(8);
lambda9 = lambda(9);
lambda10 = lambda(10);
lambda11 = lambda(11);
dlambda1 = (lambda1-lambda2)*(1-epsilonH)*(1-epsilonB)*(1-epsilonU)*u2 + (lambda1-lambda3)*(1-epsilonH)*(1-epsilonB)*epsilonU*((beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))/N1) + beta2*(x(5,:)+etaA*x(6,:)+etaBH*x(7,:)+zeta2*x(3,:))/N1) + (lambda1-lambda4)*(1-u1)*(1-epsilonH)*epsilonB*beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))/N1 + (lambda1-lambda5)*(1-u1)*epsilonH*beta2*(x(5,:)+etaA*x(6,:)+etaBH*x(7,:)+zeta2*x(3,:))/N1;
dlambda2 = lambda2*mu;
dlambda3 = (lambda1-lambda5)*(1-u1)*epsilonH*beta2*zeta2*x(1,:)/N1 + (lambda1-lambda4)*(1-u1)*(1-epsilonH)*epsilonB*beta1*zeta1*x(1,:)/N1 + (lambda1-lambda3)*(1-epsilonH)*(1-epsilonB)*epsilonU*(beta1*zeta1+beta2*zeta2)*x(1,:)/N1 + (lambda4-lambda7)*beta2*zeta2*x(4,:)/N1 + (lambda5-lambda7)*beta1*zeta1*x(5,:)/N1 + (lambda6-lambda7)*beta1*zeta1*x(6,:)/N1 + lambda3*(du+mu) + lambda11*du - 1;
dlambda4 = (lambda1-lambda3)*(1-epsilonH)*(1-epsilonB)*epsilonU*beta1*x(1,:)/N1 + (lambda1-lambda4)*(1-u1)*(1-epsilonH)*epsilonB*beta1*x(1,:)/N1 + (lambda5-lambda7)*beta1*x(5,:)/N1 + (lambda6-lambda8)*beta1*x(6,:)/N1 + (lambda4-lambda7)*beta2*(x(5,:)+etaA*x(6,:)+etaBH*x(7,:)+zeta2*x(3,:))*x(1,:)/N1 + (lambda4-lambda8)*tau1 + lambda4*(db+mu) + lambda11*db - 1;
dlambda5 = (lambda1-lambda5)*(1-u1)*epsilonH*beta2*x(1,:)/N1 + (lambda1-lambda3)*(1-epsilonH)*(1-epsilonB)*epsilonU*beta2*x(1,:)/N1 + (lambda4-lambda7)*beta2*x(4,:)/N1 + (lambda5-lambda7)*beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))*x(1,:)/N1 + (lambda5-lambda6)*rho + (lambda5-lambda9)*tau2 + lambda5*mu - 1;
dlambda6 = (lambda1-lambda5)*(1-u1)*epsilonH*beta2*etaA*x(1,:)/N1 + (lambda1-lambda3)*(1-epsilonH)*(1-epsilonB)*epsilonU*beta2*etaA*x(1,:)/N1 + (lambda4-lambda7)*beta2*etaA*x(4,:)/N1 + (lambda6-lambda7)*beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))*x(1,:)/N1 + (lambda6-lambda9)*tau3 + lambda6*(dh+mu) + lambda11*dh;
dlambda7 = (lambda1-lambda5)*(1-u1)*epsilonH*beta2*etaBH*x(1,:)/N1 + (lambda1-lambda4)*(1-u1)*(1-epsilonH)*epsilonB*beta1*xi*x(1,:)/N1 + (lambda1-lambda3)*(1-epsilonH)*(1-epsilonB)*epsilonU*beta2*etaBH*x(1,:)/N1 + (lambda4-lambda7)*beta2*etaBH*x(4,:)/N1 + (lambda5-lambda7)*beta1*xi*x(5,:)/N1 + (lambda6-lambda7)*beta1*xi*x(6,:)/N1 + (lambda7-lambda10)*tau4 + lambda7*(dbh+mu) + lambda11*dbh;
dlambda8 = (lambda8-lambda2)*r + (lambda8-lambda4)*phi;
dlambda9 = (lambda9-lambda5)*omega + (lambda9-lambda6)*kappa*rho + lambda9*mu;
dlambda10 = (lambda10-lambda7)*varphi + (lambda10-lambda9)*p + lambda10*mu;
dlambda11 = (lambda5-lambda1)*(1-u1)*epsilonH*beta2*(x(5,:)+etaA*x(6,:)+etaBH*x(7,:)+zeta2*x(3,:))*x(1,:)/N1^2 + (lambda4-lambda1)*(1-u1)*(1-epsilonH)*epsilonB*beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))*x(1,:)/N1^2 + (lambda3-lambda1)*(1-epsilonH)*(1-epsilonB)*epsilonU*(beta2*(x(5,:)+etaA*x(6,:)+etaBH*x(7,:)+zeta2*x(3,:))*x(1,:) + beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))*x(1,:))/N1^2 + (lambda7-lambda4)*beta2*(x(5,:)+etaA*x(6,:)+etaBH*x(7,:)+zeta2*x(3,:))*x(4,:)/N1^2 + (lambda7-lambda5)*beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))*x(5,:)/N1^2 + (lambda7-lambda6)*beta1*(x(4,:)+xi*x(7,:)+zeta1*x(3,:))*x(6,:)/N1^2 + lambda11*mu;
dlambda = [dlambda1; dlambda2; dlambda3; dlambda4; dlambda5; dlambda6; dlambda7; dlambda8; dlambda9; dlambda10; dlambda11];
end

%Solves adjoint system backward in time using new values of states and the controls
function lambda = RK4bwd(t,x,u1,u2,lambda,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r,N,h,h2)
for i = 1:N
    j = N + 2 - i;
    
    k1 = lambda_func(t,x(:,j),u1(j),u2(j),lambda(:,j),mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    k2 = lambda_func(t-h2,0.5*(x(:,j)+x(:,j-1)),0.5*(u1(j)+u1(j-1)),0.5*(u2(j)+u2(j-1)),lambda(:,j)-h2*k1,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    k3 = lambda_func(t-h2,0.5*(x(:,j)+x(:,j-1)),0.5*(u1(j)+u1(j-1)),0.5*(u2(j)+u2(j-1)),lambda(:,j)-h2*k2,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    k4 = lambda_func(t-h,x(:,j-1),u1(j-1),u2(j-1),lambda(:,j)-h*k3,mu,db,dh,dbh,du,beta1,beta2,xi,zeta1,zeta2,etaA,etaBH,epsilonH,epsilonB,epsilonU,rho,tau1,tau2,tau3,tau4,kappa,omega,phi,varphi,p,r);
    lambda(:,j-1) = lambda(:,j) - (h/6).*(k1 + 2.*k2 + 2.*k3 + k4);
end
end

[Torsten]

Why don't you start with a simpler problem from which you know the solution to get a feeling for what is going on ?

[John D'Errico]

Just because a large array is not mathematically singular, does not mean the same system of equations will not be numerically singular in double precision. The classic example lies in the Hilbert matrix, which is known to be non-singular, yet is also known to be highly ill-conditioned.

A = hilb(5)

A = 5×5

1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cond(A)

ans = 4.7661e+05

rank(A)

ans = 5

which is not too bad. But that is only a 5x5 matrix.

A = hilb(13)

A = 13×13

1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.0500 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.0500 0.0476 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.0500 0.0476 0.0455 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.0500 0.0476 0.0455 0.0435 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.0500 0.0476 0.0455 0.0435 0.0417 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 0.0500 0.0476 0.0455 0.0435 0.0417 0.0400

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cond(A)

ans = 1.3259e+19

rank(A)

ans = 11

How can you fix it? If the problem is numerically singular, then you work with smaller systems of equations. You solve simpler problems. Or you decide to work in higher precision arithmetic, which will never make you happy either because of the slowness of working in high precision. Or you find a different way to solve the problem.

[Gerald]

Thank you for your comments. I already tried to solve a simpler optimal control problem with 3 equations. It worked using the same method. I am able to produce graphs of the control functions that do not reach the upper bound (i.e., 1) since the adjoint functions (lambdas) are not returning NaN values. In addition, I think the reason why the adjoint functions in the more complex problem I posted above are returning NaN values is because the values are very large. For example, lambda1 values are multiples of 1.0e+275. Please see attached image.

What's your take on this?

[Torsten]

Maybe scaling your problem such that the lambda values get normalized to smaller values ?

Maybe the stepsize of your integration is too large such that the explixit method diverges ?

Maybe using a MATLAB ode integrator (ode45,ode15s) instead of the self-written RK4 code ?

Maybe using a solver for boundary value problems (bvp4c) instead of an initial value solver ?

We only have your code and neither a mathematical formulation of your problem nor a description of your solution method. So it's difficult to give advice.