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)$$?

• @copper.hat I've amended the question.
– Leo
May 14 '21 at 19:27

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 = A*x + [0;1.8*x2^2*x3];
g = x3*(1+x2^2) - 1;
V = x'*P*x;
dV = jacobian(V,x)*f;
[s,coeffs] = polynomial([x;x3],1);
Q = sdpvar(2);
solvesos([Q >= .1*eye(2), sos(-x'*Q*x + 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

• Can’t you show it analytically?
– Leo
May 15 '21 at 9:54
• A slightly easier approach (IMO) is to just show that $-\dot{V}(x)(x_2^2 + 1)$ is SOS. Then you don't need that artifical third state and you can even get a quite nice $Q$ matrix (with integer coefficients). That should also classify as shown analytically. May 15 '21 at 11:02
• Yes, much easier, so I added that. May 15 '21 at 13:51

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.