Stability via multiple Lyapunov functions

We start with an example in which we show it can happen that even if the hybrid (switched) system is stable, there is no common quadratic Lyapunov function.

Example 1 (No common quadratic Lyapunov function can be found) We consider a switched system with two modes, linear models of which are parameterized by the matrices \mathbf A_1 and \mathbf A_2:

A₁ = [-0.1 -1; 2 -0.1]
A₂ = [-0.1 -2; 1 -0.1]

Switching curve is given by x_1=0, that is, the vertical axis. State portrait is in Fig. 1.

f₁(x) = A₁*x
f₂(x) = A₂*x
f(x) = x[1] <= 0.0 ? f₁(x) : f₂(x)

using CairoMakie
fig = Figure(size = (800, 800),fontsize=20)
ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
streamplot!(ax,(x₁,x₂)->Point2f(f([x₁,x₂])), -2.0..2.0, -2.0..2.0, colormap = :magma)
vlines!(ax,0.0,ymin=-2.0,ymax=2.0, color = :red, linewidth=3)
fig

Figure 1: State portrait for the switched system with no common quadratic Lyapunov function

The equilibrium (the origin) of this switched system appears stable.

The individual systems are stable, which we can immediately see by computing the eigenvalues of the matrices A_1 and A_2:

using LinearAlgebra
eigvals(A₁), eigvals(A₂)

(ComplexF64[-0.1 - 1.4142135623730951im, -0.1 + 1.4142135623730951im], ComplexF64[-0.1 - 1.4142135623730951im, -0.1 + 1.4142135623730951im])

We now try to find a common quadratic Lyapunov function for both subsystems. We will formulate the problem as an LMI feasibility problem.

using Convex, SCS
X = Semidefinite(2)
constraint₁ = A₁'*X + X*A₁ ⪯ -Matrix{Float64}(I, 2, 2)
constraint₂ = A₂'*X + X*A₂ ⪯ -Matrix{Float64}(I, 2, 2)
constraints = [constraint₁, constraint₂]
problem = satisfy(constraints)
solve!(problem,SCS.Optimizer,silent=true)

Problem statistics
  problem is DCP         : true
  number of variables    : 1 (4 scalar elements)
  number of constraints  : 3 (12 scalar elements)
  number of coefficients : 24
  number of atoms        : 10

Solution summary
  termination status : INFEASIBLE
  primal status      : INFEASIBLE_POINT
  dual status        : INFEASIBILITY_CERTIFICATE

Expression graph
  satisfy
   └─ nothing
  subject to
   ├─ PSD constraint (convex)
   │  └─ + (affine; real)
   │     ├─ 2×2 Matrix{Float64}
   │     └─ Convex.NegateAtom (affine; real)
   │        └─ …
   ├─ PSD constraint (convex)
   │  └─ + (affine; real)
   │     ├─ 2×2 Matrix{Float64}
   │     └─ Convex.NegateAtom (affine; real)
   │        └─ …
   ├─ PSD constraint (convex)
   │  └─ 2×2 real variable (id: 636…095)
   ⋮

The solver does not find a solution. Well, perhaps trying another solver or two would make our conclusion more robust (MosekTools is another recommendable alternative to SCS). But we are now be tempted to conclude that there is no common quadratic Lyapunov function for both subsystems.

Dual LMI problem confirms the infeasibility of the primal one

It is possible to formulate a dual LMI problem. Its feasibility certifies infeasibility of the primal one. Namely, the problem is, if there are matrices \mathbf R_1\succ 0 and \mathbf R_2\succ 0 such that the following LMI holds: \mathbf R_1\mathbf A_1^\top+\mathbf A_1\mathbf R_1 + \mathbf R_2\mathbf A_2^\top+\mathbf A_2\mathbf R_2 \prec 0.

The conclusion about the imposibility to find a single quadratic Lyapunov function for both subsystems is also suppported by plotting the invariant sets for the two subsystems. First, we need to compute Lyapunov functions for the two subsystems.

X₁ = Semidefinite(2)
constraint₁ = A₁'*X₁ + X₁*A₁ ⪯ -Matrix{Float64}(I, 2, 2)
problem₁ = satisfy(constraint₁)
solve!(problem₁,SCS.Optimizer,silent=true)
X₁.value

2×2 Matrix{Float64}:
 423.71     12.5739
  12.5739  212.123

