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 that 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: 114…880)
   ⋮

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 inclined to conclude that there is no common quadratic Lyapunov function for both subsystems.

NoteDual LMI problem confirms the infeasibility of the primal one

It is possible to certify infeasibility of the primal optimization problem by certifying feasibility of the dual one. Namely, if there are matrices \bm R_1\succ 0 and \bm R_2\succ 0 such that the following LMI holds: \bm R_1\mathbf A_1^\top+\mathbf A_1\bm R_1 + \bm R_2\mathbf A_2^\top+\mathbf A_2\bm R_2 \prec 0, the primal problem is infeasible.

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 respective 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 previous chapter, wherein we considered higher degree polynomials for which we imposed the nonnegativity constraint through the SOS constraint. Here we are going to consider another approach. We are going 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 individually 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}

Before we proceed, we need the following technical results.

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 \bm P_i \bm x that qualify as Lyapunov functions only on respective subsets \Omega_i\sub \mathbb R^n:

\begin{aligned} V_i(\bm x) = \bm x^\top \bm 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 \bm P_i + \bm P_i \mathbf A_i \right) \bm x &< 0\quad \forall \;\bm x\in \Omega_i. \end{aligned}

Using comparison functions and nonstrict inequalities

We can use our good old comparison functions to formulate the conditions of positive definiteness and negative definiteness. \begin{aligned} \alpha_i \bm x^\top \mathbf I \bm x &\leq \bm x^\top \bm P_i \bm x \;{\color{grey}(\leq \beta_i \bm x^\top \mathbf I \bm x)} \quad \forall \;\bm x\in \Omega_i,\\ \bm x^\top \left( \mathbf A_i^\top \bm P_i + \bm P_i \mathbf A_i \right) \bm x &\leq -\gamma_i \bm x^\top \mathbf I \bm x \quad \forall \;\bm x\in \Omega_i. \end{aligned}

The difference now is that these conditions are only required 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,
\Omega_i = \{\bm x \in \mathbb R^n \mid (\mathbf c^\top \bm x \geq 0 \land \mathbf d^\top \bm x \geq 0) \lor (\mathbf c^\top \bm x \leq 0 \land \mathbf d^\top \bm x \leq 0)\}.

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 test for asymptotic stability, that is, we want to learn if the the following hold \begin{aligned} \alpha_i \bm x^\top \mathbf I \bm x &\leq \bm x^\top \bm P_i \bm x,\\ \bm x^\top \left( \mathbf A_i^\top \bm P_i + \bm P_i \mathbf A_i \right) \bm x &\leq -\gamma_i \bm x^\top \mathbf I \bm x, \end{aligned}

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:

\boxed {\begin{aligned} \bm P_i - \alpha_i \mathbf I - \mu_i \mathbf Q_i &\succeq 0,\quad \alpha_i > 0, \; \mu_i \geq 0,\\ \mathbf A_i^\top \bm P_i + \bm P_i \mathbf A_i + \gamma_i \mathbf I + \xi_i \mathbf Q_i &\preceq 0,\quad \gamma_i > 0, \; \xi_i \geq 0. \end{aligned}}

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 c^\top \bm x + d \geq 0,

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

that is, as an LMI with the matrix variable \bar{\mathbf{Q}}, that is, \begin{bmatrix} \mathbf Q & \mathbf c \\ \mathbf c^\top & d \end{bmatrix} \succeq 0.

Affine subspace

A special instance is that of an affine subspace characterized by \mathbf c^\top \bm x + \frac{d}{2} \geq 0, which we can reformulate as \begin{bmatrix} \bm x^\top & 1 \end{bmatrix} \begin{bmatrix} \mathbf 0 & \mathbf c \\ \mathbf c^\top & d \end{bmatrix} \begin{bmatrix} \bm x \\ 1 \end{bmatrix} \geq 0, which can be expressed as an LMI \begin{bmatrix} \mathbf 0 & \mathbf c \\ \mathbf c^\top & d \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} \bm P & \bm q \\ \bm q^\top & r \end{bmatrix}}_{\bar{\bm P}} \begin{bmatrix} \bm x \\ 1 \end{bmatrix}, and the system matrix as well: \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 \bm P_i \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right) + 2\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \bm q_i + r_i \\ &\qquad \qquad = \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \bm P_j \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right) + 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, \boxed {\begin{aligned} \mathbf F_{ij}^\top (\bm P_1 - \bm P_2) \mathbf F_{ij} &= 0, \\ \mathbf F_{ij}^\top (\bm P_1 - \bm P_2) \mathbf l_{ij} + (\mathbf q_1-\mathbf q_2) &= 0, \\ \mathbf l_{ij}^\top (\bm P_1 - \bm P_2)\mathbf l_{ij} + 2\mathbf l_{ij}^\top (\mathbf q_1-\mathbf q_2) + r_1-r_2 &= 0. \end{aligned}}

