Group Abstract Group Abstract

Message Boards Message Boards

Optimize Function Module using NMinimize

GROUPS:
I defined a simulation as a module with simulation parameters as the variables. I want to find optimal simulation parameters using NMinimize method. I noticed that my simulation module will run successfully for any sample parameters, but when used in NMinimize will not execute. Please see below a very simplified version of the code:
 demand[n_, k_] := Min[k Vf, n capacity];
 supply[n_, k_] := Min[(n Kj - k) w, n capacity];
 flo[n_, Ku_, Kd_] := Min[demand[n, Ku], supply[n, Kd]];
 dx = Vf*dt; capacity =w*Vf*Kj/(Vf + w); Kj = 150.; w = 20.; Vf = 100.;
 n = Round[Flen/dx]; m = Round[SimTime/dt]; p = Round[Rlen/dx]; RMLocation = Round[(2/3) p];
 \[Alpha][a1_] := 1800.; \[Beta][a2_] := 0.1; L = 1.; Flen = 4.; Rlen = 3.; delta = 1.; SimTime = 15./60.; dt = 6./3600.;
 f[a_] := Module[{k0 = ConstantArray[0, n],kr = Table[Table[0, {i1, 1, p}], {i2, 1, n}], \[Gamma] = ConstantArray[1, n], \[Phi]},
   Clear[j]; j = 0;
   RM[x_, t_] := 100 a; k = k0;
  For[i = 2, i < n, i++, kr[[i, 1]] = \[Alpha][i dx] delta/Vf];
  NtwrkTT = TT = Plus @@ (Plus @@ kr);
  While[TT > 0,
   For[i = 2, i < n, i++,
    FQin = If[i == 2, Min[demand[L, k0[[i - 1]]], supply[L, k0[[i]]]],FQout];
    dem = demand[L, k0[[i]]]; dem = If[dem == 0, 0.001, dem];
    \[Gamma][[i]] = Min[1, supply[L, k0[[i + 1]]]/dem];
    \[Phi] = \[Gamma][[i]] demand[1, kr[[i, p]]]/delta;
    Qr = (\[Phi] - \[Beta][i dx] FQin) dx;
    FQout = Min[demand[L, k0[[i]]], supply[L, k0[[i + 1]]]];
    k[[i]] = k0[[i]] + (FQin - FQout + Qr)/Vf; kr0 = kr[[i]];
    For[ir = 2, ir <= p, ir++,
     MR = If[ir == RMLocation + 1, RM[i dx, j dt], capacity];
     RQin = Min[MR, If[ir == 2, flo[1, kr0[[ir - 1]], kr0[[ir]]], RQout]];
     MR = If[ir == RMLocation, RM[i dx, j dt], capacity];
     RQout = Min[MR, If[ir < p,flo[1, kr0[[ir]], kr0[[ir + 1]]], \[Phi] delta]];
     kr[[i, ir]] = kr0[[ir]] + (RQin - RQout)/Vf];
    kr[[i, 1]] = If[j <= m, \[Alpha][i dx] delta/Vf, 0]];
   TT = Plus @@ (Plus @@ kr);
   TT += Plus @@ k;
   k0 = k; NtwrkTT += TT; j++];
  NtwrkTT dt]
NMinimize[{f[a], 3 <= a <= 12 && Element[a, Integers]}, a, Method -> "SimulatedAnnealing", EvaluationMonitor :> Print["a = ", a]]
POSTED BY: brama gatech
Answer
7 months ago
Numerical functions like NMinimize don't run on problems with undefined parameters.
POSTED BY: Frank Kampas
Answer
7 months ago
@Frank Kampas, Can you please elaborate on undefined parameters? I have only one parameter that is defined with constraints. 
POSTED BY: brama gatech
Answer
7 months ago
What's Vf?  What's Kj?
POSTED BY: Frank Kampas
Answer
7 months ago
They are variables defined as Vf=100 and Kj=150
POSTED BY: brama gatech
Answer
7 months ago
The warning or error looks like it might be resolved if you could fix this
Min[1800., 20. (150. - {18., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}[])]

