Lyapunov Equation: How to Compute the Solution and Origins in Stability Analysis of Dynamical Systems – Fusion of Engineering, Control, Coding, Machine Learning, and Science

In this tutorial, you will learn the following topics:

Derivation of the Lyapunov equation. You will learn to derive the Lyapunov equation by analyzing the stability of dynamical systems.
Stability test of dynamical systems that is based on solving the Lyapunov equation.
Numerical solution of the Lyapunov equation. We will use a relatively simple method that is based on vectorization, as well as a more advanced method implemented in MATLAB’s lyap() function.

Origins of the Lyapunov Equation in Stability Analysis of Dynamical Systems

First, let us talk about the origins of the Lyapunov equation. The Lyapunov equation is a fundamental equation that is used for the stability and performance analysis of dynamical systems. Apart from the stability analysis, the solution of the Lyapunov equation is also used in other control algorithms and methods. For example, in order to compute the solution of an optimal control problem, it is often necessary to solve Lyapunov equations.

In this post, for brevity and simplicity, we explain the derivation of the Lyapunov equation that is based on the stability analysis of dynamical systems. Let us consider a dynamical system of the following form:

(1) $\begin{align*}\dot{\mathbf{x}}=A\mathbf{x}\end{align*}$

where $\mathbf{x}\in \mathbb{R}^{n}$ is the system’s state vector and $A\in \mathbb{R}^{n\times n}$ is the system state transition matrix. This linear system can also come from the linearization of the nonlinear system.

Let us assume that the system (1) has a unique equilibrium point. Then, the stability analysis is concerned with the following question:

What is the condition under which the system’s equilibrium point is asymptotically stable?

To answer this question, we will briefly recall the Lyapunov stability theorem.

Lyapunov stability theorem: Let $V(\mathbf{x})$ be a function that maps $\mathbf{x}\in \mathbb{R}^{n}$ into a real variable. Then, let $\dot{V}(\mathbf{x})$ be the first derivative of this function along the state trajectories of the system (1). If there exists some subset of $\mathbb{R}^{n}$ such that:

The function $V(\mathbf{x})$ is positive definite on this subset. This means that 1) $V(\mathbf{x})=0$ if and only if $\mathbf{x}=0$ , and 2) $V(\mathbf{x})>0$ if and only if $\mathbf{x}\ne 0$ .
The function $\dot{V}(\mathbf{x})$ is negative (negative definite) on this subset.

Then, the equilibrium point of the system (1) is asymptotically stable.

The function $V(\mathbf{x})$ satisfying these two conditions is called the Lyapunov function. Here the Lyapunov function should not be confused with the Lyapunov equation that is introduced in the sequel.

Quadratic Forms, Positive Definite, Negative Definite, and Semi-Definite Matrices

The Lyapunov equation for the linear system (1), will be derived by assuming a Lyapunov function with a quadratic form. Consequently, it is important to explain quadratic forms and relate these forms with positive definite matrices. Let us consider the following function

(2) $\begin{align*}V(\mathbf{x})=\mathbf{x}^{T}P\mathbf{x}\end{align*}$

where $P$ is a symmetric matrix, that is, $P=P^{T}$ . This function represents a quadratic form written in a compact matrix form. Function $V(\mathbf{x})$ , where $\mathbf{x}\in \mathbb{R}^{2}$ for several different forms of the matrix $P$ (that are given in Eq. (8)), is shown in Figs. 1-5 below. For example, consider this function:

(3) $\begin{align*}V(x_{1},x_{2})=2x^{2}_{1}+3x^{2}_{2}\end{align*}$

where $x_{1}$ and $x_{2}$ are real variables. This function is actually a quadratic form since

(4) $\begin{align*}V(\mathbf{x})=2x^{2}_{1}+3x^{2}_{2}=\begin{bmatrix} x_{1} & x_{2} \end{bmatrix}\underbrace{\begin{bmatrix}2 & 0 \\ 0 & 3 \end{bmatrix}}_{P}\begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix} = \mathbf{x}^{T}\underbrace{\begin{bmatrix}2 & 0 \\ 0 & 3 \end{bmatrix}}_{P} \mathbf{x}\end{align*}$

where

(5) $\begin{align*}\mathbf{x}=\begin{bmatrix}x_{1}\\x_{2} \end{bmatrix}\end{align*}$

In order to better illustrate the structure of quadratic forms, let us consider a 2D case below

(6) $\begin{align*}P=\begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix}\end{align*}$

