Adjoint function in an optimal control problem returns NaN values

Question

0 개 추천

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

% 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

댓글 수: 4
이전 댓글 2개 표시 이전 댓글 2개 숨기기

John D'Errico 2025년 10월 4일

편집: John D'Errico 2025년 10월 4일

MATLAB Online에서 열기

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 2025년 10월 9일

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 2025년 10월 9일

편집: Torsten 2025년 10월 9일

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.

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

Follow Question