Min expects a comma separated list of values and it appears that you have a value separated by a comma from a list. You need to track down which of your expressions is giving you a list instead of a simple value.

If you look closely at the end of that you also see [ ] which probably indicates trying to use something incorrectly as a function.

Perhaps you could insert some Print statements to track the values of the expressions you are trying to apply Min to and see if you can pinpoint or at least narrow down where the incorrect expression is.

IF the forum software or the copy-and-paste hasn't corrupted your code then the problem appears to involve kro.
You have kr0[ ] in three places and those are a problem.

Can you correct what those should be?
POSTED BY: Bill Simpson
Answer
7 months ago
@ Bill Simpson, I did some debuging and realized that the forum software messed up the code at 3 palces. I fixed it. I believe the error you got were due to the errors in the code. Will you please look at this updated code?
POSTED BY: brama gatech
Answer
7 months ago
Your revised code seems to run in a second, repeatedly tries various values of a in f and finishes with {5.13, {a -> 3}}

Can you be much more specific about what is incorrect now?

Can you provide much more detailed diagnostics messages showing how the simulation is converging?

If 'a' must be an integer between 3 and 12 then why not just try each of the values of 'a' once and choose the best result?

There are a number of other things about the code which make me worry.

For example, you define k0 to be a ConstantArray, which usually means that k0 will never change, but then you use k0=k and then you Clear. I don't believe there is any need to Clear variables that are local to a Module, they should be automatically removed at the end of the execution of the Module.

There are other parts of the code that are written in a way and that are complicated enough that make me worry that I do not understand it and that it may be doing unexpected tihngs.

Is it possible to make the code substantially simpler?
POSTED BY: Bill Simpson
Answer
7 months ago
@Bill, I apologize for the confusion. Let me clarify things.
The original objective function I want to evaluate has several variables and has a very large solution space. However, to check that Nminimize works properly, I am testing with a smaller solution space with fewer variables and modified parameters (shown in the above code).  
For the above "test parameter set", I know that the optimal happens at a=10. you can check it with
ListPlot[Table[f[a], {a, 3, 12}]]
instead of the Nminimize function.

However, as your would notice, the NMinimize method is giving a=3 as the minimizer. That is the problem.

I do not think constantArray mean it will never change. I have checked it on the hep. I was using "Clear" to make sure the variables are cleared during every iteration. Thanks to your comment, I realize that it is absolete.   Can you please be more specific about what parts of my code are complicated? The "Plus@@" simply adds all the values of the array. Thanks for your attempts to answer my question. 
POSTED BY: brama gatech
Answer
7 months ago
I apologize for my misunderstanding of ConstantArray. I was confusing that with something else.

I understand your real problem is far larger and far more complicated than the quickly simplified version you have shown here. It is not clear to me whether some parts of your code are quickly simplified versions of much more complicated code from your real problem or whether they are written in a style that is not usual for Mathematica or whether they are written in a style that is not as simple and understandable as possible.

Your Module depending on a mixture of internally calculated values and global variables is one example. Your RM function which ignores the parameters it is given is another.

There is no way of knowing whether this applies to your actual problem, but for your simplified example this
 In[65]:= NMinimize[{f[a], 3 <= a <= 12 && Element[a, Integers]}, a,
    Method -> "SimulatedAnnealing", EvaluationMonitor :> Print[{a, f[a]}]]
 
 During evaluation of In[65]:= {3,2398.37}
 During evaluation of In[65]:= {3,2398.37}
 During evaluation of In[65]:= {6,2226.5}
 During evaluation of In[65]:= {7,2221.31}
 During evaluation of In[65]:= {9,2316.97}
 During evaluation of In[65]:= {12,3285.99}
