mathematicalprogram: add constraint naming/debugging workflow to tutorials

[RussTedrake]

I just responded to a student question, and realized that I might have been hiding a pretty useful debugging workflow. @hongkai-dai -- would you consider adding this to your MP tutorials?

But please know that there are some good tools for debugging constraints in mathematical program.  > I'm sorry we haven't advertised them better.

Take a look at the trick shot example in https://github.com/RussTedrake/underactuated/blob/master/examples/contact.ipynb e.g.

prog.AddConstraint(xdot[1] == -e * xdot[0]).evaluator().set_description('xdot_wall')
prog.AddConstraint(zdot[1] == 3.*zdot_pre/5. - 2.*r*thetadot[0]/5.).evaluator().set_description('zdot_wall')

You can add names for each constraints with that (slightly clunky) API.   then later, you can do something like:

result = Solve(prog)
if not result.is_success():
infeasible = GetInfeasibleConstraints(prog, result)
print("Infeasible constraints:")
for i in range(len(infeasible)):
print(infeasible[i])

to get a super useful display that tells you which constraints are failing.

[hongkai-dai]

Good call, I will add it to the tutorial.

Often I find that I need not only the name of the infeasible constraints, but also the Binding. I am considering to provide two APIs

1. Rename GetInfeasibleConstraints to GetInfeasibleConstraintNames which returns a std::vector containing the name of the constraint.
2. GetInfeasibleConstraints returns std::vector<Binding> instead of std::vector

[hongkai-dai]

@jwnimmer-tri may I ask you what would be a good workflow to change this API? Currently GetInfeasibleConstraints returns std::vector<string> as the names of the infeasible constraints. I would like this function to return std::vector<Binding<Constraint>>, and add a new function GetInfeasibleConstraintNames to return std::vector<string>. I imagine the workflow could be like this

1. Add a new function std::vector<string> GetInfeasibleConstraintNames(...). Add DRAKE_DEPRECATE to the existing std::vector<string> GetInfeasibleConstraints(...). In the deprecating message, we mention by certain date (say June 1), we will remove this deprecated function, and the user should all switch to GetInfeasibleConstraintNames.
2. After the deprecation date, we add the function std::vector<Binding<Constraint>> GetInfeasibleConstraint. Hopefully by the time we add this new function, our users have all switched to GetInfeasibleConstraintNames.

The alternative is to support two functions

std::vector<std::string> GetInfeasibleConstraintNames(...);
std::vector<Binding<Constraint>> GetInfeasibleBindings(...);

and deprecate std::vector<std::string> GetInfeasibleConstraint.

[jwnimmer-tri]

Somehow I lost the thread from adding debugging tips to needing a new API.

It seems to me like the only information in Binding that isn't in the error message already is the decision variable names?

As a starting point, would it be easier to add the name of the decision variable to the currently-returned string with no API changes? I think right now its returning "xdot_wall[0]: 0.0 <= 1.5 <= 1.0", for example. Is there a good place in there to add in the variable name(s) as well?

If we're going to change the API, I think we need a more substantial rethink that the proposal above. This functionality should not be in solve.h, and should probably involve new supporting classes instead of more free functions piling up. For example an InfeasibilityReport class that can provide the list of bindings, or string, or constraint+vars+failed_indices+failed_values tuples, or etc.

[hongkai-dai]

@jwnimmer-tri sorry I have been silent for two weeks.

I think the information needed for the infeasibility report include

1. A vector of Binding<Constraint> containing the unsatisfied constraint.
2. The string description of these constraints.
3. The value of the unsatisfied constraints evaluated at the solution.

We can obtain (2) (the string description) from (1) (the Binding). So this InfeasibilityReport class could be a container of std::vector<std::pair<Binding<Constraint>, Eigen::VectorXd>>. Each pair contains the unsatisfied constraint, together with the value of the constraint evaluated at the solution.

[jwnimmer-tri]

I was trying to nudge you to a simpler solution. The complaint was that GetInfeasibleConstraints returns some strings, and the strings don't have enough information. They lack the binding-specific details (which variables they applied to).

