Finding a region of attraction using sum of squares optimization

Specific Domains Optimization (Mathematical)
,
1

Out of interest I decided to explore sum of squares optimization in the context of Lyapunov functions. Right now I’d like to use it to find the region of attraction for the value function (viewed as a Lyapunov function) of the LQR. As an example I have a single pendulum which I am keeping in its unstable equilibrium using the regulator. As for the approach I am trying to use the “The equality-constrained formulation” in section 9.2.3 of the Underactuated MIT course.

I’ve written the problem using SumOfSquares.jl but even though the simulation shows that the controlled system is stable the region found is practically zero. The script is not too long (83 lines) and hopefully commented thoroughly enough. Any help would be appreciated.

It is quite possible that there is some theoretical mistake based on the fact that I am searching over a region that is a cylinder but to my understanding that shouldn’t be the case.

2

This is a question for @blegat

3

I suspect (but I am not sure) that the culprit might lie hidden in the way the lifted state variables (c_\theta, s_\theta, \omega) are transformed into the physical state variables (\theta, \omega) around the upper (unstable) equilibrium. Try to run your simulation a bit longer, say, set tspan = (0.0, 10) on line 47 of your code. You will observe that once the pendulum is around the upper position, even tiny oscillations of c_\theta and s_\theta around their equilibrium values lead to a jump in the angle \theta.

By the way, you may have wanted to put -1 instead of 1 on line 54: plt = plot(ts, mapreduce(t -> state_difference(sol(t), [0, 1, 0])', vcat, ts)).

4

I think that should only effect the plotting, where I have theta=0 at the stable equilibrium. As I linearize the system around the upright equilibrium the controller does behave as it should. I could Indeed somewhat fix it by selecting theta=0 to be at the unstable equilibrium for plotting.

I stripped down the script down a bit. The lyapunov function is now -potential_energy + kinetic_energy. For it, I created the controller ad-hoc. Its first and third term would stabilize the system from any state beside theta=0 so I added the second term to limit the maximum potential ROA. I still run into the same problem.

5

Mosek terminates early with the status SLOW_PROGRESS for some reason which is probably why the region it too large. I wonder what makes it get stuck and if it could be helped in any way.

6

From my experience, nonlinear conic optimization has its intrinsic severer numerical issues than linear optimization. Therefore it is not unusual to see Mosek slowing its progress.

7

By the way, didn’t you mean
@constraint(model, pos * (V - ρ) + μ * (s^2 + c^2 - 1) + λ * dVdt >= 0)
instead of
@constraint(model, pos * (V - ρ) + μ * (s^2 + c^2) + λ * dVdt >= 0) ?

True, this does not solve the original trouble. But anyway.

8

Yeah, I also noticed it and fixed it.

9

Well, I admit did not have prior experience with that equality variant of the ROA optimization problem described in the Tedrake’s notes. But have you tried the more “traditional” inequality one (also described in the notes) for which the constraint would be -\dot V(x) + \lambda(x)(V(x)-\rho) \in \text{SOS}, and \lambda(x)\in\text{SOS}? Of course, you will have to proceed by line-searching for \rho (fixing \rho to some positive value), but that is not a big deal. I am just curious if at least the other approach worked for this system.

10

Have you tried another solver? I tried CSDP, this is what I got

CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 8.32e-01 Pobj: -7.7148538e+01 Ad: 7.66e-01 Dobj:  2.5109508e+01 
Iter:  2 Ap: 8.34e-01 Pobj: -1.3583956e+02 Ad: 8.37e-01 Dobj:  1.4435212e+01 
Iter:  3 Ap: 9.15e-01 Pobj: -8.4140204e+01 Ad: 8.92e-01 Dobj:  2.2222230e+01 
:
:
:
Iter: 40 Ap: 4.99e-01 Pobj:  4.9164478e-03 Ad: 8.28e-01 Dobj:  2.4595803e-03 
Iter: 41 Ap: 6.91e-01 Pobj:  4.9436046e-03 Ad: 6.61e-01 Dobj:  2.4595351e-03 
Iter: 42 Ap: 8.80e-01 Pobj:  4.9130371e-03 Ad: 9.80e-01 Dobj:  2.4595125e-03 
Iter: 43 Ap: 8.09e-02 Pobj:  4.9143295e-03 Ad: 8.77e-01 Dobj:  2.4594533e-03 
Stuck at edge of primal feasibility, giving up. 
Partial Success: SDP solved with reduced accuracy
Primal objective value: 4.9157212e-03 
Dual objective value: 2.4594527e-03 
Relative primal infeasibility: 4.85e-08 
Relative dual infeasibility: 2.66e-10 
Real Relative Gap: -2.44e-03 
XZ Relative Gap: 9.62e-09 
DIMACS error measures: 7.93e-08 0.00e+00 4.13e-10 0.00e+00 -2.44e-03 9.62e-09

ALMOST_OPTIMAL::TerminationStatusCode = 7

0.004915721213810453
11

Yes, I tried alternating between searching for \lambda(x) and \rho but it again failed with SLOW PROGRESS (not only Mosek but also CSDP), searching for \lambda(x), in the first iteration (I’ve added the script to the gist).

I also tried Hypatia.jl which gave a similar result to Mosek but claimed to be successful. Regarding your result it’s way smaller than I would expect. I also tried searching for a region of attraction when I have the globally stabilizing controller with CSDP and the result was similarly small.

Now that I’ve been thinking about this topic for a bit, I’m not sure I see utility in quantifying the region of attraction using a single value, especially if my Lyapunov function is fixed. It feels like it is restricted by the how the coefficients scale different polynomials. Maybe this also causes some of the numerical issues.

The reason I bring this up is that for the same system I was able to do some cool stuff regarding optimizing a controller for the same Lyapunov function.