During evaluation of In[65]:= {11,2790.2}
During evaluation of In[65]:= {11,2790.2}
During evaluation of In[65]:= {10,2469.33}
During evaluation of In[65]:= {10,2469.33}
During evaluation of In[65]:= {7,2221.31}
During evaluation of In[65]:= {10,2469.33}
During evaluation of In[65]:= {12,3285.99}
During evaluation of In[65]:= {12,3285.99}
During evaluation of In[65]:= {12,3285.99}
During evaluation of In[65]:= {11,2790.2}
Out[65]= $Aborted
shows that your present use of NMinimize is going to be a horrible way of finding what you want. Clearly NMinimize in the way you are using it was never intended to minimize a list of discrete values.

There also appears to be yet more problems, because you said above that you know the optimal value happen at a=10, but we can see from this list that the minimum value of f[ a ] happens at a=7. I cannot know what that problem is. This also demonstrates that
In[66]:= ListPlot[Table[{a, f[a]}, {a, 3, 12}], AxesOrigin -> {2, 2000}]

Out[66]= ...PlotSnipped...

This
In[67]:= Sort[Table[{a, f[a]}, {a, 3, 12}],
OrderedQ[{#1[[2]], #2[[2]]}] &][[1]]

Out[67]= {7, 2221.31}
will be ten times faster than the way you are trying to use NMinimize in your example.

But again this may have nothing to do with what your real problem is.
POSTED BY: Bill Simpson
Answer
7 months ago
Comment: "It is not clear to me whether some parts of your code are quickly simplified versions of much more complicated code from your real problem or whether they are written in a style that is not usual for Mathematica or whether they are written in a style that is not as simple and understandable as possible. Your Module depending on a mixture of internally calculated values and global variables is one example. Your RM function which ignores the parameters it is given is another. "
Response
: This code is basically a simulation of traffic flow on a network. I load the network with vehicles for "m" timeperiods and then then let the network empty. The objective is to find the optimal control variables (here "a") that minimize the output, NtwrkTT, that is dependent on the control value used. The simulation depends on some global traffic flow parameters such as freeflow speed, capacity, etc. to calculate the traffic densities from which the NetworkTT is calulated.  Thus, the simulation depends on both global parameters and calculated values. I am not sure what you mean by "RM function ignores the parameters". RM function is supposed to be 100 times the "a" as mentioned in the code. 

This is the actual code. I did not simplify the code; only used fewer variables for the objective function and smaller parameter values. 

Comment: There is no way of knowing whether this applies to your actual problem, but for your simplified example this shows that your present use of NMinimize is going to be a horrible way of finding what you want. Clearly NMinimize in the way you are using it was never intended to minimize a list of discrete values.

Response: That is what bothers me.  What I have found is if I donot print the function value (f) during every iteration using the EvaluationMonitor, NMinimize runs 1000s of iterations extremely quickly without running the simulation (because if it runs simulations, it will need few seconds for each simulation run and can not do 1000s of simulations in no time). NMinimize only outputs a bunch of possible values for "a". However, when I include "f" in the evaluation monitor, NMinimize runs the simulation disconnected to the optimization process (I suspect NMinimize barely runs the simulation for every value of "a" it outputs). I found that NMinimize does not converge even after 1000s of iterations, when there are only 10 possible values "a" can take.
Comment:There also appears to be yet more problems, because you said above that you know the optimal value happen at a=10, but we can see from this list that the minimum value of f[ a ] happens at a=7. I cannot know what that problem is.
Response:  I apologize for the software issue where I entered  "RMLocation = Round[(2/3) p];" and it showed "RMLocation = Round[(2/13) p];". with my numbers, you should see that the minimum happens at a=10.
POSTED BY: brama gatech
Answer
7 months ago
I am guessing there must be a good reason why it is not written more simply as RM=100 a; Or perhaps you did not simplify the code and just used fewer variables and it is important that you write MR = If[ ir == RMLocation+1, RM[ i dx, j dt ], capacity ] instead of just MR = If[ ir == RMLocation+1, RM, capacity ].

It is perhap not important, but it might use the same style to do similar things and thus write kr = ConstantArray[ 0, { p, n } ].

I am guessing there must be a good reason why some variables in your Module are listed inside the { } at the beginning of the module, but others are not. The general idea of a Module is that everything it uses is hidden inside and it depends on nothing outside, other than the parameter 'a' that is given to it. In that way someone who is trying to understand it can know they are able to ignore everything else and concentrate only on what is inside the Module. Otherwise it would be as if I had to read the entire program carefully to try to decide whether there might be anything anywhere else that is changing what the value of Sin[ a ] is other than the value of 'a'.

If I use RMLocation = Round[ (2/3) p ]; and
In[ 42 ]:= Table[ { a, f[ a ] }, { a, 3, 12 } ]
Out[ 42 ]= { { 3, 19796.5 } , { 4, 15470.9 } , { 5, 12806. } , { 6, 10969.3 } , { 7, 9603.36 } , { 8, 8528.46 } , { 9, 7858.54 } , { 10, 7609.2 } , { 11, 7771.77 } , { 12, 8102.1 } }
I do get a minimum at n=10

If the real problem you were originally trying to ask was why NMinimize doesn't seem to be iterating correctly that seems like something specific we can actually try to track down. If you can post the current form of your code with all the corrections and with two different uses of NMinimize at the bottom, one "fast" and the other "slow" then perhaps we can figure out what the problem is. It is true that pausing the calculations to display some value can slow things down, but that should not make changes in how it converges. Hopefully we can get through all the small errors and misunderstandings and get to the real problem and solve that.
POSTED BY: Bill Simpson
Answer
7 months ago
Yes. My objective function is complex with multiple variables. therefore, I need to use RM[ i dx, j dt ] instead of RM = 100 a. kr is supposed to be nested tables. So it needs to be kr = Table[Table[0, {i1, 1, p}], {i2, 1, n}] and not kr = ConstantArray[ 0, { p, n } ]. 
I list some variables inside the Module and others outside the Module because the ones outside are global and does not change with changing control, They only need to be defined once and need not be defined every iteration. Hoever, as per your suggestion, I moved everything inside the Module, but it still does not help. 
The code I posted above is the final version of the code and we can use that for debugging. May be we can talk on phone or Gtalk/Skype call, if you think it will help.
POSTED BY: brama gatech
Answer
7 months ago
That was why I kept saying that because I did not know what your real code was that I did not know how much of what I was saying applied to the real problem.
In[1]:= kr = Table[Table[0, {i1, 1, 2}], {i2, 1, 3}]
Out[1]= {{0, 0}, {0, 0}, {0, 0}}

In[2]:= kr = ConstantArray[0, {3, 2}]
Out[2]= {{0, 0}, {0, 0}, {0, 0}}

In[3]:= kr = Table[0, {3}, {2}]
Out[3]= {{0, 0}, {0, 0}, {0, 0}}

There is no difference between a table of table and a two dimensional table in Mathematica, they are exactly the same.

Please post the exact two different NMinimize[ ,,, ] that you see behave differently when used on exactly the code you have posted at the top of this. Please describe exactly what difference you see when you run each of those. Then perhaps we can track down the problem and find a solution.
Thank you
POSTED BY: Bill Simpson
Answer
7 months ago
I found the solution on another website "http://eternaldisturbanceincosmos.wordpress.com/2011/04/27/nminimize-in-mathematica-could-drive-you-insane/" which says "It turns out that NMinimize does not hold its arguments. This means that as the list of arguments is read from left to right, each argument is evaluated and replaced by the result of the evaluation". So I used `Hold[]` in the `NMinimize` function that fixed the problem.I found the solution on another website "http://eternaldisturbanceincosmos.wordpress.com/2011/04/27/nminimize-in-mathematica-could-drive-you-insane/" which says "It turns out that NMinimize does not hold its arguments. This means that as the list of arguments is read from left to right, each argument is evaluated and replaced by the result of the evaluation". So I used `Hold[]` in the `NMinimize` function that fixed the problem.
POSTED BY: Updating Name
Answer
7 months ago
Please see the responses and comments on :

http://mathematica.stackexchange.com/questions/42835/simulated-annealing-convergence
POSTED BY: brama gatech
Answer
7 months ago