Maximize does not return the global maximum. Why?

Question

Bug introduced in between 10.0.1 and 10.4.1, and fixed in 11.1

---

I have the following code

Maximize[{-(((1-3 p0^2-3 p1^2-3 p2^2) (18 p0 p1-27 p0^2 p1-27 p0 p1^2+18 p0    p2-27 p0^2 p2+18 p1 p2-27 p1^2 p2-27 p0 p2^2-27 p1 p2^2))/(9 Sqrt[3]))+((-2+9 p0^2-9 p0^3+9 p1^2-9 p1^3+9 p2^2-9 p2^3) (-p1 p2-p0 (p1+p2)))/(3 Sqrt[3]),p0>=0&&p0<=1&&p1>=0&&p1<=1&&p2>=0&&p2<=1&&p0+p1+p2==1&&1/27 (-2+9 p0^2-9 p0^3+9 p1^2-9 p1^3+9 p2^2-9 p2^3)<=0},{p0,p1,p2}]

The result I get is zero. But this is wrong. This triple

p0->0.284088318172241,p1->0.4318233636555181,p2->0.284088318172241

satisfies all the constraints and the objective expression is positive. It equals 0.0000921367. Curiously, if I multiply the objective function by $\sqrt{3}$ I get the presumed global maximum

$-{16\over(-50013-6049\sqrt{69})}>0.$

The documentation to Maximize says If f and cons are linear or polynomial, Maximize will always find a global maximum. But $\sqrt{3}$ is just a number in my expression and it should not matter that it is not rational.

Answer 1

Amplifying on comment by Michael E2

$Version

(*  "11.1.0 for Mac OS X x86 (64-bit) (March 16, 2017)"  *)

sys = {-(((1 - 3 p0^2 - 3 p1^2 - 3 p2^2) (18 p0 p1 - 27 p0^2 p1 - 
           27 p0 p1^2 + 18 p0 p2 - 27 p0^2 p2 + 18 p1 p2 - 27 p1^2 p2 - 
           27 p0 p2^2 - 27 p1 p2^2))/(9 Sqrt[3])) + ((-2 + 9 p0^2 - 9 p0^3 + 
         9 p1^2 - 9 p1^3 + 9 p2^2 - 9 p2^3) (-p1 p2 - p0 (p1 + p2)))/(3 Sqrt[
        3]), p0 >= 0 && p0 <= 1 && p1 >= 0 && p1 <= 1 && p2 >= 0 && p2 <= 1 &&
     p0 + p1 + p2 == 1 && 
    1/27 (-2 + 9 p0^2 - 9 p0^3 + 9 p1^2 - 9 p1^3 + 9 p2^2 - 9 p2^3) <= 0};

The exact solution is

max1 = Maximize[sys, {p0, p1, p2}]

(*  {-((22228 Sqrt[3])/1953125) + (24196 Sqrt[23])/
  5859375, {p0 -> 
   1 - Root[{-23 + #1^2 &, -3 + #2^2 &, 
      2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 + 
        7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 - 
        71191406250 #2 #3^5 + 35595703125 #2 #3^6 &}, {2, 2, 2}] - 
    Root[{-23 + #1^2 &, -3 + #2^2 &, 
      2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 + 
        7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 - 
        71191406250 #2 #3^5 + 35595703125 #2 #3^6 &, 
      72588 #1 - 200052 #2 + 7812500 #2 #3 - 25390625 #2 #3^2 + 
        35156250 #2 #3^3 - 17578125 #2 #3^4 + 7812500 #2 #4 - 
        148437500 #2 #3 #4 + 439453125 #2 #3^2 #4 - 562500000 #2 #3^3 #4 + 
        263671875 #2 #3^4 #4 - 25390625 #2 #4^2 + 439453125 #2 #3 #4^2 - 
        843750000 #2 #3^2 #4^2 + 527343750 #2 #3^3 #4^2 + 35156250 #2 #4^3 - 
        562500000 #2 #3 #4^3 + 527343750 #2 #3^2 #4^3 - 17578125 #2 #4^4 + 
        263671875 #2 #3 #4^4 &}, {2, 2, 2, 1}], 
  p1 -> Root[{-23 + #1^2 &, -3 + #2^2 &, 
     2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 + 
       7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 - 
       71191406250 #2 #3^5 + 35595703125 #2 #3^6 &}, {2, 2, 2}], 
  p2 -> Root[{-23 + #1^2 &, -3 + #2^2 &, 
     2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 + 
       7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 - 
       71191406250 #2 #3^5 + 35595703125 #2 #3^6 &, 
     72588 #1 - 200052 #2 + 7812500 #2 #3 - 25390625 #2 #3^2 + 
       35156250 #2 #3^3 - 17578125 #2 #3^4 + 7812500 #2 #4 - 
       148437500 #2 #3 #4 + 439453125 #2 #3^2 #4 - 562500000 #2 #3^3 #4 + 
       263671875 #2 #3^4 #4 - 25390625 #2 #4^2 + 439453125 #2 #3 #4^2 - 
       843750000 #2 #3^2 #4^2 + 527343750 #2 #3^3 #4^2 + 35156250 #2 #4^3 - 
       562500000 #2 #3 #4^3 + 527343750 #2 #3^2 #4^3 - 17578125 #2 #4^4 + 
       263671875 #2 #3 #4^4 &}, {2, 2, 2, 1}]}}  *)

Simplifying the forms of p0, p1, and p2

max2 = max1 // RootReduce // ToRadicals

(*  {8/3 Sqrt[2/(837672973 + 100842879 Sqrt[69])], {p0 -> 1/75 (49 - 2 Sqrt[69]), 
  p1 -> 1/75 (13 + Sqrt[69]), p2 -> 1/75 (13 + Sqrt[69])}}  *)

Verifying that all of the conditions are met

verify = sys /. max2[[2]] // Simplify

(*  {(4 (-50013 + 6049 Sqrt[69]))/(5859375 Sqrt[3]), True}  *)

The form of the maximum given in verify is simpler than the ones in max1 or max2

LeafCount /@ First /@ {max1, max2, verify}

(*  {19, 21, 18}  *)

The different forms of the maximum are equal

max1[[1]] == max2[[1]] == verify[[1]] // FullSimplify

(*  True  *)

Answer 2

Mathematica does not give you the whole truth.

Maximize does not recognize, that the system is totaly symmetric and there are three equal maxima,

corresponding to the six permutations of {p0,p1,p2} divided by two, because of two p's are equal (I corrected my first version).

Redefine sys, calculate the maximum and take the permutations (using MMA 8.0)

    sys[p0_, p1_, p2_] = sys

    max1 = Maximize[sys[p0, p1, p2], {p0, p1, p2}]

    (*   {-((22228 Sqrt[3])/1953125) + (24196 Sqrt[23])/
         5859375, {p0 -> 1 + 2/75 (-13 - Sqrt[69]), 
          p1 -> 1/75 (13 + Sqrt[69]), p2 -> 1/75 (13 + Sqrt[69])}}    *)

    sys[Sequence @@ #] & /@ Permutations[{p0, p1, p2}] /. 
        max1[[2]] // N

    (*   {{0.0000921367, True}, {0.0000921367, True}, {0.0000921367, 
          True}, {0.0000921367, True}, {0.0000921367, True}, {0.0000921367, 
          True}}    *)