Why is FindInstance failing when I relax a set of constraints?

Question

I'm attempting to use FindInstance to generate coordinate sets for plausible triangles with edge length distance constraints. E.g.:

PlaneDimensionXYnm = 1000;
TriangleEdgeDistance = 400;

TriangleSolution = {};
While[Length[TriangleSolution] == 0,
  x1 = RandomReal[{200, PlaneDimensionXYnm - 200}];
  y1 = RandomReal[{200, PlaneDimensionXYnm - 200}];

TriangleSolution = FindInstance[  
    EuclideanDistance[{x1, y1}, {x2, y2}] == TriangleEdgeDistance &&  
    EuclideanDistance[{x1, y1}, {x3, y3}] == TriangleEdgeDistance &&  
    EuclideanDistance[{x2, y2}, {x3, y3}] == TriangleEdgeDistance &&  
    x2 >= 200 && y2 >= 200 &&  
    x3 >= 200 &&  y3 >= 200 &&  
    x2 <= PlaneDimensionXYnm - 200 &&  
    y2 <= PlaneDimensionXYnm - 200 &&  
    x3 <= PlaneDimensionXYnm - 200 &&  
    y3 <= PlaneDimensionXYnm - 200, {x2, y2, x3, y3}];
    ];

TriangleCoordinates = N[{x2, y2, x3, y3} /. TriangleSolution[[1]]];

While the above script works fine (albeit inefficiently, but it doesn't matter for me), the below script fails:

TriangleSolution = {};
TriangleEdgeDistance = 400;

While[Length[TriangleSolution] == 0,
  x1 = RandomReal[{200, PlaneDimensionXYnm - 200}];
  y1 = RandomReal[{200, PlaneDimensionXYnm - 200}];  

TriangleSolution = FindInstance[  
    EuclideanDistance[{x1, y1}, {x2, y2}] > TriangleEdgeDistance &&  
    EuclideanDistance[{x1, y1}, {x3, y3}] > TriangleEdgeDistance &&  
    EuclideanDistance[{x2, y2}, {x3, y3}] > TriangleEdgeDistance &&  
    x2 >= 200 && y2 >= 200 &&  
    x3 >= 200 && y3 >= 200 &&  
    x2 <= PlaneDimensionXYnm - 200 &&  
    y2 <= PlaneDimensionXYnm - 200 &&  
    x3 <= PlaneDimensionXYnm - 200 &&  
    y3 <= PlaneDimensionXYnm - 200, {x2, y2, x3, y3}];
    ];

TriangleCoordinates = N[{x2, y2, x3, y3} /. TriangleSolution[[1]]];

All I've done here is relax the constraint that the triangle edges should be the same length, and required that the edges all be greater than a certain length. Why does FindInstance now fail?

Answer 1

For nonlinear systems of equations and inequalities FindInstance uses the cylindrical algebraic decomposition (CAD) algorithm. The algorithm may use equational constraints to simplify the computations by "eliminating" variables. Hence replacing equations with inequalities can make the problem harder.

Another problem in this example is that EuclideanDistance returns a quadratic radical. The CAD algorithm deals only with polynomial systems. FindInstance polynomializes the input by replacing radicals with new variables and adding the corresponding quadratic equations. The CAD algorithm is very sensitive to adding new variables - its worst-case complexity in the number of variables is doubly-exponential. Since EuclideanDistance is a square root of an expression that is always nonnegative, in this example the quadratic radicals may be removed without adding new variables, namely by squaring both sides of each equation or inequality. If you use

EuclideanDistance[{x1, y1}, {x2, y2}]^2 == TriangleEdgeDistance^2

EuclideanDistance[{x1, y1}, {x2, y2}]^2 > TriangleEdgeDistance^2

in you code the example runs twice faster in the equation case and ~100 times faster in the inequality case.