and consequently, our original function takes the following form

(7) $\begin{align*} & V(\mathbf{x}) =\begin{bmatrix}x_{1} & x_{2} \end{bmatrix}\begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix}x_{1} \\ x_{2} \end{bmatrix} \\& = \begin{bmatrix}x_{1} & x_{2} \end{bmatrix} \begin{bmatrix}p_{11}x_{1} +p_{12}x_{2} \\ p_{12}x_{1} +p_{22}x_{2} \end{bmatrix} \\& = p_{11}x_{1}^2+p_{22}x_{2}^{2}+2p_{12}x_{1}x_{2} \end{align*}$

Let us visualize these functions, for several values of the matrix $P$ . Consider the following five $P$ matrices

(8) $\begin{align*}& P_{1}=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \\& P_{2}=\begin{bmatrix} 1 & -2 \\ 1 & 5 \end{bmatrix} \\& P_{3}=\begin{bmatrix} -1 & -2 \\ 1 & -5 \end{bmatrix} \\& P_{4}=\begin{bmatrix} -1 & -2 \\ 1 & 2 \end{bmatrix} \\& P_{5}=\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} \\\end{align*}$

The function $V(\mathbf{x})$ for these five cases is shown in Figs 1-5 below.

Figure 1: $V(\mathbf{x})$ for the $P_{1}$ matrix in (8).

Figure 2: $V(\mathbf{x})$ for the $P_{2}$ matrix in (8).

Figure 3: $V(\mathbf{x})$ for the $P_{3}$ matrix in (8).

Figure 4: $V(\mathbf{x})$ for the $P_{4}$ matrix in (8).

Figure 5: $V(\mathbf{x})$ for the $P_{5}$ matrix in (8).

The MATLAB code used to generate these plots is given below.

clear,pack,clc
x1=-50:2:50;
x2=-50:2:50;

[X1,X2]=meshgrid(x1,x2);

syms x1s
syms x2s

% P1 - positive definite completely symmetric with respect to the zero axis
%P=[1 0; 0  1];
% P2 - positive definite not completely symmetric
% P=[1 -2;1  5];
% P3 - negative definite 
%P=[-1 -2;1  -5];
%P4 - indefinite 
%P=[-1 -2; 1  2];
%P5 - positive semi-definite
P=[1 1; 1 1];

expression= [x1s x2s]*P*[x1s; x2s];

% print this expression in the command window to see the expanded form of
% the Z function

% use the expanded expression to define a new function
expand(expression)

V= X1.^2+X2.^2+(X1.*X2).*2

figure(1)
surf(V)
xlabel('x1')
ylabel('x2')
zlabel('V(x1,x2)')

Here, we need to mention one important fact about quadratic forms. Notice that in the definition of the quadratic form, given by Eq. (2), it is assumed that the matrix $P$ is symmetric. However, many matrices are not symmetric. Is it possible to define quadratic forms for such matrices? The answer is yes, since

(9) $\begin{align*}\mathbf{x}^{T}A\mathbf{x}=\mathbf{x}^{T}\underbrace{\Big(\frac{1}{2} \big( A+A^{T} \big)\Big)}_{\text{symmetric matrix}} \mathbf{x}=\mathbf{x}^{T}A_{1}\mathbf{x} \\A_{1}=\frac{1}{2} \big( A+A^{T} \big), \;\; A_{1} \;\; \text{is symmetric} \end{align*}$

Consequently, any non-symmetric matrix defines a quadratic form with a symmetric matrix. That is, for any non-symmetric matrix, we can find a quadratic form with a symmetric matrix.

From (7), we can see that the function $V(\mathbf{x})$ is quadratic. For this function to be positive definite, the matrix $P$ has to be a positive definite matrix. The symmetric real matrix $P$ is said to be a positive definite matrix if and only if, the following condition is satisfied

(10) $\begin{align*} \mathbf{x}^{T}P\mathbf{x}>0, \;\; \forall \mathbf{x}\in \mathbb{R}^{n} \setminus \{0\}\end{align*}$

Positive definite quadratic forms corresponding to the positive definite matrices $P_{1}$ and $P_{2}$ are illustrated in Figs. 1 and 2. Let us now discuss a test for testing positive definiteness of a matrix. Let us consider the 2D case given by Eq.(7). From this equation, we have:

