# What is the best optimization method for the following function?

I have calculated the energy of the electrons as a function of their angles. and I want to find the angles for which my function (energy) is minimized. I have used NMinimize, but the angles I get from it doesn't seem to be right as I expect a smooth change in the angles (My expectation may be wrong but I want to make sure that the problem is not coding or in other words the minimum is the Global minimum).

I have used Monte Carlo as well (I can post that if necessary but it doesn't find the global at all so I thought may be it's not the best way of optimization in this case). My first question is how to make sure the minimum energy I get is the Global minimum?

The second question is if you know any other better way to minimize my function so that I am sure the answer is the Global Min?

I post the code here (the angle for the first electron should be 0 and the last one should ne pi/2, so I fixed those in the begining):

 Clear["Global*"]
g[\[ScriptL]0_Integer?NonNegative] :=
g[\[ScriptL]0] =
Module[{ne = 2 \[ScriptL]0 + 2, vars, cons, \[Gamma], factorFxn,
ene}, vars = Array[Symbol["x$" <> ToString[#]] &, ne]; vars[[1]] = 0; vars[[ne]] = \[Pi]/2; cons = Apply[And, Thread[0 <= vars[[2 ;; ne - 1]] <= \[Pi]/2]]; \[Gamma] = Join[Table[{Riffle[Range[0, -\[ScriptL]0, -1], Range[\[ScriptL]0]][[i]], 1}, {i, 1, 2 \[ScriptL]0 + 1}], Table[{Riffle[Range[0, -\[ScriptL]0, -1], Range[\[ScriptL]0]][[i]], -1}, {i, 1, 2 \[ScriptL]0 + 1}]]; factorFxn[\[ScriptL]_, m1_, m2_, p1_, p2_] := (If[\[Gamma][[m1, 2]] == \[Gamma][[p1, 2]] && \[Gamma][[m2, 2]] == \[Gamma][[p2, 2]], If[\[Gamma][[p1, 1]] - \[Gamma][[m1, 1]] == \[Gamma][[m2, 1]] - \[Gamma][[p2, 1]], Sum[(2 \[ScriptL] + 1)^2 Sum[ If[\[Gamma][[p1, 1]] - \[Gamma][[m1, 1]] == mval && \[Gamma][[m2, 1]] - \[Gamma][[p2, 1]] == mval, (-1)^(\[Gamma][[m1, 1]] + \[Gamma][[m2, 1]] + mval) ThreeJSymbol[{\[ScriptL], -\[Gamma][[m1, 1]]}, {\[ScriptL], \[Gamma][[p1, 1]]}, {\[ScriptL]temp, -mval}] ThreeJSymbol[{\ \[ScriptL]temp, mval}, {\[ScriptL], -\[Gamma][[m2, 1]]}, {\[ScriptL], \[Gamma][[p2, 1]]}] ThreeJSymbol[{\[ScriptL], 0}, {\[ScriptL], 0}, {\[ScriptL]temp, 0}]^2, 0], {mval, -\[ScriptL]temp, \[ScriptL]temp}], \ {\[ScriptL]temp, 0, 2 \[ScriptL]}], 0] , 0] - If[\[Gamma][[m1, 2]] == \[Gamma][[p2, 2]] && \[Gamma][[m2, 2]] == \[Gamma][[p1, 2]], If[-\[Gamma][[p1, 1]] + \[Gamma][[m2, 1]] == -\[Gamma][[m1, 1]] + \[Gamma][[p2, 1]], Sum[(2 \[ScriptL] + 1)^2 Sum[ If[-\[Gamma][[p1, 1]] + \[Gamma][[m2, 1]] == mval && -\[Gamma][[m1, 1]] + \[Gamma][[p2, 1]] == mval, (-1)^(\[Gamma][[m1, 1]] + \[Gamma][[m2, 1]] + mval) ThreeJSymbol[{\[ScriptL], -\[Gamma][[m1, 1]]}, {\[ScriptL], \[Gamma][[p2, 1]]}, {\[ScriptL]temp, -mval}] ThreeJSymbol[{\ \[ScriptL]temp, mval}, {\[ScriptL], -\[Gamma][[m2, 1]]}, {\[ScriptL], \[Gamma][[p1, 1]]}] ThreeJSymbol[{\[ScriptL], 0}, {\[ScriptL], 0}, {\[ScriptL]temp, 0}]^2, 0], {mval, -\[ScriptL]temp, \[ScriptL]temp}], \ {\[ScriptL]temp, 0, 2 \[ScriptL]}], 0] , 0] ); ene[th1_?VectorQ] :=(*(2 \[ScriptL]0 +1)^2*) Sum[ (* Find out which states we're calculating the matrix element of \ *) (Cos[th1[[m2]]] Cos[th1[[m1]]] + Cos[th1[[m2]]] Sin[th1[[m1 + 1]]] + Cos[th1[[m1]]] Sin[th1[[m2 + 1]]] + Sin[th1[[m2 + 1]]] Sin[th1[[m1 + 1]]] + If[m1 == m2, Cos[th1[[m1]]] Sin[th1[[m2 + 1]]], 0])*(Cos[th1[[p2]]] Cos[th1[[p1]]] + Cos[th1[[p2]]] Sin[th1[[p1 + 1]]] + Cos[th1[[p1]]] Sin[th1[[p2 + 1]]] + Sin[th1[[p2 + 1]]] Sin[th1[[p1 + 1]]] + If[p1 == p2, Cos[th1[[p1]]] Sin[th1[[p1 + 1]]], 0]) factorFxn[\[ScriptL]0, m1, m2, p1, p2] , {p1, 1, 2 \[ScriptL]0 + 1}, {m1, 1, 2 \[ScriptL]0 + 1}, {p2, 1, p1}, {m2, 1, m1}]; {ne, #} & /@ NMinimize[{ene[vars], cons}, vars[[2 ;; ne - 1]], Method -> "SimulatedAnnealing"]]  I would appreciate any suggestions to go about this. • just for you or anyone looking to answer this, I think a strategy would be to look at the derivative of your function over the domain and have mathematica Reduce to find the regions over which the derivative is positive and negative. If you get a result, this can be used to ensure that you know where all possible minima are, and then you can check each one. You may have to do some more math to ensure you have a global minimum if mathematica can't Reduce the appropriate inequalities... May 21 at 22:57 • For comparison: With Method -> "DifferentialEvolution" g[11]//AbsoluteTiming produces {924.257, {{24, 302.918}, {24, {x$2 -> 0., x$3 -> 0., x$4 -> 0., x$5 -> 0., x$6 -> 0., x$7 -> 0., x$8 -> 0., x$9 -> 0., x$10 -> 0., x$11 -> 0., x$12 -> 0., x$13 -> 1.5708, x$14 -> 1.5708, x$15 -> 1.5708, x$16 -> 1.5708, x$17 -> 1.5708, x$18 -> 1.5708, x$19 -> 1.5708, x$20 -> 1.5708, x$21 -> 1.5708, x$22 -> 1.5708, x$23 -> 1.5708}}}} May 22 at 5:24 • With Method -> "SimulatedAnnealing" g[11]//AbsoluteTiming produces {677.309,{{24,393.869},{24,{x$2->1.5708,x$3->1.5708,x$4->1.5708,x$5->1.5708,x$6->1.5708,x$7->1.5708,x$8->1.5708,x$9->0.,x$10->0.,x$11->0.,x$12->0.,x$13->1.5708,x$14->1.5708,x$15->1.5708,x$16->1.57077,x$17->0.0000281708,x$18->0.,x$19->0.,x$20->0.,x$21->0.,x$22->1.5708,x$23->1.5708}}}}. May 22 at 5:26 • With Method -> "NelderMead" g[11]//AbsoluteTiming produces {672.662, {{24, 347.683}, {24, {x$2 -> 0., x$3 -> 1.5708, x$4 -> 1.5708, x$5 -> 0., x$6 -> 0., x$7 -> 0., x$8 -> 0., x$9 -> 1.5708, x$10 -> 1.5708, x$11 -> 1.5708, x$12 -> 1.5708, x$13 -> 1.5708, x$14 -> 1.5708, x$15 -> 1.5708, x$16 -> 1.5708, x$17 -> 1.5708, x$18 -> 1.5708, x$19 -> 1.5708, x$20 -> 1.5708, x$21 -> 1.5708, x$22 -> 1.5708, x\$23 -> 1.5708}}}}. Method->"RandomSearch"` is behind. May 22 at 5:30
• @user64494 Thanks a lot, I have already tried those, Do you think a MonteCarlo routine would be more effective in some cases than the built-in Mathematica optimization functions? May 22 at 19:05