Examples

Example 2 (Piecewise quadratic Lyapunov function for the flower example) We consider the switched system with two modes, linear models that are parameterized by the matrices \mathbf A_1 and \mathbf A_2, for which their domains are given by the sectors defined by the lines x_2=x_1 and x_2=-x_1. State portrait is in Fig. 3.

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

using Convex, SCS

c = [1, 1]
d = [-1, 1]
Q₁ = (c*d' + d*c')/2
f(x) = dot(x,Q₁,x) >= 0 ? A₁*x : A₂*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, alpha = 0.6)
lines!(ax,[-2,2],[-2,2], color = :red, linewidth=3)
lines!(ax,[-2,2],[2,-2], color = :red, linewidth=3)
fig

Figure 3: State portrait for the switched system with piecewise-defined dynamics

The implementation of the stability analysis based on formulating linear matrix inequalities is as follows

P₁ = Semidefinite(2)
P₂ = Semidefinite(2)
α₁ = 1e-2
α₂ = 1e-2
γ₁ = Variable(Positive())
γ₂ = Variable(Positive())
μ₁ = Variable(Positive())
μ₂ = Variable(Positive())
ξ₁ = Variable(Positive())
ξ₂ = Variable(Positive())

Q₂ = -Q₁

using LinearAlgebra
constraint₁ = P₁ - α₁*Matrix{Float64}(I, 2, 2) - μ₁*Q₁ ⪰ 0
constraint₂ = A₁'*P₁ + P₁*A₁ + γ₁*Matrix{Float64}(I, 2, 2) + ξ₁*Q₁ ⪯ 0

constraint₃ = P₂ - α₂*Matrix{Float64}(I, 2, 2) - μ₂*Q₂ ⪰ 0
constraint₄ = A₂'*P₂ + P₂*A₂ + γ₂*Matrix{Float64}(I, 2, 2) + ξ₂*Q₂ ⪯ 0

# Continuity condition at the switching surface
constraint₅ = dot(c, P₁-P₂, c) == 0
constraint₆ = dot(d, P₁-P₂, d) == 0
constraints = [constraint₁, constraint₂, constraint₃, constraint₄, constraint₅, constraint₆]

problem = satisfy(constraints)
solve!(problem,SCS.Optimizer,silent=true)

P₁, P₂ = P₁.value, P₂.value

([0.11467296044627187 -1.0040614304026844e-12; -1.1668125205840983e-12 0.023007007749945577], [0.023007007749945296 1.004072143028592e-12; 1.1668232180154956e-12 0.11467296044627218])

Several sublevel sets for the computed piecewise quadratic Lyapunov function are superposed onto the state portrait in Fig. 4.

V(x) = dot(x,Q₁,x) >= 0 ? x'*P₁*x : x'*P₂*x

x₁ = x₂ = range(-2.0, 2.0; length=100)
X₁ = [i for i in x₁, j in x₂]
X₂ = [j for i in x₁, j in x₂]
Z1 = [V([X₁[i,j], X₂[i,j]]) for i in 1:length(x₁), j in 1:length(x₂)]
levels = logrange(1e-3,1e1,10)
contour!(ax, x₁, x₂, Z1; levels=levels, linewidth=2)
fig

Figure 4: State portrait with the piecewise quadratic Lyapunov function’s sublevel sets superposed

The exotic shape of this piecewise quadratic Lyapunov function is worth visualizing as a surface in the three-dimensional space as in Fig. 5.

