# Finding a region of attraction using sum of squares optimization

**URL:** <https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, controlsystems\
**Created:** [July 25, 2025, 9:05pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067 "2025-07-25T21:05:25Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![lieskjur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lieskjur/32/27611_2.png) [@lieskjur](https://discourse.julialang.org/u/lieskjur)\
**Post date:** [July 25, 2025, 9:05pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/1 "2025-07-25T21:05:25Z")

</div>

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](https://underactuated.csail.mit.edu/lyapunov.html#optimization).

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](https://gist.github.com/JurajLieskovsky/f67379121925634c38c2d3d3aaa260c6) 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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 27, 2025, 3:20am UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/2 "2025-07-27T03:20:08Z")

</div>

This is a question for @blegat

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 28, 2025, 8:01am UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/3 "2025-07-28T08:01:14Z")

</div>

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))`.

---

<div class="post-metadata">

**Author:** ![lieskjur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lieskjur/32/27611_2.png) [@lieskjur](https://discourse.julialang.org/u/lieskjur)\
**Post date:** [July 28, 2025, 1:26pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/4 "2025-07-28T13:26:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lieskjur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lieskjur/32/27611_2.png) [@lieskjur](https://discourse.julialang.org/u/lieskjur)\
**Post date:** [July 28, 2025, 1:34pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/5 "2025-07-28T13:34:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 28, 2025, 1:40pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/6 "2025-07-28T13:40:06Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 28, 2025, 2:13pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/7 "2025-07-28T14:13:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lieskjur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lieskjur/32/27611_2.png) [@lieskjur](https://discourse.julialang.org/u/lieskjur)\
**Post date:** [July 28, 2025, 2:15pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/8 "2025-07-28T14:15:30Z")

</div>

Yeah, I also noticed it and fixed it.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 28, 2025, 2:44pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/9 "2025-07-28T14:44:41Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [July 28, 2025, 3:03pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/10 "2025-07-28T15:03:27Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![lieskjur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lieskjur/32/27611_2.png) [@lieskjur](https://discourse.julialang.org/u/lieskjur)\
**Post date:** [July 28, 2025, 5:46pm UTC](https://discourse.julialang.org/t/finding-a-region-of-attraction-using-sum-of-squares-optimization/131067/11 "2025-07-28T17:46:18Z")

</div>

> [@zdenek\_hurak](#):
>
> But have you tried the more “traditional” _inequality_ one

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](https://gist.github.com/JurajLieskovsky/f67379121925634c38c2d3d3aaa260c6)).

> [@zdenek\_hurak](#):
>
> Have you tried another solver?

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.