Drake: Mosek/MosekSolver numerical issues lead to false confidence in RegionOfAttraction estimates

[RussTedrake]

In adding the following example to underactuated:

To demonstrate that, let's make a system with a known, non-convex region of attraction. We'll do this by taking some interesting potential function 𝑈(𝑥)∈𝑆𝑂𝑆 and setting the dynamics to be 𝑥˙=(𝑈(𝑥)−1)∂𝑈∂𝑥𝑇 , which has 𝑈(𝑥)<=1 as the region of attraction.

import numpy as np
import matplotlib.pyplot as plt
from pydrake.all import (Variable, Jacobian, SymbolicVectorSystem, 
                         RegionOfAttraction, RegionOfAttractionOptions, 
                         plot_sublevelset_expression)

# Construct a non-convex 2D level set.
x = np.array([Variable("x"), Variable("y")]).reshape((2,))
A1 = np.array([[1, 2], [3, 4]])
A2 = A1 @ np.array([[-1, 0], [0, 1]])  # mirror about y-axis
U = (x.T @ A1.T @ A1 @ x) * (x.T @ A2.T @ A2 @ x)

fig, ax = plt.subplots()

dUdx = U.Jacobian(x)
sys = SymbolicVectorSystem(state=x, dynamics=(U-1)* dUdx.T)
context = sys.CreateDefaultContext()

options = RegionOfAttractionOptions()
options.lyapunov_candidate = x.dot(x)
options.state_variables = x
V = RegionOfAttraction(sys, context, options)
plot_sublevelset_expression(ax, V)
plot_sublevelset_expression(ax, U, 101, linewidth=3, fill=False);

# Note: fine to remove this part if you want to reproduce w/o underactuated as a dependency
from underactuated import plot_2d_phase_portrait
plot_2d_phase_portrait(sys, (-.8, .8), (-.6, .6))

I was disappointed to see the solver happily report that it has discovered the RoA, but that is signficantly outside the known RoA. (Black line is known, blue ellipse is estimated).

I've also constructed the C++ unit test that captures the issue as a reproduction. It currently fails on master if we remove the guard disabling it for Mosek.

@shensquared took a quick look and agrees with me that the formulation is correct, it is only numerical issues getting involved. She observed that setting d=1 in my RoA code resolves this issue. Taking that as an option could make the test pass, but it doesn't address my fundamental concern (that I've been lied to by the solver!).

Is there a general way to know, even after the solve has finished, that my certificates were invalid? Is there a solver parameter that could have improved the results. My quick playing with mosek parameters was not encouraging.

@frankpermenter -- i wonder if you have any thoughts?
cc @hongkai-dai , @TobiaMarcucci

[hongkai-dai]

@RussTedrake do you still have the matlab code with this example? I am wondering if using spotless would give us the right result. I suspect it could be an issue when I convert MathematicalProgram to mosek format in mosek_solver.cc.

[RussTedrake]

No, I don't have a matlab version... though I suppose it wouldn't be too hard to make. I made this one up fresh. The recipe is the same one that Mark T used for his example that I wanted to recreate, but when I checked a few years back, Mark said that he couldn't find the code. But the particular form -- and therefore the numerics -- are new.

[hongkai-dai]

I solved the problem with CSDP instead of Mosek, and here is the plot

So the problem is in either our mosek_solver.cc or in the mosek itself. I will try to create a matlab version with spotless, when I figure out how to run mosek solver on my VPN (it complains "Could not acquire a mosek license when I use TRI VPN connection).

[shensquared]

@hongkai-dai I have spotless and Mosek all setup on my local machine; I'll take care of checking matlab/spotless result.

[shensquared]

For this example, on matlab, when all parameters are set the same way as the C++ version, including all degrees (pre-multiplier d degree is 7, L degree is 8), and how the multiplier is setup (declared Free and then add SOS constraint, instead of NewSos directly), Mosek verbose prints Unknown error.

When keeping the L degree fixed at 8, and reduce degree d=6, solver reports primal and dual feasible, but still gives wrong result, plot below.

When stepping d further down from 5 to 1, Mosek gives unknown error except for when d = 3 and d =2. In these two cases, Mosek reports primal and dual feasible, yet the verified ROA still goes beyond the ground truth (slightly less wrong than the above plot where d =6).

Note that, the default quadratic balancing has to be disabled. This is because the dynamic linearization A at the origin is the all-zero matrix. This causes the 2nd argument in https://github.com/RobotLocomotion/drake/blob/d7f3c011d37d471d7b9293ecf2066c98d88b2a05/drake/matlab/systems/%40PolynomialSystem/regionOfAttraction.m#L159 to be exactly zero matrix too, which then causes https://github.com/RobotLocomotion/drake/blob/d7f3c011d37d471d7b9293ecf2066c98d88b2a05/drake/matlab/util/balanceQuadForm.m#L16 to fail.

[RussTedrake]

Thanks Shen! Ironically, I just ported the balancing code this weekend (see #12875), and realized the oversight in my old code when I re-derived the equations. My PR is careful to NOT do balancing in the indefinite case.

[shensquared]

Just tried Russ’s script, using python binding with the new C++ ROA code. The MOSEK problem size looks off, (too small, and no sdp variable?), see below.

Moreover, if put in van-der-pol dynamics with linearized/quadratic V instead, MOSEK reports the same problem size, but it shouldn’t - the two examples have very different degf, degVdot, degL. Strangely, VDP solution looks about right despite the off-size.

---

If directly writing the ROA formulation in python, the problem sizes, below, are much more sensible.

These sizes are consistent with spotless outputs; the solution are ‘better’ than spotless. Here, Mosek reports numerical unknown error for all d except for d =1 and 2. In these two `solved’ cases, Mosek returns numerically-zero rho, so at least it's correct.

[frankpermenter]

@RussTedrake A standard error measure for SDPs are the DIMACS errors:
http://plato.asu.edu/dimacs/node3.html

They are actually computed in spotless (https://github.com/spot-toolbox/spotless/blob/master/spotopt/util/spot_sdp_dimacs.m).