Region of Attraction Estimation for Polynomial Systems

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:

\[\begin{aligned} -\dot V(\mathbf{x}) + \lambda(\mathbf{x})\bigl(V(\mathbf{x}) - \rho\bigr) &\text{ is SOS},\\ V(\mathbf{x}) - \epsilon\,\mathbf{x}^\top\mathbf{x} &\text{ is SOS},\\ \lambda(\mathbf{x}) &\text{ is SOS}, \epsilon>0, \end{aligned}\]

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

Bi-level Optimization

Primal problem:

\[\begin{aligned} \max_{\rho,\,V,\,\lambda}\quad & \rho \\[3pt] \text{s.t.}\quad & -\dot V(\mathbf{x}) + \lambda(\mathbf{x})\bigl(V(\mathbf{x})-\rho\bigr)\ \text{is SOS},\\ & \lambda(\mathbf{x})\ \text{is SOS},\\ & V(\mathbf{x})-\epsilon\,\mathbf{x}^{\!\top}\mathbf{x}\ \text{is SOS},\quad \epsilon>0 . \end{aligned}\]

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\))

\[\begin{aligned} \max_{\lambda(\mathbf{x}),\,s}\quad & s \\[4pt] \text{s.t.}\quad & s>0,\\ & -\dot V(\mathbf{x}) + \lambda(\mathbf{x})\bigl(V(\mathbf{x})-\rho\bigr) - s \ \text{is SOS},\\ & \lambda(\mathbf{x})\ \text{is SOS}. \end{aligned}\]

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

\[\begin{aligned} \max_{V(\mathbf{x}),\,\rho}\quad & \rho \\[4pt] \text{s.t.}\quad & \rho>0,\\ & -\dot V(\mathbf{x}) + \lambda(\mathbf{x})\bigl(V(\mathbf{x})-\rho\bigr) \ \text{is SOS},\\ & V(\mathbf{x})-\epsilon\,\mathbf{x}^{\!\top}\mathbf{x} \ \text{is SOS}. \end{aligned}\]

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.

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. We implement the SOS programs in MATLAB using SPOT (Systems Polynomial Optimization Tools) and solve the resulting SDPs with SeDuMi.

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