Improve Drake support of SOS programming over Ideals

[AlexandreAmice]

Is your feature request related to a problem? Please describe. Many programs for dynamical system naturally have constraints of the form "f(x) is SOS for all g(x) = 0". For example, when searching for a Lyapunov function for pendulum dynamics, we require a Lyapunov function that is SOS for all s^2 + c^2 = 1.

Currently, the way users handle this is via the multiplier method which very quickly becomes intractable due to the bloat of the program. The correct way to handle these constraints is to operate on the coordinate ring of the variety associated to the equality constraint. In general, this is challenging but Drake should make it easier for us to do this for at least commonly encountered varieties.

Describe the solution you'd like The two known ways to operate on the coordinate ring is to either consider the Groebner basis of the associated variety or by sampling the coordinate ring. We should implement both approaches as sometimes a Groebner basis is easy to obtain, while sampling from the coordinate ring is hard, while in other instances sampling from the coordinate ring is easy, but a Groebner basis will take exponentially long to compute.

For both approaches, we should implement the class PolynomialIdeal

• [ ] A polynomial ideal is internally just a set of generator polynomials.
• [ ] Add to MathematicalProgram the function IsSosModuloIdeal(p, ideal, options) which will add the constraint that p is SOS via the multiplier method. Options should include things such as the maximum degree of the multiplier polynomials.

For the Groebner Basis approach:

• [ ] Subclass from PolynomialIdeal to make the GroebnerBasisIdeal
• [ ] Implement is_groebner_basis() which will check whether a generating set is a Groebner Basis. We should not (at least immediately) implement the construction of a Groebner basis from an arbitrary generating set.
• [ ] Implement normal_form(p) which given a polynomial p will compute its normal form modulo the GroebnerBasis
• [ ] Implement GetStandardMonomialBasis(degree) which returns a MonomialBasis of all the standard monomials in the ideal up to a certain degree.
• [ ] Modify IsSosModuloIdeal(p, ideal, options) to add the SOS constraint modulo the GroebnerBasisIdeal. A reference for GroebnerBasis algorithms is "Becker, T., Weispfenning, V., Becker, T., & Weispfenning, V. (1993). Gröbner bases (pp. 187-242). Springer New York."

For the Sampling approach:

• [ ] Implement SamplingVariety. Formally, this will be a collection of polynomials along with a function that is capable of sampling from the variety. I believe the best design pattern here is SamplingVariety has-a PolynomialIdeal for representing the polynomial collection and has-a sampling function to do the sampling. While ideals and varieties are very related they are NOT formally the same thing.
• [ ] Implement is_poised_sample_set. See sections 3 and 5 of this paper.
• [ ] Add IsSosOnVariety(p, sampling_variety) which would add the necessary constraints.

To achieve these goals, we likely want to close out or make progress on #13803 and #13602

[hongkai-dai]

The plan sounds good. Who is the immediate user of these approaches? Do we have a project that would benefit from the SOS over Ideals implementation? Answering those would help decide the priority of the task.

[AlexandreAmice]

There have been various efforts in our lab thinking about doing SOS over equality constraints. See Lu's work, Shen's work, etc. I may be doing more in this direction soon. I believe the priority can remain low for now. I just wanted to put out some design ideas before starting slowly.