PointInSet for VPolytope with tolerance=0 gives unexpected results

[RussTedrake]

This has long been discussed, but now may be a blocker for #17167, so I'm raising it to an issue.

Paraphrasing a discussion I've had with @umenberger :

I really need a better solution for PointInSet for VPolytopes. The one I have looks strange because I was trying hard to make the tolerance play well with the solver tolerances. But it has the quirky property that right now if I have points that are safely in the interior of the polytope, and ask for a certification of that with tolerance 0, it will often fail.

The implementation is here, with a small pre-amble here. We need to decide on a better solution.

cc @TobiaMarcucci .

[RussTedrake]

For a victory condition: we should be able to set many of the tolerances in vpolytope_test.cc to zero. For instance this line that tests that the center of a triangle is in the triangle.

[TobiaMarcucci]

I think that checking if a point is in the convex hull of other points is necessarily an LP: I didn't prove it, but I'd imagine that any LP can be written in this form. So I think calling an LP solver is unavoidable, as well as having to deal with its tolerances (otherwise we'd have an LP solver with infinite precision). Certain LP algorithms are more precise than others (simplex gives almost machine precision, IP tolerance is user defined). But I guess we don't want to go in these weeds.

If we think that our check can never be more precise than the LP solver we use, I imagine that it'd be fine to keep the check we have but instead of minimizing the norm of the residual x - sum_i a_i v_i just enforce x = sum_i a_i v_i explicitly. If the problem is infeasible, then x is not in the convex hull of v1, ..., vn. If we do want to also have a user-defined tolerance for this check than I'd do -1*tol <= x - sum_i a_i v_i <= 1*tol or something like that?

I hope I got your point and I didn't say anything too obvious.

[AlexandreAmice]

Checking membership in a VPolytope indeed requires solving an LP, or requires converting it to an HPolyhedron using Fourier-Motzkin (if you want to avoid numerical precision issues). Two potential areas of improvement: 1) Use a relative tolerance rather than an absolute tolerance. 2) Perhaps we should add to MathematicalProgram a way to query each solver's constraint tolerance. Then when the tolerance is set to 0, we default to this value.

[cohnt]

If the simplex method will give higher precision, we could query the solver, and set the appropriate solve option to force it to use the simplex method. I know Gurobi and Mosek have the option to specify simplex vs interior point. It looks like CLP does as well, but do we have a way of setting those options in SolverOptions?

[jwnimmer-tri]

Aside: the SolverOptions is a mapping from SolverId to its associated options, so you don't necessarily need to "query the solver" -- you can blindly set options for Gurobi, Mosek, CLP, etc. all simultaneously. Whichever solver gets chosen for the Solve will just use its associated options and ignore the other solvers.