NMinimize stops with no more memory?

Question

I am using NMinimize function for simulation based optimization. So my objective function is a simulation that runs for every combination of variable values evaluated by NMinimize function. However, the problem I have is the Nminimize function ceases after the first run (I am printing the time stamp for every iteration) and eventually after a long time gives me out of memory error. I even tried with different methods such as "RandomSearch" and "SimulatedAnnealing" with custom method parameter values, but in vain. Can some one pinpoint where I am going wrong?

edit: My code is long, but as requested is given below:

f[a1_, a2_, a3_] := Module[{b1 = a1, b2 = a2, b3 = a3, L = 3, Flen = 1, Rlen = 1,SimTime = 60, Kj = 150,w = 20,Theta = 5, dt = 6,delta = 1,DemandDuration = 10,RMstart = 1,RMLocation = 3, TT = 0}, 

Print[DateString[]]; Vf = Theta w; dx = Vf dt/3600; capacity = w*Vf*Kj/(Vf + w);n = Round[Flen/dx];m = Round[SimTime/dt];p = Round[Rlen/dx];Rdensity = Table[0*i, {i, p}, {i, m}, {i, n}];Rflow = Table[0*i, {i, p}, {i, m}, {i, n}];Fdensity = Table[0*i, {i, n}, {i, m}];Fflow = Table[0*i, {i, n}, {i, m}];demand[n_, k_] := Min[k*Vf, n*capacity];supply[n_, k_] := Min[(n*Kj - k)*w, n*capacity];flo[demand_, supply_] := Min[demand, supply];den[k_, qin_, qout_] := k + (qin - qout)/Vf;
merge[n_, Fu_, Fd_, Rd_] := Min[1, supply[n, Fd]/(demand[n, Fu] + 0.01)]*demand[1, Rd]/delta;Nsupply[n_, k_, qsum_] := Min[(n*Kj - k)*Vf - qsum, n*capacity];
RM[x_, t_] := N[b1 x^2 + b2 x + b3];alpha[a_] := 1500*a/Flen;beta[a_] := 0.1*a/Flen;

For[k = 1, k <= n, k++, 
For[i = 1, i <= p, i++, Rdensity[[i, 1, k]] = 0;];
For[j = 1, j <= DemandDuration, j++, Rdensity[[p, j, k]] = alpha[n*dx]*delta/Vf;
TT = TT + Rdensity[[p, j, k]]];];

For[j = 1, j <= 4, j++, 
For[k = 1, k <= n, k++, 
If[k == 1, 
  Rflow[[1, j, k]] = 
   merge[L, Fdensity[[k, j]], Fdensity[[k, j]], 
    Rdensity[[1, j, k]]], 
  Rflow[[1, j, k]] = 
   merge[L, Fdensity[[k, j]], Fdensity[[k - 1, j]], 
    Rdensity[[1, j, k]]]];
 If[k == 1, 
  Fflow[[k, j]] = 
   flo[demand[L, Fdensity[[k, j]]], supply[L, Fdensity[[k, j]]]], 
  Fflow[[k, j]] = 
   flo[demand[L, Fdensity[[k, j]]], 
    supply[L, Fdensity[[k - 1, j]]] - Rflow[[1, j, k - 1]]*dx]];
 If[k > 1 && j < m, 
  Fdensity[[k - 1, j + 1]] = 
   den[Fdensity[[k - 1, 
     j]], (Rflow[[1, j, k - 1]] - beta[(n)*dx]*Fflow[[k, j]])*dx +
      Fflow[[k, j]], Fflow[[k - 1, j]]];
  TT = TT + Fdensity[[k - 1, j + 1]];];
 If[k == n && j < m, 
  Fdensity[[k, j + 1]] = 
   den[Fdensity[[k, j]], Rflow[[1, j, k]]*dx, Fflow[[k, j]]];
  TT = TT + Fdensity[[k, j + 1]];];
 For[i = 2, i <= p, i++, 
  If[i == RMLocation && j >= RMstart, 
   Rflow[[i, j, k]] = 
    Min[RM[k dx, j dt], 
     flo[demand[1, Rdensity[[i, j, k]]], 
      supply[1, Rdensity[[i - 1, j, k]]]]], 
   Rflow[[i, j, k]] = 
    flo[demand[1, Rdensity[[i, j, k]]], 
     supply[1, Rdensity[[i - 1, j, k]]]]];
  If[j < m, 
   Rdensity[[i - 1, j + 1, k]] = 
    den[Rdensity[[i - 1, j, k]], Rflow[[i, j, k]], 
     Rflow[[i - 1, j, k]]];
   TT = TT + Rdensity[[i - 1, j + 1, k]];]];];];
For[j = 5, j <= m, j++, 
For[k = 1, k <= n, k++, 
If[k == 1, 
 Rflow[[1, j, k]] = 
  merge[L, Fdensity[[k, j]], Fdensity[[k, j]], 
   Rdensity[[1, j, k]]], 
 Rflow[[1, j, k]] = 
  merge[L, Fdensity[[k, j]], Fdensity[[k - 1, j]], 
   Rdensity[[1, j, k]]]];
FQsum = 0;
For[r = 1, r <= Theta - 1, r++, FQsum = FQsum + Fflow[[k, j - r]]];
If[k == 1, 
 Fflow[[k, j]] = 
  flo[demand[L, Fdensity[[k, j]]], supply[L, Fdensity[[k, j]]]], 
 Fflow[[k, j]] = 
  flo[demand[L, Fdensity[[k, j]]], 
   Nsupply[L, Fdensity[[k - 1, j - Theta + 1]], FQsum] - 
    Rflow[[1, j, k - 1]]*dx]];
If[k > 1 && j < m, 
 Fdensity[[k - 1, j + 1]] = 
  den[Fdensity[[k - 1, 
    j]], (Rflow[[1, j, k - 1]] - beta[(n)*dx]*Fflow[[k, j]])*dx + 
    Fflow[[k, j]], Fflow[[k - 1, j]]];
 TT = TT + Fdensity[[k - 1, j + 1]];];
If[k == n && j < m, 
 Fdensity[[k, j + 1]] = 
  den[Fdensity[[k, j]], Rflow[[1, j, k]]*dx, Fflow[[k, j]]];
 TT = TT + Fdensity[[k, j + 1]];];
For[i = 2, i <= p, i++, RQsum = 0;
 For[r = 1, r <= Theta - 1, r++, 
  RQsum = RQsum + Rflow[[i, j - r, k]]];
 If[i == RMLocation && j >= RMstart, 
  Rflow[[i, j, k]] = 
   Min[RM[k dx, j dt], 
    flo[demand[1, Rdensity[[i, j, k]]], 
     Nsupply[1, Rdensity[[i - 1, j - Theta + 1, k]], RQsum]]], 
  Rflow[[i, j, k]] = 
   flo[demand[1, Rdensity[[i, j, k]]], 
    Nsupply[1, Rdensity[[i - 1, j - Theta + 1, k]], RQsum]]];
 If[j < m, 
  Rdensity[[i - 1, j + 1, k]] = 
   den[Rdensity[[i - 1, j, k]], Rflow[[i, j, k]], 
    Rflow[[i - 1, j, k]]];
  TT = TT + Rdensity[[i - 1, j + 1, k]];]];];];

TT]

NMinimize[{f[x, y, z], {x, y,z} \[Element] Integers}, {{x , 27, 30}, {y, 797, 800}, {z, 2497, 2500}}];

Ps. Also any suggestions to improve the performance of this code will be greatly appreciated!!

Answer 1

The solution to release memory space is redefine functions (Fdensity, Fflow, etc) conveniently in order to reduce the number of excessive functional calls in the execution. Using the Alexey Popkov's routine ("Profiling memory usage in Mathematica") we have the following result for the allocated bytes:

 - 6739593360    Fdensity
 - 1984112096    Fflow
 - 1679274128    Rdensity
 - 1453549616    Rflow
 - 37327744      FQsum
 - 977560        RQsum

Also, it is necessary the use of _?NumericQ restriction in the functions to avoid inconvenient symbolic processing.