Isn't the simplest solution to just enhance GetInfeasibleConstraints to also incorporate the variables into the strings its already returning?

InfeasibilityReport is a lot of work, if the strings that users are already printing can just be amended to be sufficient on their own.

[hongkai-dai]

@jwnimmer-tri Sorry I didn't make it clear. From my experience on debugging nonlinear optimization, I think I sometimes need more than returning the string. I want the Binding<Constraint> object that are not satisfied at the infeasible solution. I sometimes evaluate this binding to check if it has the correct gradient (some constraint don't compute gradient through AutoDiff, such as the kinematic constraint in multibody/inverse_kinematics folder).

[jwnimmer-tri]

If you need the full report class, so be it. But if there is a conventional trick to debugging, maybe the current helper (strings) could incorporate that information as well. Report is a lot of work. If there is one best answer (the variables, their gradients, etc.) then maybe just print that into the strings instead? We don't want to have users have to copy-pasta a whole debugging interactions with the Report class, if the suggestions can be computed automatically and by default.

[hongkai-dai]

Thanks Jeremy, I think I want to propose this API

std::vector<std::pair<Binding<Constraint>, Eigen::VectorXd>> MathematicalProgramResult::GetInfeasibleConstraints(std::optional<double> tolerance) const;

So I return a pair of Bindings, together with the values of the Binding at the solution. We could extract the string description from each Binding object.

[jwnimmer-tri]

To avoid API explosion, would it be simpler to only return the infeasible bindings? The user can still evaluate them if they want details on the values. (If the x_val is not set, I don't think we should be returning evaluated bindings anyway.)

[hongkai-dai]

To avoid API explosion, would it be simpler to only return the infeasible bindings? The user can still evaluate them if they want details on the values. (If the x_val is not set, I don't think we should be returning evaluated bindings anyway.)

We get the infeasible bindings by evaluating each binding with x_val, and then check if the value of the binding is within the bounds. So if x_val is not set, we should throw an error when calling GetInfeasibleConstraints.

The reason why I propose to return the value of the Binding is that computing some of the binding can be time consuming (if the evaluator uses some very complicated function), and we already have the Binding value to obtain the infeasible constraint.

[jwnimmer-tri]

Why does the computational expense matter? We are producing output for user eyeballs. The constraints were already evaluated hundreds to thousands of times. One more can't possibly matter?

[hongkai-dai]

Makes sense. I can then call result.EvalBinding(binding) to return the value of the binding very easily.

I realize this GetInfeasibleConstraint function cannot belong to MathematicalProgramResult. This function will need to access all constraints inside MathematicalProgram, but MPResult doesn't link to MP (otherwise it would cause circular dependence). So I still need to put this function in solve.h (which includes both mathematical_program.h and mathematical_program_result.h.

In solve.h we already have a function GetInfeasibleConstraint that returns std::vector<std::string>. I cannot add another GetInfeasibleConstraint function but just change its return types. What would be a good process to deprecate the old GetInfeasibleConstraint function (returning string) and add the new GetInfeasibleConstraint function?

[jwnimmer-tri]

The only functions that should be in solve.h are those that require a dependency on all possible solvers, e.g. ChooseBestSolver, Solve, then maybe GetSolverIdFromName in the future, etc.

Characterizing a solution does not meet that bar. The current GetInfeasibleConstraints is already misplaced. So any change here should also fix the file.

The new function name could be GetInfeasibleConstraintBindings. Once the old one is gone, we could rename it back to be shorter.

[RussTedrake]

Just discussed with @hongkai-dai again about landing this tutorial (the original topic of the post). It would be very valuable.

Here is a quick list of debugging techniques that we came up with:

• name your costs/constraints/variables (via set_description)
• call GetInfeasibleConstraints
• use the visualization_callback
• set solver verbosity level
• use EvalBinding to test individual constraints

plus some general advice about commenting out some costs/constraints to understand the simplest version of the problem first.