Region of Attraction Estimation for Polynomial Systems

← Back

TL;DR

For a nonlinear polynomial dynamical system \(\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}),\quad \mathbf{f}:\mathbb{R}^n\!\to\!\mathbb{R}^n,\) this project seeks a local Lyapunov function \(V(\mathbf{x})\succ 0\) and a scalar \(\rho>0\) whose sublevel set \(B(\rho)=\{\mathbf{x}\mid V(\mathbf{x})\le\rho\}\) is an inner estimate of the Region of Attraction (ROA) to the origin. The goal is to maximize \(\rho\) while guaranteeing the implication \(V(\mathbf{x})\le\rho\;\Longrightarrow\;\dot{V}(\mathbf{x})<0.\)

S-Procedure and Sum-of-Squares Relaxation

Define the candidate Lyapunov sublevel set where maxmizing largest feasible \(\rho>0\): \(B := \Bigl\{\mathbf{x}\,\Bigl|\,V(\mathbf{x}) \le \rho\Bigr\}, \dot V(\mathbf{x}) = \frac{\partial V}{\partial \mathbf{x}}\,\mathbf{f}(\mathbf{x}).\) Our objective is to certify the implication: whenever \(V(\mathbf{x}) \le \rho\), we have \(\dot V(\mathbf{x}) \le 0\), \(V(\mathbf{x}) \le \rho \;\Longrightarrow\; \dot V(\mathbf{x}) \le 0.\) S-procedure: Introduce a polynomial multiplier \(\lambda(\mathbf{x}) \ge 0\). If \(\dot V(\mathbf{x}) \;\le\; \lambda(\mathbf{x})\bigl(V(\mathbf{x})-\rho\bigr),\) then the desired implication holds for all \(\mathbf{x}\). SOS relaxation: Cast the conditions as sum-of-squares (SOS) constraints:

Each SOS constraint ensures global non-negativity; the last enforces \(V(\mathbf{x}) \succ 0\).

Bi-level Optimization

Primal problem:

Primal problem is non-convex; the following two sub-problems decouple the bilinear term and are solved in alternation.

Step A — Multiplier update (fixed \(V(x),\rho\))

Step B — Lyapunov update (fixed \(s\))

Initialization.
Choose a feasible quadratic seed \(V_0(\mathbf{x})=\mathbf{x}^{\!\top}P\mathbf{x}\) from a quadratic Lyapunov equation or LQR, then alternate Step A/B until \(\rho\) converges.

We already saw the standard regional SOS certificate $-\dot V+\lambda(V-\rho)\in\text{SOS},\ \lambda\in\text{SOS}$, whose nuisance is the bilinear $\lambda\rho$ (and $\lambda V$ if $V$ is variable). In [1] Section 9.2.3 and UR Deepnote notebook is to move from a semi-algebraic set $V\le\rho$ to the algebraic variety ${\dot V=0}$, and encode the “except $x=0$” clause via a vanishing factor.

For equality constraints, $\exists\,\lambda(x)$ with $p(x)+\lambda(x)g(x)\in\text{SOS}\Rightarrow p(x)\ge0$ on ${g=0}$; $\lambda$ need not be SOS.

Naively taking $g=\dot V,\ p=V-\rho$ fails because $0\in{\dot V=0}$ and $V(0)=0\Rightarrow \rho\le0$. Fix by multiplying $p$ with $(x^\top x)^d$ ($d\in\mathbb N_+$):

On $\dot V=0$: $(x^\top x)^d(V-\rho)\ge0$, i.e. $V\ge\rho$ for all $x\neq0$ with $\dot V=0$; the origin is automatically exempt since $(x^\top x)^d=0$ at $x=0$. With the usual local condition ($\nabla^2\dot V(0)\prec0$), this excludes any nontrivial $\dot V=0$ crossing inside $V<\rho$, hence yields $\dot V<0$ on $V\le\rho\setminus{0}$.

Crucially, in (★) we have $(x^\top x)^d(V-\rho)=(x^\top x)^dV-\rho(x^\top x)^d$: $\rho$ enters linearly (no $\lambda\rho$).

With $V$ fixed, solve one SOS program:

No line search; $\lambda$ unconstrained-sign; $d$ is chosen for degree parity/matching.

Searching $V$: bilinear alternation

If $V$ is also variable, $\lambda\,\dot V$ is bilinear in $(\lambda,V)$. Use block-coordinate convexification:

1. Multiplier/level update (fix $V$)

   \[\max_{\rho,\lambda}\ \rho\ \ \text{s.t.}\ \ (x^\top x)^d(V-\rho)+\lambda\,\dot V\in\text{SOS}.\]

   Rescale $V\leftarrow V/\rho$ so the certified set is $V\le1$.

2. Lyapunov update (fix $\lambda$)
   For quadratic $V=x^\top P x$:

   \[\min_{P\succeq0}\ \operatorname{tr}(\hat P^{-1}P)\quad \text{s.t.}\quad (x^\top x)^d(V-1)+\lambda\,\dot V\in\text{SOS},\]

   where $\operatorname{tr}(\hat P^{-1}P)$ is the first-order surrogate for $\log\det P$ (ellipsoid volume). Iterate 1 $\leftrightarrow$ 2 to convergence.

Implementation on Van der Pol

The nonlinear system: \(\dot x_1 = -x_2,\\ \dot x_2 = x_1-(1-x_1^{2})x_2\) has a limit cycle, making inner ROA enlargement non-trivial. The SOS programs is implemented in MATLAB using SPOT and solve the resulting SDPs with SeDuMi. UR Deepnote notebook use Mosek (in pydrake).

The certified region stabilizes after \(\sim 23\) alternations and remains strictly inside the true ROA.

References

[1]. Russ Tedrake. Underactuated Robotics Ch.9 Lyapunov Analysis Course Notes for MIT 6.832, 2023. https://underactuated.csail.mit.edu/lyapunov.html#roa_cubic_system