(11) $\begin{align*}& V(\mathbf{x})=p_{11}x_{1}^2+p_{22}x_{2}^{2}+2p_{12}x_{1}x_{2} \\& =p_{11}\Big(x_{1}^2+\frac{p_{22}}{p_{11}}x_{2}^{2}+2\frac{p_{12}}{p_{11}}x_{1}x_{2} \Big) \\& = p_{11} \Big(x_{1}^2+2x_{1}\frac{p_{12}}{p_{11}}x_{2}+\frac{p_{12}^{2}}{p_{11}^{2}}x_{2}^{2}-\frac{p_{12}^{2}}{p_{11}^{2}}x_{2}^{2}+\frac{p_{22}}{p_{11}}x_{2}^{2} \Big) \\& = p_{11} \Big( \big(x_{1}+\frac{p_{12}}{p_{11}}x_{2} \big)^{2}+x_{2}^{2}\big(\frac{p_{22}}{p_{11}}-\frac{p_{12}^{2}}{p_{11}^{2}} \big) \Big) \end{align*}$

Let us analyze the last expression. Let us assume that $p_{11}$ is positive. Then, it is easy to see that the first term is positive for any $x_{1}$ and $x_{2}$ different from zero, that is

(12) $\begin{align*} p_{11}\big(x_{1}+\frac{p_{12}}{p_{11}}x_{2} \big)^{2}> 0,\;\; \text{for all} \;\; x_{1}\ne 0,\;\; x_{2}\ne 0\end{align*}$

Let us consider the second term

(13) $\begin{align*} p_{11}x_{2}^{2}\big(\frac{p_{22}}{p_{11}}-\frac{p_{12}^{2}}{p_{11}^{2}} \big) \end{align*}$

The second term is positive for all $x_{2}\ne 0$ if

(14) $\begin{align*}\frac{p_{22}}{p_{11}}-\frac{p_{12}^{2}}{p_{11}^{2}}>0\end{align*}$

The last equation is equivalent to

(15) $\begin{align*}p_{22}p_{11}-p_{12}^{2}>0\end{align*}$

To summarize, the quadratic form (7) is positive definite if

(16) $\begin{align*}p_{11}>0,\;\;\; p_{22}p_{11}-p_{12}^{2}>0 \end{align*}$

On the other hand, consider again the original matrix $P$

(17) $\begin{align*}P=\begin{bmatrix}p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix}\end{align*}$

The number $p_{11}$ is the leading principal minor of order $1$ (or briefly, the first-order leading principal minor) of the matrix $P$ and the expression $p_{22}p_{11}-p_{12}^{2}$ is the leading principal minor of order $2$ (or briefly, the second-order leading principal minor) of the matrix $P$ . The leading principal minors of certain order are determinants of the upper-left submatrices of the matrix. Generally speaking, we can say that the leading principal minor of order $k$ is the determinant of a matrix obtained by erasing the last $n-k$ columns and $n-k$ rows of the matrix $P\in \mathbb{R}^{n\times n}$ . Thus, the second-order matrix $P$ is positive definite if and only if, its first-order and second-order leading principal minors are positive.

By using this test, we can determine that the matrices $P_{1}$ and $P_{2}$ defined in (8) are positive. Consequently, the corresponding quadratic forms that are shown in Fig. 1. and Fig. 2. are positive definite.

Condition for positive definiteness:
Generally speaking, the n-dimensional matrix $P\in \mathbb{R}^{n\times n}$ :

(18) $\begin{align*}P=\begin{bmatrix}p_{11} & p_{12} & p_{13} & \ldots & p_{1n} \\p_{21} & p_{22} & p_{23} & \ldots & p_{2n} \\ p_{31} & p_{32} & p_{33} & \ldots & p_{3n} \\\vdots & \vdots & \vdots & \ldots & \vdots \\p_{n1} & p_{n2} & p_{n3} & \ldots & p_{n3n} \end{bmatrix}\end{align*}$

is positive definite if and only if

the first-order leading principal minor is positive, that is, the determinant of the upper left 1-by-1 corner of P is positive:
(19) $\begin{align*}D_{1}=p_{11},\;\;\; \text{det}\big(D_{1} \big)=\text{det}(p_{11})>0\end{align*}$
the second-order leading principal minor is positive, that is, the determinant of the upper left 2-by-2 corner of P is positive:
(20) $\begin{align*}D_{2}=\begin{bmatrix}p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix},\;\;\;\text{det}\big(D_{2} \big)>0\end{align*}$
the third-order leading principal minor is positive, that is, the determinant of the upper left 3-by-3 corner of P is positive:
(21) $\begin{align*}D_{3}=\begin{bmatrix}p_{11} & p_{12} & p_{13}\\ p_{21} & p_{22} & p_{23} \\ p_{31} & p_{32} & p_{33} \\\end{bmatrix},\;\;\;\text{det}\big(D_{3} \big)>0\end{align*}$
…
the determinant of the matrix P is positive (this is the $n$ -th order leading principal minor)
(22) $\begin{align*}\text{det}\big(P \big)>0\end{align*}$

