NMaximize - global maximum fail - How to?

Question

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[1]:= 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[1]= {0.226147,{x->0.,y->0.455858}}
(*Failed*)

In[2]:= 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[2]= {0.290656,{x->0.754955,y->0.,z->0.245045}}
(*Wow, it works for 3 variables! The original problem*)

In[3]:= 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[3]= {0.290656,{x->0.754974,y->0.}}
(*Changing 0 to 0.01 works! *)

In[4]:= 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[4]= {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 :

Answer 1

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[10]), {x -> 1 + 1/15 (-10 + 2 Sqrt[10]), y -> 0, z -> 1/15 (10 - 2 Sqrt[10])}}

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[10]), {x -> 1/15 (5 + 2 Sqrt[10]), 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.}}

Answer 2

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[10], f[11], f[20]}, {f[29], f[30], f[50]}}]