Quadratic Lyapunov function

Question

Consider the nonlinear model of a chemical process:

$$ \begin{align} \dot{x}_1 &= x_2\\ \dot{x}_2 &=-5x_1 -3x_2\left(1 - \frac{0.6x_2}{1 + x_2^2}\right). \end{align} $$

Writing this system in terms of

$$\dot{x}_1 = Ax +g(x) $$

We have:

\begin{align} \dot{x}= \begin{pmatrix} 0 & 1 \\ -5 & -3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} + \begin{pmatrix} 0 \\ \frac{1.8\,x_2^2}{1 + x_2^2} \end{pmatrix} \end{align}

Now solving the Lyapunov equation

$$ A^\top P + P\,A =-I $$

where $I$ is the identity matrix we have $P$ as:

\begin{pmatrix} 13/10 & 1/10 \\ 1/10 & 1/5 \end{pmatrix}

Since we have $P$, how can we use the quadratic Lyapunov function

$$ V(x) = x^\top P\,x $$

to show asymptotic global stability using the quadratic Lyapunov function $V(x)$?

Answer 1

One popular computational approach is to use a sum-of-squares (SOS) strategy.

SOS typically works with polynomial systems, but you can easily outer approximate your dynamics by introducing a new fictious state $x_3$ and study stability of

\begin{align} \dot{x}= \begin{pmatrix} 0&1\\ -5 &-3\end{pmatrix} \begin{pmatrix} x_1\\x_2\end{pmatrix}+\begin{pmatrix} 0\\1.8x_3\end{pmatrix} \end{align}

on the set $(1+x_2^2)x_3=1$ which we denote $g(x)=0$

Hence, we want to prove that $\dot{V}(x_1,x_2) \leq -q(x_1,x_2)$ on $g(x)$. A sufficient condition for this is the existence of a polynomial $s(x)$ such that $\dot{V}(x) \leq -q(x) + s(x)g(x)$. In other words, if we can prove $-q(x) + s(x)g(x)-\dot{V}(x)$ is positive semidefinite we are done. Replace non-negativity with being sum-of-squares and you are done.

All you need is a method to compute a SOS-decomposition, and there are many alternatives for that. Here is a complete implementation in the MATLAB toolbox YALMIP. I search for a linear multiplier $s(x)$ and use a quadratic positive definite bound $q(x_1,x_2)$ on the derivative. The problem is feasible thus proving feasibility.