This is the Sylvester condition for testing the positive definiteness of a matrix.

Let us apply this test on the following matrix

(23) $\begin{align*}P=\begin{bmatrix} 3 & -1 & 0\\ -1 & 2 & -1 \\ 0 & -1 & 3 \end{bmatrix}\end{align*}$

We have

(24) $\begin{align*}& \text{det}\big(D_{1} \big)=3>0 \\& \text{det}\big(D_{2} \big)=5>0 \\& \text{det}\big(D_{3} \big)=12>0 \\ \end{align*}$

Consequently, the matrix $P$ is positive-definite.

Let us now state definitions of negative definite, negative semidefinite, positive semidefinite, and indefinite matrices.

Negative definite matrix: A matrix $P$ is negative definite if and only if

(25) $\begin{align*}\mathbf{x}^{T}P\mathbf{x}<0, \;\; \forall \mathbf{x}\in \mathbb{R}^{n} \setminus \{0\}\end{align*}$

Negative semi-definite matrix: A matrix $P$ is negative definite if and only if

(26) $\begin{align*}\mathbf{x}^{T}P\mathbf{x}\le 0, \;\; \forall \mathbf{x}\in \mathbb{R}^{n} \end{align*}$

Positive semi-definite matrix: A matrix $P$ is negative definite if and only if

(27) $\begin{align*}\mathbf{x}^{T}P\mathbf{x}\ge 0, \;\; \forall \mathbf{x}\in \mathbb{R}^{n} \end{align*}$

Indefinite matrix: A matrix $P$ is indefinite if and only if $\mathbf{x}^{T}P\mathbf{x}>0$ for some $\mathbf{x}$ , and $\mathbf{x}^{T}P\mathbf{x}<0$ for some $\mathbf{x}$ .

The matrix $P_{3}$ in (8) is negative definite and the matrix $P_{4}$ is indefinite. The quadratic form defined for the negative definite matrix $P_{3}$ is shown in Fig. 3, whereas the quadratic form defined for the indefinite matrix $P_{4}$ is shown in Fig. 4. The matrix $P_{5}$ in (8) is positive semi-definite. The quadratic form defined for this matrix is shown in Fig. 5. The quadratic form for the matrix $P_{5}$ is

(28) $\begin{align*}V(\mathbf{x})=x_{1}^{2}+2x_{1}x_{2}+x_{2}^{2}=(x_{1}+x_{2})^{2} \end{align*}$

We see that this quadratic form is positive semi-definite since for $x_{1}=-x_{2}$ , the quadratic form is zero, and otherwise it is positive (except at $(0,0)$ ). That is, there is some point $(x_{1},x_{2})$ that is different from $(0,0)$ for which this quadratic form is zero. Next, we state tests for negative definiteness, semi-definiteness, and indefiniteness.

Condition for negative definiteness: A matrix $P$ in (18) is negative definite if and only if:

(29) $\begin{align*}\text{det}\big(D_{1}\big)<0 \\\text{det}\big(D_{2}\big)>0 \\\text{det}\big(D_{3}\big)<0 \\\ldots\end{align*}$

That is, the matrix $P$ is negative definite if and only if its leading principal minors alternate the sign, starting from the first-order minor that should be negative. Or in other words, the matrix $P$ is negative definite if and only if its leading principal minors of odd degree are negative and leading principal minors of even degree are positive.

Condition for indefiniteness: A matrix $P$ is indefinite if and only if some leading principal minor is non-zero, but its sign does not follow the pattern for either positive definite or negative definite matrices.

The conditions for positive and negative semidefiniteness are different from the conditions for strict positive and negative definiteness since they involve all principal minors, and not only leading principal minors. Principal minors are simply obtained by erasing certain columns and rows of a matrix and computing the determinant. Here, the process of erasing the columns and rows is not only restricted to the last rows and columns such as in the case of leading principal minors. The column and row numbers should be identical. That is, the indices of deleted rows must be equal to the indices of deleted columns. For example, consider the 3×3 matrix $P$

