NMaximize fails for linear program?

Question

Consider the following code, maximizing the variable diff subject to some constraints. Notice that all constraints are linear (in)equalities -- the integrality constraints are commented out. (Apologies for the length of the code, but removing any constraints no longer allows reproducing the issue.)

Clear["Global`*"]

n = 9;
r = 6;

(*allow writing e.g. d[2,1] to refer to variable d12*)
d[i_, j_] := ToExpression["d" ~~ ToString@Row@Sort[{i, j}]];
w[i_, j_, l_] := ToExpression["w" ~~ ToString@Row@Sort[{i, j, l}]];

vars := Flatten@Join[
    Table[d[i, j], {i, n - 1}, {j, i + 1, n}],
    Table[w[i, j, l], {i, n - 2}, {j, i + 1, n - 1}, {l, j + 1, n}],
    Table[m3[i, j],
     {i, DeleteCases[Range[1, n], r]},
     {j, DeleteCases[Range[1, n], r | i]}],
    Table[m2[i], {i, DeleteCases[Range[1, n], r]}],
    {m1, phi, phiAlt, diam, diff} 
    ];

NMaximize[{diff,
  (*Table[Element[d[i,j],Integers], {i,n-1},{j,i+1,n}],*)
  (*Table[Element[w[i,j,l],Integers],{i,n-2},{j,i+1,n-1},{l,j+1,n}],*)
  Table[d[i, j] <= d[i, l] + d[j, l] &&
    d[i, l] <= d[i, j] + d[j, l] &&
    d[j, l] <= d[i, j] + d[i, l],
   {i, 1, n - 2}, {j, i + 1, n - 1}, {l, j + 1, n}],
  Table[
   -d[i, j] <= w[i, l, m] - w[j, l, m] <= d[i, j] &&
    -d[i, l] <= w[i, j, m] - w[l, j, m] <= d[i, l] &&
    -d[i, m] <= w[i, j, l] - w[m, j, l] <= d[i, m] &&
    -d[j, l] <= w[j, i, m] - w[l, i, m] <= d[j, l] &&
    -d[j, m] <= w[j, i, l] - w[m, i, l] <= d[j, m] &&
    -d[l, m] <= w[l, i, j] - w[m, i, j] <= d[l, m],
   {i, 1, n - 3},
   {j, i + 1, n - 2},
   {l, j + 1, n - 1},
   {m, l + 1, n}
   ],
  Table[m3[p3, p2] <= w[p2, j, r] - d[j, p3],
   {p3, DeleteCases[Range[1, n], r]},
   {p2, DeleteCases[Range[1, n], r | p3]},
   {j, DeleteCases[Range[1, n], r | p3 | p2]}
   ],
  Table[m2[p2] <= w[j, l, r] - d[j, p2] - d[l, p2],
   {j, DeleteCases[Range[1, n - 1], r]},
   {l, DeleteCases[Range[j + 1, n], r]},
   {p2, DeleteCases[Range[1, n], j | l | r]}
   ],
  Table[m1 <= w[j, l, r] - d[j, r] - d[l, r],
   {j, DeleteCases[Range[1, n - 1], r]},
   {l, DeleteCases[Range[j + 1, n], r]}
   ],
  Table[phiAlt <= w[p3, p2, r] + m3[p3, p2] + m2[p2] + m1,
   {p3, DeleteCases[Range[1, n], r]},
   {p2, DeleteCases[Range[1, n], r | p3]}
   ],
  Table[w[i, j, l] >= 0, {i, n - 2}, {j, i + 1, n - 1}, {l, j + 1, n}],
  Table[0 <= d[i, j] <= diam, {i, n - 1}, {j, i + 1, n}],
  w128 == w126 + d68,
  w349 == w346 + d69,
  w578 == w567 + d68,
  w589 == w568 + d69,
  w789 == w679 + d68,
  d34 + d89 <= d38 + d49 == d39 + d48,
  d27 + d59 <= d29 + d57 == d25 + d79,
  d27 + d69 <= d29 + d67 == d26 + d79,
  d27 + d56 <= d25 + d67 == d26 + d57,
  w156 + w268 >= w168 + w256 == w126 + w568,
  w267 + w569 >= w269 + w567 == w256 + w679,
  w569 + w678 >= w567 + w689 == w568 + w679,
  phi == w789 + w589 - d57 + w349 - d38 - d48 + w126 - d19 - d29 - d69,
  diff <= phiAlt - phi,
  diff <= 1,
  diam <= 3
  }, vars]

If I run this on Mathematica 12.0 on MacOS, it outputs "There are no points that satisfy the constraints" and "The following constraints are not valid [...]". If I now add the constraints

Table[Element[d[i,j],Integers], {i,n-1},{j,i+1,n}]

(commented out above), turning the linear program into a mixed integer linear program, Mathematica now finds a feasible solution with value 0.333333. If I further add the constraints

Table[Element[w[i,j,l],Integers],{i,n-2},{j,i+1,n-1},{l,j+1,n}]

it finds an even better solution with value 1. Clearly, adding constraints should neither turn an infeasible problem into a feasible one nor improve the objective value. I speculate that the observed infeasibility might be due to numerical errors and the 0.333333 might be a local rather than global maximum. However, the documentation of NMaximize states that "If f and cons are linear, NMaximize can always find global maxima, over both real and integer values", which seems incorrect in light of this example. Can anyone shed light onto what is going on here? How can I know whether I can trust the output of Mathematica, or how can I avoid getting incorrect results on purely linear programs in the first place?

Note: If instead of Mathematica 12.0 on a Mac I run it on the old Mathematica 11.3 on a Linux system, it outputs objective value 1 in all cases.

Answer 1

In the new Mathematica 12.1, the above code now runs correctly, so I am hopeful that the bug has been fixed.