# NMaximize - global maximum fail - How to?

Problem :

Let x,y,z be non-negative real number whose sum is 1.

Find maximum of 0.2x + 0.2y + x^2*z + y*z^2


What I did first was reducing variable.

0.2x + 0.2y + x^2(1-x-y) + y(1-x-y)^2 with constraint
x>=0, y>=0, x+y<=1


Mathematica code is

NMaximize[{0.2x+0.2y+x^2(1-x-y)+y(1-x-y)^2,0<=x<=1&&0<=y<=1&&x+y<=1},{x,y}]


and I got

{0.226147,{x->0.,y->0.455858}}


But it is not a global maximum.

In fact, global maximum is

{0.290656, {x->0.754955, y->0.}}


I know that NMaximize tries to give us global maximum, but it can fail. And if fail, mathematica gives us local maximum.
But, I didn't expect NMaximize could fail to find the global maximum for such a easy problem.

Q1) Is it possible to output an additional message:

{0.226147,{x->0.,y->0.455858}}
Bad luck for NMaximize : Failed to ensure this is the global maximum


Q2) Is there a sure way even if it takes more time ? (finding global maximum)

Below is a helpful code... but they are like seeing the results first and fitting things

In:= NMaximize[{0.2x+0.2y+x^2(1-x-y)+y (1-x-y)^2,0<=x<=1&&0<=y<=1&&x+y<=1},{x,y}]
Out= {0.226147,{x->0.,y->0.455858}}
(*Failed*)

In:= NMaximize[{0.2x+0.2y+x^2z+y z^2,0<=x<=1&&0<=y<=1&&0<=z<=1&&x+y+z==1},{x,y,z}]
Out= {0.290656,{x->0.754955,y->0.,z->0.245045}}
(*Wow, it works for 3 variables! The original problem*)

In:= NMaximize[{0.2x+0.2y+x^2(1-x-y)+y (1-x-y)^2,0.01<=x<=1&&0<=y<=1&&x+y<=1},{x,y}]
Out= {0.290656,{x->0.754974,y->0.}}
(*Changing 0 to 0.01 works! *)

In:= NMaximize[{0.2x+0.2y+x^2(1-x-y)+y (1-x-y)^2,0.001<=x<=1&&0<=y<=1&&x+y<=1},{x,y}]
Out= {0.225852,{x->0.001,y->0.455708}}
(*0.001 not works*)


======================================= After a getting a comment from from Akku14 :

Thank you, WorkingPrecision works, but surprisingly it is
not true
that "the higher, the better". See below. I recommend to try with some range of numbers (for WorkingPrecision) Below is a photo magnified : • Often higher WorkingPrecision helps. (but if higher than \$MachinePrecision, you have to Rationalize your function first). NMaximize[{0.2 x + 0.2 y + x^2 (1 - x - y) + y (1 - x - y)^2, 0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1}, {x, y}, WorkingPrecision -> 15]  yields {0.290656125830494, {x -> 0.754985264404849, y -> 0}}  Oct 24, 2021 at 15:21
• Thank you! I added some surprising phenomenon in the question article. Oct 24, 2021 at 22:48
• Plot3D[ConditionalExpression[0.2 x + 0.2 y + x^2 (1 - x - y) + y (1 - x - y)^2, 0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1], {x, 0, 1}, {y, 0, 1}] show why it's reasonable that NMaximize could find either point. Also, it suggests that WorkingPrecision is irrelevant. More relevant are the random points it starts with, as well as the method and the method's tuning parameters. Compare Method -> {"NelderMead", "RandomSeed" -> 1} with the method used by default here, Method -> {"NelderMead"}.... Oct 24, 2021 at 23:53
• Comparing SeedRandom; RandomReal[1, 3] and SeedRandom; RandomReal[1, 3, WorkingPrecision -> 15] shows that WorkingPrecision affects random number generation, which accounts for it affecting the result. Oct 24, 2021 at 23:55
• Thank you very much ! The technique Method -> {"NelderMead", "RandomSeed" -> WithANumber} is just invaluable. One of the most important tip I have ever learned! Oct 25, 2021 at 7:01

Maximize[{x/5 + y/5 + x^2*z + y*z^2, x + y + z == 1&& x >= 0&& y >= 0&& z >= 0},{x,y,z}]


{1/675 (95 + 32 Sqrt), {x -> 1 + 1/15 (-10 + 2 Sqrt), y -> 0, z -> 1/15 (10 - 2 Sqrt)}}

N[%]


{0.290656, {x -> 0.75497, y -> 0., z -> 0.24503}}

Maximize[{x/5 + y/5 + x^2 (1 - x - y) + y (1 - x - y)^2, 0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1}, {x, y}]


{1/675 (95 + 32 Sqrt), {x -> 1/15 (5 + 2 Sqrt), y -> 0}}

NMaximize[{x/5 + y/5 + x^2 (1 - x - y) + y (1 - x - y)^2, 0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1}, {x, y}, Method -> "DifferentialEvolution"]


{0.290656, {x -> 0.75497, y -> 0.}}

This is just an extended comment to show what is happening in terms of the iterative nature as opposed to "why" it's doing what it's doing and that it's always helpful to plot one's function when one can.

Below are the iterative processes for WorkingPrecision values of 10, 11, 20, 29, 30, and 50. The green dots are the starting values (all the same), the red dots are the ending values, and the black line segments show the path of the iterative process.

Note that the starting value is pretty close to a saddle point. That probably spells trouble for iterative methods.

f[wp_] := Module[{a}, a = Reap[NMaximize[{x/5 + y/5 + x^2 (1 - x - y) + y (1 - x - y)^2,
0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1}, {x, y}, WorkingPrecision -> wp,
EvaluationMonitor :> Sow[{x, y}]]];
Show[ContourPlot[x/5 + y/5 + x^2 (1 - x - y) + y (1 - x - y)^2, {x, 0, 1}, {y, 0, 1 - x},
PlotLabel -> "WorkingPrecision -> " <> ToString[wp],
PlotPoints -> 200, Contours -> 100, PlotRange -> All, ImageSize -> Medium],
ListPlot[{a[[2, 1]], {a[[2, 1, 1]]}, {a[[2, 1, Length[a[[2, 1]]]]]}},
PlotStyle -> {Black, {Green, PointSize[0.03]}, { Red, PointSize[0.03]}, Black},
Joined -> {True, False, False}, PlotRange -> {{0, 1}, {0, 1}},
PlotRangeClipping -> False]]]

GraphicsGrid[{{f, f, f}, {f, f, f}}] 