2
$\begingroup$

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

$\endgroup$
1
  • 1
    $\begingroup$ @copper.hat I've amended the question. $\endgroup$
    – Leo
    May 14 '21 at 19:27
2
$\begingroup$

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
$\endgroup$
3
  • $\begingroup$ Can’t you show it analytically? $\endgroup$
    – Leo
    May 15 '21 at 9:54
  • $\begingroup$ 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. $\endgroup$
    – SampleTime
    May 15 '21 at 11:02
  • $\begingroup$ Yes, much easier, so I added that. $\endgroup$ May 15 '21 at 13:51
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.