using CairoMakie
fig = Figure(size = (800, 800),fontsize=20)
ax = Axis3(fig[1, 1], xlabel = "x₁", ylabel = "x₂", zlabel = "V(x)")
x₁ = x₂ = range(-2.0, 2.0; length=100)
X₁ = [i for i in x₁, j in x₂]
X₂ = [j for i in x₁, j in x₂]
Z1 = [V([X₁[i,j], X₂[i,j]]) for i in 1:length(x₁), j in 1:length(x₂)]
surface!(ax, x₁, x₂, Z1; colormap = :magma)
fig

Figure 5: 3D surface plot of the piecewise quadratic Lyapunov function

This piecewise quadratic Lyapunov function certifies asymptotic stability of the switched system.

Example 3 (Piecewise quadratic Lyapunov function for Example 1) We have already found that a common Lyapunov function could not be found among quadratic functions. Here we attempt to find a piecewise quadratic Lyapunov function. Namely, we are going to search for two quadratic Lyapunov-like functions, each of which is valid on one side of the switching curve, which is the vertical axis in this case. And finally we stitch the two together via a continuity condition. The difference with respect to the previous example is that since the two domains are half-plances, we need to embed the system matrices and Lyapunov-like function matrices into matrices of size 3x3, as explained in the text above.

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

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

using Convex, SCS

Ā₁ = zeros(Float64, 3, 3)
Ā₁[1:2, 1:2] = A₁         # A₁ is embedded in the top-left corner of a 3x3 matrix

Ā₂ = zeros(Float64, 3, 3)
Ā₂[1:2, 1:2] = A₂         # A₂ is embedded in the top-left corner of a 3x3 matrix

P̄₁ = Semidefinite(3)
P̄₂ = Semidefinite(3)
α₁ = 1e-2
α₂ = 1e-2
γ₁ = Variable(Positive())
γ₂ = Variable(Positive())
μ₁ = Variable(Positive())
μ₂ = Variable(Positive())
ξ₁ = Variable(Positive())
ξ₂ = Variable(Positive())

Q₁ = zeros(Float64, 3, 3) # Region x₁ ≤ 0  =>  -x₁ ≥ 0
Q₁[1,3] = -1.0
Q₁[3,1] = -1.0

Q₂ = -Q₁  # Region x₁ ≥ 0

using LinearAlgebra
constraint₁ = P̄₁ - α₁*Matrix{Float64}(I, 3, 3) - μ₁*Q₁ ⪰ 0
constraint₂ = Ā₁'*P̄₁ + P̄₁*Ā₁ + γ₁*Matrix{Float64}(I, 3, 3) + ξ₁*Q₁ ⪯ 0

constraint₃ = P̄₂ - α₂*Matrix{Float64}(I, 3, 3) - μ₂*Q₂ ⪰ 0
constraint₄ = Ā₂'*P̄₂ + P̄₂*Ā₂ + γ₂*Matrix{Float64}(I, 3, 3) + ξ₂*Q₂ ⪯ 0

# Continuity condition at the switching surface x₁ = 0
constraint₅ = P̄₁[2,2] == P̄₂[2,2]
constraints = [constraint₁, constraint₂, constraint₃, constraint₄, constraint₅]

problem = satisfy(constraints)
solve!(problem,SCS.Optimizer,silent=true)

P₁, P₂ = P̄₁.value[1:2, 1:2], P̄₂.value[1:2, 1:2]

([0.04385265895052883 0.00031044970469989373; 0.0003104497560347673 0.023985039700069443], [0.011614022462915962 -0.0005879313889673506; -0.00058794336958938 0.023985040899855464])

Several sublevel sets for the computed piecewise quadratic Lyapunov function are superposed onto the state portrait in Fig. 6.

V(x) = x[1] <= 0.0 ? x'*P₁*x : x'*P₂*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, alpha = 0.6)
vlines!(ax,0.0,ymin=-2.0,ymax=2.0, color = :red, linewidth=3)
fig 
x₁ = x₂ = range(-2.0, 2.0; length=100)
X₁ = [i for i in x₁, j in x₂]
X₂ = [j for i in x₁, j in x₂]
Z1 = [V([X₁[i,j], X₂[i,j]]) for i in 1:length(x₁), j in 1:length(x₂)]
levels = logrange(1e-2,1e1,10)
contour!(ax, x₁, x₂, Z1; levels=levels, linewidth=2)
fig

Figure 6: State portrait with the piecewise quadratic Lyapunov function’s sublevel sets superposed