(30) $\begin{align*}P=\begin{bmatrix}p_{11} & p_{12} & p_{13} \\ p_{21} & p_{22} & p_{23} \\ p_{31} & p_{32} & p_{33} \end{bmatrix}\end{align*}$

Its second order principal minors are

(31) $\begin{align*}D_{21}=\begin{bmatrix}p_{22} & p_{23} \\ p_{32} & p_{33} \end{bmatrix}, D_{22}=\begin{bmatrix}p_{11} & p_{13} \\ p_{31} & p_{33} \end{bmatrix},\; D_{23}=\begin{bmatrix}p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix}\end{align*}$

Condition for positive semi-definiteness: A matrix $P$ is positive semi-definite if and only if its all principal minors are $\ge 0$ .

Condition for negative semi-definiteness: A matrix $P$ is negative semi-definite if and only if its principal minors of odd degree are $\le 0$ , and principal minors of even degree are $\ge 0$ .

Asymptotic stability and the Lyapunov equation

Let us assume that the function (2) is positive definite, and this is equivalent to assuming that the matrix $P$ is positive definite. In the 2D case, as well as in the n-dimensional case, we can always find the parameters $p_{11},p_{12}$ and $p_{22}$ to make the matrix $P$ positive definite. This means that condition (1) of the Lyapunov stability theorem is satisfied.

Let us now compute the derivative of the function $V(\mathbf{x})$ along the trajectories of the system:

(32) $\begin{align*} \dot{V}(\mathbf{x})=\dot{\mathbf{x}}^{T}P\mathbf{x}+\mathbf{x}^{T}P\dot{\mathbf{x}} \end{align*}$

Next, from (1), we have

(33) $\begin{align*}\dot{\mathbf{x}}^{T}=\mathbf{x}^{T}A^{T}\end{align*}$

By substituting (1) and (33) in (32), we obtain

(34) $\begin{align*} \dot{V}(\mathbf{x}) & =\mathbf{x}^{T}A^{T}P\mathbf{x}+\mathbf{x}^{T}PA\mathbf{x} \\& = \mathbf{x}^{T}\Big( A^{T}P+PA\Big)\mathbf{x} \end{align*}$

Let us define the matrix $Q$ :

(35) $\begin{align*}A^{T}P+PA=-Q\end{align*}$

Note since the matrix $P$ is symmetric, the matrix $Q$ is also symmetric. By substituting (35) in (34), we obtain

(36) $\begin{align*} \dot{V}(\mathbf{x}) & =- \mathbf{x}^{T}Q\mathbf{x} \end{align*}$

Since the matrix $Q$ is symmetric, we see that if the matrix $Q$ is positive definite, then the function $\dot{V}(\mathbf{x})$ is negative definite and the equilibrium point is asymptotically stable.

This derivation motivates the introduction of the following theorem.

Theorem about the system stability and the solution of the Lyapunov equation: The real parts of all the eigenvalues of the matrix $A$ of the system (1) are negative (the matrix A is the Hurwitz matrix and the equilibrium point is asymptotically stable), if and only if for any positive definite symmetric matrix $Q$ there exists a positive definite symmetric matrix $P$ that satisfies the Lyapunov equation

(37) $\begin{align*}A^{T}P+PA=-Q\end{align*}$

Furthermore, if the matrix A is Hurwitz, then $P$ is the unique solution of (37).

From everything being said, we can deduce a test for testing the stability of the system:

Stability test of the system:

Select a symmetric positive definite matrix $Q$ (for example, a diagonal matrix with positive entries).
Solve the Lyapunov equation (37)
Test the definiteness of the matrix P. If the matrix $P$ is positive definite, the system is asymptotically stable. If the matrix $P$ is not positive definite, the system is not asymptotically stable.

Basically, this test enables us to investigate the system stability without computing the eigenvalues of the matrix $A$ .

Origins of the Lyapunov Equation in Stability Analysis of Dynamical Systems

Quadratic Forms, Positive Definite, Negative Definite, and Semi-Definite Matrices

Asymptotic stability and the Lyapunov equation

admin

You might also like

Simulation of Transfer Function Response in Python – Control Engineering Tutorials

How to Check Observability of a Dynamical System in Python

Z-Transform Introductory Tutorial with Examples of Z-Transform of Unit Impulse Sequence and Unit Step Sequence