How can I compute the maximum value of a ConditionalExpression?

Question

If we use GeometricSolveValues in version 14.0, we can use this code to get a ConditionalExpression expr:

RandomInstance[
 sence = GeometricScene[{{P -> {0, 0}, X, Y, Z}, {p -> 3, 
     q -> 2}}, {circle = Circle[P, {p, q}], X ∈ circle, 
    Y ∈ circle, Z ∈ circle, 
    GeometricAssertion[Triangle[{X, Y, Z}], "Equilateral"]}], 
 RandomSeeding -> 1]
expr = GeometricSolveValues[sence, EuclideanDistance[X, Y]]

How can I compute the maximum value of the expr in MMA?

Answer 1

The InputForm of the output is

ConditionalExpression[EuclideanDistance[X, Y], (Indexed[X, {1}] - Indexed[Y, {1}])^2 + (Indexed[X, {2}] - 
   Indexed[Y, {2}])^2 == EuclideanDistance[X, Y]^2 && 
  EuclideanDistance[X, Y] >= 0 && p == 3 && q == 2 && 
  Indexed[P, {1}] == 0 && Indexed[P, {2}] == 0 && 
  p*Cos[C[1]] + Indexed[P, {1}] == Indexed[X, {1}] && 
  Indexed[P, {2}] + q*Sin[C[1]] == Indexed[X, {2}] && 
  p*Cos[C[2]] + Indexed[P, {1}] == Indexed[Y, {1}] && 
  Indexed[P, {2}] + q*Sin[C[2]] == Indexed[Y, {2}] && 
  p*Cos[C[3]] + Indexed[P, {1}] == Indexed[Z, {1}] && 
  Indexed[P, {2}] + q*Sin[C[3]] == Indexed[Z, {2}] && 
(Indexed[X, {1}] - 
   Indexed[Y, {1}])^2 + (Indexed[X, {2}] - 
   Indexed[Y, {2}])^2 == (Indexed[Y, {1}] - 
   Indexed[Z, {1}])^2 + (Indexed[Y, {2}] - 
   Indexed[Z, {2}])^2 == (Indexed[X, {1}] - 
   Indexed[Z, {1}])^2 + (Indexed[X, {2}] - Indexed[Z, {2}])^2 && 
  p > 0 && q > 0 && 0 <= C[1] <= 2*Pi && 0 <= C[2] <= 2*Pi &&  0 <= C[3] <= 2*Pi]

Making use of the canonical equation of the ellipse Circle[{0,0}, {3, 2}] instead of its parametric equation in the above, the optimization problem to find an equilateral inscribed triangle can be formulated in WL as follows (The notations x1, x2 ... are used instead of subscripts and (x1 - y1)^2 + (x2 - y2)^2 instead of EuclideanDistance.)

Maximize[{(x1 - y1)^2 + (x2 - y2)^2, 
x1^2/9 + x2^2/4 == 1 && y1^2/9 + y2^2/4 == 1 && 
z1^2/9 + z2^2/4 == 1 && (x1 - z1)^2 + (x2 - z2)^2 == 
(y1 - z1)^2 + (y2 - z2)^2 && (x1 - z1)^2 + (x2 - z2)^2 == 
(y1 - x1)^2 + (y2 -  x2)^2}, {x1, x2, y1, y2, z1, z2}]

The above code is running without any response for a long time. Applying NMaximize, one obtains

16.1831, {x1 -> 2.01347, x2 -> -1.48263, y1 -> -2.00936, y2 -> -1.48511, z1 -> -0.0000915262, z2 -> 2.}}

Edit. I stated the InputForm of the output in the question for the user's convenience. After that I solved the asked optimization problem and added that.