Example 4 (Higher-degree polynomial Lyapunov function for Example 1) In this example we demonstrate that it is possible to combine the results from this chapter with those of the previous one. Namely, instead of searching for several quadratic Lyapunov-like functions that are then stitched together to form a piecewise quadratic Lyapunov function, we can also try to find a single higher-degree SOS polynomial function that qualifies as a Lyapunov function for both subsystems only on their respective domains.

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

f(x) = x[1] <= 0.0 ? A₁*x : A₂*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, alpha = 0.6)
vlines!(ax,0.0,ymin=-2.0,ymax=2.0, color = :red, linewidth=3)
fig

using DynamicPolynomials
using SemialgebraicSets
using SumOfSquares
using CSDP   

solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
@polyvar x[1:2]

d = 8
monos = monomials(x,2:d)
@variable(model, V_sos, SumOfSquares.Poly(monos))

using LinearAlgebra # Needed for `dot`

dV1dt_sos = dot(differentiate(V_sos, x), A₁*x)
dV2dt_sos = dot(differentiate(V_sos, x), A₂*x)

ϵ = 1e-3
φ = dot(monomials(x,0:(d÷2)),monomials(x,0:(d÷2)))

Q₁ = @set x[1] <= 0.0 
Q₂ = @set x[1] >= 0.0

@constraint(model, V_sos >= ϵ*φ)
@constraint(model, dV1dt_sos <= -ϵ*φ, domain = Q₁)
@constraint(model, dV2dt_sos <= -ϵ*φ, domain = Q₂)

JuMP.optimize!(model)
JuMP.primal_status(model)
L_sos = value(V_sos)

0.003952682168346655x₂² + 0.0009619573720556218x₁x₂ + 0.005052837664152321x₁² - 0.0008835307071421994x₂³ - 0.00428159348484769x₁x₂² - 0.0016380519282392925x₁²x₂ - 0.018771479943097802x₁³ + 0.022288429163381807x₂⁴ + 0.0023429662296621245x₁x₂³ + 0.05045960880124767x₁²x₂² + 0.0016680790613463614x₁³x₂ + 0.04185683650393912x₁⁴ - 0.0016599901209701784x₂⁵ - 0.002075666364362405x₁x₂⁴ - 0.007522912513195479x₁²x₂³ - 0.03277394879842177x₁³x₂² - 0.002694867373975285x₁⁴x₂ - 0.023748388578496815x₁⁵ + 0.001783181311111548x₂⁶ + 0.00039192032727441983x₁x₂⁵ + 0.00506562890677742x₁²x₂⁴ + 0.0025091236375374137x₁³x₂³ + 0.010926253350589832x₁⁴x₂² + 0.0009061110440597986x₁⁵x₂ + 0.005737503322961857x₁⁶ + 5.890828470000997e-5x₂⁷ - 1.2867503755842336e-5x₁x₂⁶ + 0.0001809388058973127x₁²x₂⁵ - 2.0811858121305704e-6x₁³x₂⁴ + 0.00016238967873505317x₁⁴x₂³ - 0.0005696839834854472x₁⁵x₂² + 1.7584510715096258e-5x₁⁶x₂ - 0.0003698326527228346x₁⁷ + 1.2178697033959907e-5x₂⁸ - 0.00010888981887546834x₁x₂⁷ + 4.843364422413288e-5x₁²x₂⁶ - 4.735463699034881e-5x₁³x₂⁵ + 5.027391762268962e-5x₁⁴x₂⁴ + 6.256856977415737e-5x₁⁵x₂³ + 1.6007582416932564e-5x₁⁶x₂² + 0.00010169200686505064x₁⁷x₂ + 2.0688508811872452e-5x₁⁸

x₁ = x₂ = range(-2.0, 2.0; length=100)
X₁ = [i for i in x₁, j in x₂]
X₂ = [j for i in x₁, j in x₂]

Z1 = [L_sos([X₁[i,j], X₂[i,j]]) for i in 1:length(x₁), j in 1:length(x₂)]

levels = logrange(1e-2,1e1,10)
contour!(ax, x₁, x₂, Z1; levels=levels, linewidth=2)
fig

Figure 7: State portrait with the higher-degree SOS polynomial Lyapunov function’s sublevel sets superposed