X₂ = Semidefinite(2)
constraint₂ = A₂'*X₂ + X₂*A₂ ⪯ -Matrix{Float64}(I, 2, 2)
problem₂ = satisfy(constraint₂)
solve!(problem₂,SCS.Optimizer,silent=true)
X₂.value

2×2 Matrix{Float64}:
 212.123   -12.5739
 -12.5739  423.71

Generally, a Lyapunov function has the property that its sublevel set \{\bm x \mid V(\bm x) \leq \alpha\} is forward (also positive) invariant. We plot invariant sets corresponding to some particular value for both subsystems superposed on the state portrait in Fig. 2.

x1s = LinRange(-2, 2, 100)
x2s = LinRange(-2, 2, 100)
V₁(x) = x'*X₁.value*x
V₂(x) = x'*X₂.value*x
V1s = [V₁([x₁,x₂]) for x₁ in x1s, x₂ in x2s]
V2s = [V₂([x₁,x₂]) for x₁ in x1s, x₂ in x2s]
contour!(x1s, x2s, V1s, levels=[300.0], linewidth=3, color=:blue)
contour!(x1s, x2s, V2s, levels=[300.0], linewidth=3, color=:green)
fig

Figure 2: Invariant ellipses for the two subsystems superposed on the state portrait

We can see in Fig. 2 that none of the two ellipses works as an invariant set for the switched system – a state trajectory entering the set leaves it afterwards. No way to to come up with a single ellipse (hence a single quadratic Lyapunov function) that would work here.

Good, we have seen in the example that it is not always possible to find a common quadratic Lyapunov function for a switched system, even it it is stable. We need to expand the set of functions in which we search for a Lyapunov function. We have proposed one way to do it in the previouse chapter wherein we considered higher degree polynomials on which we imposed the nonnegativity constraint through in the form of SOS constraint. Here we are going to consider another approach. We are goint to stitch together several Lyapunov-like functions, each of which is a Lyapunov function only on some subset of the state space (that is why we call them just Lyapunov-like and not Lyapunov). This approach is sometimes called the Multiple Lyapunov Function (MLF) approach, or Piecewise Lyapunov Function approach.

Multiple Lyapunov Function (MLF) approach to analysis of stability

Instead of just a single common Lyapunov function V(\bm x), we are now going to consider several Lyapunov-like functions V_i(\bm x),\; i=1,\ldots,r, that qualify as Lyapunov function only on some subsets of the state space \Omega_i. And we “stitch” them together to form a piecewise Lyapunov function V(\bm x): V(x) = \begin{cases} V_1(\bm x) & \text{if } \bm x\in \Omega_1, \\ \vdots\\ V_r(\bm x) & \text{if } \bm x\in \Omega_r. \end{cases}

S-procedure

In order to restrict the requirement of positive definiteness of the Lyapunov function to some region in the state space, and similarly for the regurement of negative definiteness of its time derivative, we need to introduce the the S-procedure. This is a result about solvability of two or more quadratic inequalities, not necessarily convex ones (for convex problems we have nonlinear Farkas’ lemma). Origins of this result can be found in the control theory (analysis of stability of nonlinear systems, hence the S letter) with the first rigorous result provided by Yakubovich in 1970s.

It gives conditions under which (satisfaction of) one quadratic inequality follows from (satisfaction of) another one (or more). Namely, it gives a condition under which the following implication holds:

\boxed {\text{Quadratic inequality \#1 satisfied by some}\; \bm x \Rightarrow \text{Quadratic inequality \#0 satisfied by the same}\; \bm x.}

In other words, it gives a condition under which the solution set of the inequality #1 denoted as \mathcal X_1 is included in the solution set \mathcal X_0 of the inequality #0.

S-procedure with nonstrict inequalities

Consider the general quadratic functions F_i(\bm x) = \bm x^\top \mathbf A_i \bm x + 2\mathbf b_i^\top \bm x + c_i, \; i=0,\ldots, p.

The question is: under which conditions it holds that F_0(\bm x) \geq 0 for all \bm x satisfying F_i(\bm x)\geq 0,\; i=1,\ldots,p ?

In other words, we are looking for conditions under which the implication F_i(\bm x) \geq 0,\; i=1,\ldots,p \quad \Rightarrow \quad F_0(\bm x) \geq 0
holds.

In the simplest (yet relevant) case p=1 we search for conditions under which F_0(\bm x) \geq 0 for all \bm x satisfying F_1(\bm x)\geq 0, that is, conditions under which the implication F_1(\bm x) \geq 0 \quad \Rightarrow \quad F_0(\bm x) \geq 0
holds.

Sufficient conditions

The existence of \alpha_i\geq 0,\; i=1,\ldots,p such that F_0(\bm x)-\sum_{i=1}^p \alpha_i F_i(\bm x) \geq 0 is sufficient for the original implication to hold. Generally, it is not necessary; the condition is conservative.

It can be formulated as an LMI \begin{bmatrix} \mathbf A_0 & \mathbf b_0 \\ \mathbf b_0^\top & c_0 \end{bmatrix} - \sum_{i=1}^p \alpha_i \begin{bmatrix} \mathbf A_i & \mathbf b_i \\ \mathbf b_i^\top & c_i \end{bmatrix} \succeq 0
where \alpha_i \geq 0,\; i=1,\ldots,p.

Sufficient and also necessary

It is nontrivial that for p=1 it is also necessary, provided that there is some \bm x_0 such that F_1(\bm x_0)>0. Then we have the following equivalence between the two constraints: \begin{aligned} F_0(\bm x) &\geq 0 \; \forall \bm x \;\mathrm{satisfying}\; F_1(\bm x)\geq 0 \\ &\Longleftrightarrow \\ F_0(\bm x)-\alpha F_1(\bm x) &\geq 0,\;\text{for some}\; \alpha\in\mathbb{R}, \; \alpha\geq 0, \end{aligned}
which again can be formulated as an LMI, namely
\begin{bmatrix} \mathbf A_0 & \mathbf b_0 \\ \mathbf b_0^\top & c_0 \end{bmatrix} - \alpha \begin{bmatrix} \mathbf A_1 & \mathbf b_1 \\ \mathbf b_1^\top & c_1 \end{bmatrix} \succeq 0,\quad \alpha\geq 0.

More on S-procedure

There are several variants

• strict, nonstrict or mixed inequalities,
• just two or more,
• some of the constraints can be equations.

Piecewise quadratic Lyapunov function

We now restrict ourselves to quadratic Lyapunov-like functions, that is, quadratic functions V_i(\bm x) = \bm x^\top \mathbf P_i \bm x that qualify as Lyapunov functions only on respective subsets \Omega_i\sub \mathbb R^n:

V_i(\bm x) = \bm x^\top \mathbf P_i \bm x > 0\quad \forall \;\bm x\in \Omega_i,

\dot V_i(\bm x) = \bm x^\top \left( \mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i \right) \bm x < 0\quad \forall \;\bm x\in \Omega_i.

Using comparison functions and nonstrict inequalities

We can use our good old comparison functions to formulate the conditions of positive definiteness and negative definiteness. \alpha_1 \bm x^\top \mathbf I \bm x \leq \bm x^\top \mathbf P_i \bm x \leq \alpha_2 \bm x^\top \mathbf I \bm x \quad \forall \;\bm x\in \Omega_i,

\bm x^\top \left( \mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i \right) \bm x \leq -\alpha_3 \bm x^\top \mathbf I \bm x\quad \forall \;\bm x\in \Omega_i.

The difference now is that these conditions are only reuired to hold on some state regions, some subsets of the state space. It is now time to discuss how to characterize those regions.

Characterization of subsets of state space using LMI

Some subsets \Omega_i\sub \mathbb R^n characterized using linear and quadratic inequalities can be formulated within the LMI framework as \bm x^\top \mathbf Q_i \bm x \geq 0.

In particular, centered ellipsoids and cones.

For example,
\begin{aligned} \Omega_i &= \{\bm x \in \mathbb R^n \mid (\mathbf c^\top \bm x \geq 0 \land \mathbf d^\top \bm x \geq 0) \\ & \qquad \qquad \qquad \lor (\mathbf c^\top \bm x \leq 0 \land \mathbf d^\top \bm x \leq 0)\}. \end{aligned}

This constraint can be reformulated as (\mathbf c^\top \bm x) (\mathbf d^\top \bm x) \geq 0,
which can be reformatted to \bm x^\top \mathbf c \mathbf d^\top \bm x \geq 0,
which can further be symmetrized to \bm x^\top \left(\frac{\mathbf c \mathbf d^\top + \mathbf d \mathbf c^\top}{2}\right) \bm x \geq 0.

More general sets (general polyhedra, noncentered ellipsoids) can also be modelled using LMI too… We are going to have a look at them, but first we hurry to show how to combine the subset characterization and Lyapunov-ness using the S-procedure.

Combining the subset characterization and Lyapunov-ness using the S-procedure

We want to learn if the the following hold \alpha_i \bm x^\top \mathbf I \bm x \leq \bm x^\top \mathbf P_i \bm x,

\bm x^\top \left( \mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i \right) \bm x \leq -\gamma_i \bm x^\top \mathbf I \bm x,

but not for all \bm x, but only for \bm x\in \Omega_i, that is, all \bm x satisfying \bm x^\top \mathbf Q_i \bm x \geq 0. But this is now a perfect opportunity for application of the S-procedure:

\mathbf P_i - \alpha_i \mathbf I - \mu_i \mathbf Q_i \succeq 0,\quad \mu_i \geq 0,\; \alpha_i > 0,

\mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i + \gamma_i \mathbf I + \xi_i \mathbf Q \preceq 0,\quad \mu_i \geq 0,\; \gamma_i > 0.

More general sets using LMI

How can we model more general sets using LMI?

The inequality \bm x^\top \mathbf Q \bm x + 2\mathbf r^\top \bm x + s \geq 0,

can be reformulated as \begin{bmatrix} \bm x^\top & 1 \end{bmatrix} \underbrace{ \begin{bmatrix} \mathbf Q & \mathbf r \\ \mathbf r^\top & s \end{bmatrix}}_{\bar{\mathbf{Q}}} \underbrace{ \begin{bmatrix} \bm x \\ 1 \end{bmatrix}}_{\bar{\bm x}} \geq 0,

that is, as an LMI \begin{bmatrix} \mathbf Q & \mathbf r \\ \mathbf r^\top & s \end{bmatrix} \succeq 0.

Affine subspace

\mathbf c^\top \bm x + d \geq 0,

\begin{bmatrix} \bm x^\top & 1 \end{bmatrix} \begin{bmatrix} \mathbf 0 & \mathbf c \\ \mathbf c^\top & 2d \end{bmatrix} \begin{bmatrix} \bm x \\ 1 \end{bmatrix} \geq 0,

\begin{bmatrix} \mathbf 0 & \mathbf c \\ \mathbf c^\top & 2d \end{bmatrix} \succeq 0.

But then the Lyapunov-like functions and system matrices must also be extended V(\bm x) = \begin{bmatrix} \bm x^\top & 1 \end{bmatrix} \underbrace{ \begin{bmatrix} \mathbf P & \mathbf P_{12} \\ \mathbf P_{12} & P_{22} \end{bmatrix}}_{\bar{\mathbf P}} \begin{bmatrix} \bm x \\ 1 \end{bmatrix},

\bar{\mathbf{A}} = \begin{bmatrix} \mathbf A & \mathbf 0 \\ \mathbf 0 & 0 \end{bmatrix}.

Continuity conditions

The boundary between the regions \Omega_i and \Omega_j described by \Omega_{ij} = \{\bm x \in \mathbb R^n \mid \mathbf F_{ij} \bm z + \mathbf{l}_{ij}\}, where \bm z\in \mathbb R^p, \mathbf F_{ij}\in \mathbb R^{n\times p}, and \mathbf l_{ij}\in \mathbb R^{n}.

The continuity conditions are V_i(\bm x) = V_j(\bm x) \quad \forall \bm x \in \Omega_{ij}, which can be reformulated as \begin{aligned} &\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \mathbf P_i \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right) \\ & \qquad + 2\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \bm q_i + r_i \\ &= \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \mathbf P_j \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right) \\ & \qquad + 2\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \bm q_j + r_j \end{aligned}.

Collecting the terms with equal powers of \bm z. \begin{aligned} \mathbf F_{ij}^\top (\mathbf P_1 - \mathbf P_2) \mathbf F_{ij} &= 0, \\ \mathbf F_{ij}^\top (\mathbf P_1 - \mathbf P_2) \mathbf l_{ij} + (\mathbf q_1-\mathbf q_2) &= 0, \\ \mathbf l_{ij}^\top (\mathbf P_1 - \mathbf P_2)\mathbf l_{ij} + 2\mathbf l_{ij}^\top (\mathbf q_1-\mathbf q_2) + r_1-r_2 &= 0. \end{aligned}