Why I receive a NaN value

Hi everyone, i met a problem when i try to check my following code which is a function file, when i try with some inputs, Matlab generates NaN which made me confused.

%set objective functions
function f =totalcost(y,yi,sij)
%set constants for the equations
r1=0.05;%annual demand growth rate
r2=0.1;%annual demand growth rate caused by the completion of rail routes
d=[0 20 30 35].';%di distance beween links 1-4. d1=0 since link 0-1 does not exist.
ub=10;%cost of buses in $/hour
nc=6;%number of cars per train
uc=50;%cost of rail cars in $/car
L=500;%rail maintenance cost in $/mile
k=1000;%rail construction cost in $/mile
td=0.08;%dwell time for bus and rail in hours
vb=40;% bus operating speed
vr=60;% train operating speed
u=10; %user value of time in dollars per hour
N=4;
qijTemp=randi(10,4);
qij0=tril((qijTemp-diag(diag(qijTemp))),-1)+((tril(qijTemp-diag(diag(qijTemp)),-1))).' ;
qij=zeros(4,4,4);
hb1=zeros(4,1);
hb2=zeros(4,1);
hb=zeros(4,1);
hr=zeros(4,1);
cu=zeros(4,1);
ci=zeros(4,1);
rb=zeros(4,1);
rr=zeros(4,1);
cv=zeros(4,1);
cm=zeros(4,1);
cc=zeros(4,1);
f=zeros(4,1);
for t=1:N
    if t==1
        qij(:,:,t)=qij0;
        yi(:,:,t)=zeros(4);
        sij(:,:,t)=zeros(4);
        cc(t)=y(:,1,t).'*d*k;
    else
        qij(:,:,t)=(qij(:,:,t-1).*(1+r1).^t).*((1/8)*(1+sij(:,:,t-1).*r2).^t);
        cc(t)=(y(:,1,t)-y(:,1,t-1)).'*d*k;%construction cost
    end
    hb1(t)=sum(2*(1-yi(1,:,t))*d./vb)+sum((1-yi(1,:,t))*td)*ub; %bus headway
    hb2(t)=sum(sum((1-yi(:,:,t)).*qij(:,:,t)));
    hb(t)=2*sqrt(hb1(t)+hb2(t));
    hr(t)=2*sqrt(((2*yi(1,:,t)*d)/vr+sum(yi(1,:,t)*td))*nc*uc/sum(sum((0.1+y(:,:,t)).*qij(:,:,t))));
    cu(t)=sum(sum(qij(:,:,t).*(1-yi(:,:,t))))*(hb(t)/2)*(u/4)+sum(sum(qij(:,:,t).*yi(:,:,t)))*(hr(t)/2)*(u/4); %user wait cost
    ci1(t)=((diag((1-yi(:,:,t))*qij(:,:,t))).'*(d+td))*(u/vb);
    ci2(t)=((diag(yi(:,:,t)*qij(:,:,t))).'*(d+td)).*(u+vb);
    ci(t)=ci1(t)+ci2(t);
    rb(t)=2*(1-yi(1,:,t))*d/vb+2*sum((1-yi(1,:,t))*td);%bus round trip time
    rr(t)=2*((yi(1,:,t)*d)/vr)+sum((yi(1,:,t))*td); %rail round trip time
    cv(t)=(rb(t)/hb(t))*ub+(rr(t)/hr(t))*nc*uc;%vehicle operating speed
    cm(t)=(yi(1,:,t)*d)*L;%maintenance cost
    f(t)=cu(t)+ci(t)+cv(t)+cm(t)+cc(t);
end
f = sum(f);
end

The input i tried with is

totalcost(zeros(4,1,4),zeros(4,4,4),zeros(4,4,4))

Thanks for helping me!

[dpb]

Set a breakpoint in the function and step through to see where your logic error lies...

[BOWEN LI]

Thank you for noticing me about this!

[Joel Handy]

Even better than manually placing breakpoints....

dbstop if naninf

That will automatically pause execudtion of the code whenever a NaN or Inf is detected. To turn off this functionality:

dbclear if naninf

dbstop may be the single most helpful function built into Matlab.

[BOWEN LI]

Thank you so much!

[Walter Roberson]

    cv(t)=(rb(t)/hb(t))*ub+(rr(t)/hr(t))*nc*uc;%vehicle operating speed

has potential division by 0 in two places.

[BOWEN LI]

yeah I think rr(t)/hr(t) is a 0/0.

However for

hr(t)=2*sqrt(((2*yi(1,:,t)*d)/vr+sum(yi(1,:,t)*td))*nc*uc/sum(sum((0.1+y(:,:,t)).*qij(:,:,t))));

where

sum(sum((0.1+y(:,:,t)).*qij(:,:,t))

probably has the potential to be 0, but i have not found the solution yet.

Thank you!

[Walter Roberson]

In the context of the problem you are trying to solve, what result would make sense if those denominators were 0?

Sometimes it makes sense to replace 0/0 with 1 because it might make sense to do so under the theory of limits.

[BOWEN LI]

Hi for the denomiator which is

sum(sum((0.1+y(:,:,t)).*qij(:,:,t))

I am trying to multiply two matrices 'y' and 'qij' as my denominator, where 'y' is a binary desicion variable (0-1) matrix has a 4x4x4 dimension (first two dimension are just rows and columns, the third dimension is time which has four time dimensions), for qij is a scalar matrix with the dimension same as 'y'.

so here I guess for each time period t, 'y' probably all zeros, but qij definitely not. So that's why i add 0.1 before y(:,:t) to make sure that the denominator will not become 0. But the result is still NaN.

[Walter Roberson]

The denominator is not the problem here. The problem is that (2*yi(1,:,t)*d)/vr+sum(yi(1,:,t)*td) is 0 so the numerator comes out as 0, making hr(1) be 0.

[BOWEN LI]

Yeah, you are right. Can't thank you enough.