P = [13/10 1/10;1/10 1/5];
A = [0 1;-5 -3];
sdpvar x1 x2 x3
x = [x1;x2];
f = Ax + [0;1.8x2^2x3];
g = x3(1+x2^2) - 1;
V = x'Px;
dV = jacobian(V,x)f;
[s,coeffs] = polynomial([x;x3],1);
Q = sdpvar(2);
solvesos([Q >= .1eye(2), sos(-x'Qx + g*s - dV)],[],[],coeffs) 

https://yalmip.github.io/tutorial/sumofsquaresprogramming

Then you can be clever like SampleTime below (nice numbers so this can likely be done by hand)

[~,v,Q] = solvesos(sos(-x'*P*(A*x*(1+x2^2) + [0;1.8*x2^2])));
sdisplay(v{1}'*Q{1}*v{1})
0.5*x2^2+0.5*x1^2-0.36*x2^3-0.18*x1*x2^2+0.5*x2^4+0.5*x1^2*x2^2
min(eig(Q{1}))
ans =
0.3094

Answer 2

In this one can use a conservative way of showing stability, using the circle criterion. A more detailed description of this criterion can also be found in Nonlinear Systems by Hassan K. Khalil. Namely, your system can be written in the form

\begin{align} \dot{x} &= A\,x + B\,u, \tag{1a} \\ y &= C\,x + D\,u, \tag{1b} \\ u &= -\psi(t,y), \tag{1c} \end{align}

where $\psi(t,y)$ is sector bounded by $[K_1, K_2]$, such that $K = K_2 - K_1 = K^\top \succ 0$ and

$$ \left(\psi(t,y) + K_1\,y\right)^\top \left(\psi(t,y) + K_2\,y\right) \leq 0. \tag{2} $$

This system is stable for any $\psi(t,y)$ that satisfies $(2)$ if the following transfer function is strictly positive real

$$ Z(s) = \left(I + K_2\,G(s)\right) \left(I + K_1\,G(s)\right)^{-1}. \tag{3} $$

Where $G(s) = C\,(I\,s - A)^{-1}\,B + D$ is the transfer function associated with the linear time invariant state space model from $(1)$. It can be noted that strictly positive real means that $\text{Re}(Z(j\,\omega)) > 0$ for all $-\infty < \omega < \infty$.

If this condition, of strictly positive realness, is satisfied guarantees that there is a solution to the Kalman–Yakubovich–Popov equations

\begin{align} P\,A + A^\top P &= -L^\top L - \epsilon\,P, \tag{4a} \\ P\,B &= C^\top - L^\top W, \tag{4b} \\ W^\top W &= D + D^\top, \tag{4c} \end{align}

with $P = P^\top \succ 0$, $\epsilon > 0$ and $(A,B,C,D)$ matrices corresponding to a state space model corresponding to the transfer function $Z(s)$ from $(3)$. A solution to $(4)$ ensures that the following quadratic Lyapunov equation shows exponential stability

$$ V(x) = \frac{1}{2} x^\top P\,x, \tag{5} $$

since it can be shown that

$$ \dot{V} \leq -\frac{1}{2} \epsilon\,x^\top P\,x. \tag{6} $$

---

Your system can be written in the form of $(1)$ using

$$ A = \begin{bmatrix} 0 & 1 \\ -5 & -3 \end{bmatrix}, \quad B = \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \quad C = \begin{bmatrix} 0 & 1 \end{bmatrix}, \quad D = 0, \tag{7} $$

for the linear state space model and for the nonlinearity

$$ \psi(t,y) = -\frac{1.8\,y^2}{1 + y^2}. \tag{8} $$

It can be shown that $(8)$ satisfies $(2)$ using $[K_1, K_2] = [-0.9, 0.9]$. Substituting this into $(3)$ yields

$$ Z(s) = \frac{s^2 + 3.9\,s + 5}{s^2 + 2.1\,s + 5}, $$

which can be shown to be positive real by looking at its Bode plot, since its phase remains between -90 and 90 degrees. Thus it should be possible to find a quadratic Lyapunov function which shows exponential stability.

In order to keep the same state $x$ from $(1)$ using $(7)$ as for the state space representation of $Z(s)$ one can use the following state space model

\begin{align} \dot{x} &= (A - B\,K_1\,C)\,x + B\,u, \\ y &= K\,C\,x + u. \end{align}

It can be noted that if in $(1)$ the matrix $D \neq 0$ this transformation would look a little more complicated.

Using these resulting matrices in $(4)$ allows one to obtain the following solution

\begin{align} P &= \begin{bmatrix} 4.867 & 0.8117 \\ 0.8117 & 0.8121 \end{bmatrix}, \\ L &= \begin{bmatrix} -0.5739 & 0.6986 \end{bmatrix}, \\ W &= 1.4142, \\ \epsilon &= 1.6. \end{align}

Thus using this $P$ in $(5)$ would give a quadratic Lyapunov function for the initial nonlinear system, for which $(6)$ is satisfied.

In order to solve $(4)$ I wrote some quick and ugly code, but I believe there are dedicated solvers made specifically for this. Also note that showing that $Z(s)$ is strictly positive real is already sufficient to show stability, but if you want an explicit description of a Lyapunov function one would still have to solve $(4)$. Furthermore, the circle criterion is a sufficient but not necessary condition. So if $Z(s)$ is not strictly positive real does not imply the system is unstable.