Lyapunov Differential Equation Hierarchy and Polynomial Lyapunov Functions for Switched Linear Systems

TL;DR

This paper introduces a hierarchy of Lyapunov differential equations that simplifies the search for polynomial Lyapunov functions in stable switched linear systems, offering a new computational approach and comparison with traditional methods.

Contribution

It establishes an equivalence between polynomial and quadratic Lyapunov functions via a hierarchy of differential equations, enabling efficient stability analysis and high-order polynomial function generation.

Findings

The proposed method effectively generates high-order polynomial Lyapunov functions.

It provides a computationally competitive alternative to sum-of-squares approaches.

The approach offers an intuitive procedure for stability verification of switched linear systems.

Abstract

This work studies the problem of searching for homogeneous polynomial Lyapunov functions for stable switched linear systems. Specifically, we show an equivalence between polynomial Lyapunov functions for systems of this class and quadratic Lyapunov functions for a related hierarchy of Lyapunov differential equations. This creates an intuitive procedure for checking the stability properties of switched linear systems and a computationally competitive algorithm is presented for generating high-order homogeneous polynomial Lyapunov functions in this manner. Additionally, we provide a comparison between polynomial Lyapunov functions generated with our proposed approach and polynomial Lyapunov functions generated with a more traditional sum-of-squares based approach.

Equations

$$\dot{x}=A(t)x,$$

$$\begin{array}[]{rcl}V(x)&>&0\\ \dot{V}(x)=\langle\nabla V,A_{i}x\rangle&<&0\\ \forall x&\neq&0\\ \forall i&\in&\{1,\,\cdots,\,N\}.\end{array}$$

$$V(x)=x^{T}Px$$

$$A_i^{T}P+P{A}_i<0$$

$$A_{i}Q+Q{A}_i^T<0$$

$$V(x)=x^T{Q}^{-1}x$$

$$\dot{X}=A(t)X+XA(t{)}^T,$$

$$\begin{array}{rl}\dot{X}& =\dot{x}{x}^T+x{\dot{x}}^T\\ & =A(t)x{x}^T+x{x}^{T}A(t{)}^T\\ & =A(t)X+XA(t{)}^T.\end{array}$$

$$X_0=\sum _{j=1}^N{p}_{j,0}{q}_{j,0}^T.$$

$$\frac{d}{dt}{p}_j(t)=A(t)p_j(t),p_j(0)=p_{j,0},$$

$$\frac{d}{dt}{q}_j(t)=A(t)q_j(t),q_j(0)=q_{j,0}.$$

$$\begin{split}\dot{X}(t)&=\sum_{j=1}^{N}\Big{(}\dot{p}_{j}q_{j}^{T}+p_{j}\dot{q}_{j}^{T}\Big{)}\\ &=A(t)\bigg{(}\sum_{j=1}^{N}p_{j}q_{j}^{T}\bigg{)}+\bigg{(}\sum_{j=1}^{N}p_{j}q_{j}^{T}\bigg{)}A(t)^{T}\\ &=A(t)X(t)+X(t)A(t)^{T}.\end{split}$$

$$\begin{array}[]{rllc}\otimes^{1}x&=x&\in\mathbb{R}^{n},&\\ \otimes^{k}x&=x\otimes(\otimes^{k-1}x)&\in\mathbb{R}^{n^{k}},&k\geq 2.\end{array}$$

$$\dot{\vec{X}}=\mathcal{A}(t)\vec{X}$$

$${\mathcal{A}}_i:=I_n\otimes A_i+A_i\otimes I_n$$

$${\mathcal{A}}_i^{T}P+P{\mathcal{A}}_i<0,$$

$$\mathcal{V}(\vec{X}):={\vec{X}}^{T}P\vec{X}>0$$

$$\dot{\mathcal{V}}(\vec{X})\;=\;\vec{X}^{T}\Big{(}\mathcal{A}(t)^{T}P+P\mathcal{A}(t)\Big{)}\vec{X}.$$

$$\begin{array}{rl}\mathcal{A}(t{)}^{T}P& =(I_n\otimes A(t{)}^T+A(t{)}^T\otimes I_n)Q\otimes Q\\ & =Q\otimes (A(t{)}^{T}Q)+(A(t{)}^{T}Q)\otimes Q\end{array}$$

$$\begin{array}{rl}P\mathcal{A}(t)& =Q\otimes Q(I_n\otimes A(t)+A(t)\otimes I_n)\\ & =Q\otimes (QA(t))+(QA(t))\otimes Q.\end{array}$$

$$\begin{gathered} \mathcal{A}(t{)}^{T}P+P\mathcal{A}(t)=Q\otimes (A(t{)}^{T}Q+QA(t))+\cdots \\ +(A(t{)}^{T}Q+QA(t))\otimes Q. \end{gathered}$$

$$\vec{X}^{T}\Big{(}Q\otimes(A(t)^{T}Q+QA(t))\Big{)}\vec{X}=\\ \begin{array}[]{l}\cdots=\vec{X}^{T}\mbox{vec}\Big{(}(A(t)^{T}Q+QA(t))XQ\Big{)}\\ \cdots=\mbox{trace}\Big{(}X^{T}(A(t)^{T}Q+QA(t))XQ\Big{)}\\ \cdots=\mbox{trace}\Big{(}Q^{1/2}X^{T}(A(t)^{T}Q+QA(t))XQ^{1/2}\Big{)},\end{array}$$

$$\vec{X}^{T}\Big{(}(A(t)^{T}Q+QA(t))\otimes Q\Big{)}\vec{X}=\\ \begin{array}[]{l}\cdots=\vec{X}^{T}\mbox{vec}\Big{(}QX(A(t)^{T}Q+QA(t))\Big{)}\\ \cdots=\mbox{trace}\Big{(}X^{T}QX(A(t)^{T}Q+QA(t))\Big{)}\\ \cdots=\mbox{trace}\Big{(}Q^{1/2}X^{T}(A(t)^{T}Q+QA(t))XQ^{1/2}\Big{)}\end{array}$$

$$\frac{d}{dt}\xi =(I\otimes \mathcal{A}(t)+\mathcal{A}(t)\otimes I)\xi ,$$

$$V(\xi )={\xi }^T(Q\otimes Q\otimes Q\otimes Q)\xi$$

$${\mathcal{A}}_{c,i}^T{P}_c+P_c{\mathcal{A}}_{c,i}<0$$

$${\mathcal{A}}_{c,i}:=\sum _{j=0}^{c-1}{I}_{n^j}\otimes A_i\otimes I_{n^{c-1-j}}.$$

$$\dot{\vec{X}}={\mathcal{A}}_c(t)\vec{X}$$

$$V_{c}(x)=\big{(}\otimes^{c}x(t)^{T}\big{)}P_{c}\big{(}\otimes^{c}x(t)\big{)}$$

$$\vec{X}={\left[\begin{array}{cccc}{x}_1^2& x_1{x}_2& x_1{x}_2& x_2^2\end{array}\right]}^T.$$

Full text

I Introduction

Switched dynamical system models appear throughout the field of control theory, and the structure of such models has been widely explored and exploited in order to the analyze stability and performance of real-world systems [1, 2]. In turn, such results have inspired the use of switched systems as a modeling tool for many challenging analysis problems. For example, hybrid dynamical systems can be represented as switched systems, as can some stochastic systems [3, 4]. Certain nonlinearities such as saturation and mechanical backlash can be modeled using switched linear systems [1, 5, 6, 7], as can random noise [8]. Additionally, switched linear systems can be used as an over-approximating abstraction for more general nonlinearities [5, 9] and, for this reason, switched linear system models appear widely in robustness analysis literature [6, 10]. Further, the consistent use of switched system models in safety-critical applications has facilitated the need for computationally efficient analysis tools.

Stability-type proofs for switched dynamical systems often require the construction of polynomial Lyapunov functions. Such proofs guarantee system stability by associating a global energy field with the system state space and then showing that energy is decreasing for all initial conditions and all switched modes. The simplest class of polynomial Lyapunov function is the class of quadratic Lyapunov functions and, as such, the search for quadratic Lyapunov functions has computational advantages in comparison to other methods of stability analysis. Numerous works, including [7, 11] explore the guarantees attainable when solely searching for quadratic Lyapunov functions, however, recent progress in sum-of-squares based techniques have shown that higher-order polynomial Lyapunov functions can be calculated as well with more accurate stability guarantees [12, 13]. In general, sum-of-squares based techniques require little machinery to implement; these methods cast the search for a polynomial Lyapunov function as a convex feasibility problem, and many efficient solvers exist to solve such problems [12]. Additionally, for systems which are known to be stable, the computation of high-order polynomial Lyapunov functions has the ability to help characterize invariant regions of the state space with complex geometries; this is not possible when computing quadratic Lyapunov functions.

This work provides an algorithm for constructing homogeneous polynomial Lyapunov functions for switched linear systems. The aforementioned algorithm searches for polynomial Lyapunov functions through a convex feasibility problem, however, the structure of our algorithm differs significantly from traditional sum-of-squares formulations. Specifically, we encode the search for polynomial Lyapunov functions as a search for quadratic Lyapunov functions for a related hierarchy of Lyapunov differential equations. This creates an intuitive procedure for checking the stability properties of switched linear systems and enables new applications as well [8]. Moreover, we show that every homogeneous sum-of-squares polynomial Lyapunov function for a given initial system can be transformed to a quadratic polynomial Lyapunov function for a system in the related hierarchy; this procedure can also be conducted in the reverse order, allowing one to generate sum-of-squares polynomial Lyapunov functions for an initial system through the identification of a quadratic polynomial Lyapunov function for a related system.

This paper is organized in the following way. We review common analysis tools for assessing the stability of switched linear systems in Section II. Specifically, we introduce a time-varying Lyapunov differential equation, which we define in reference to an initial switched linear system. Using the time-varying Lyapunov differential equation as an initial case, we then form a hierarchy of Lyapunov differential equations in Section III. Quadratic Lyapunov functions for differential equations in this hierarchy are shown to correspond to homogeneous polynomial Lyapunov functions for the initial switched system later in the same section. Section IV explores the relation between quadratic Lyapunov functions for the aforementioned hierarchy of Lyapunov differential equations and homogeneous sum-of-squares polynomial Lyapunov functions for the initial switched linear system. Finally, we provide an algorithm, formulated as a convex optimization problem, for computing high-order homogeneous polynomial Lyapunov functions for switched linear systems in Section V; this algorithm is presented in conjunction with a numerical example.

II Stability and Switched Linear Systems

II-A Preliminaries

Consider the linear time-variant system

[TABLE]

where $x(t)\in\mathbb{R}^{n}$ denotes the system state, and $A(t)\in\mathbb{R}^{n\times n}$ evolves nondeterministically inside a finite set of switched linear modes $A(t)\in\left\{A_{1},\,\cdots,\,A_{N}\right\}$. We assume that each of the switched modes $\dot{x}=A_{i}x$, with $i\in\{1,\cdots,N\}$, converges asymptotically to the origin for all initial conditions $x(0)\in\mathbb{R}^{n}$.

Importantly, the asymptotic stability of each mode does not, by itself, imply the asymptotic stability of the system (1) under arbitrary switching; see [1] Chapter 2 for further details. As such, more complex techniques are required to analyze the stability of (1).

In this work, we consider a traditional approach for stability analysis for switched linear systems, involving the search for a common polynomial Lyapunov function that stabilizes each switched mode (Definition 1).

Definition 1**.**

A common Lyapunov function for the system (1) is a mapping $V:\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that

[TABLE]

It is well known that the system (1) is stable if and only if there exists a Lyapunov function $V(x)$ which satisfies (2). Moreover, the authors of [14] show that (1) is stable if and only if there exists a common homogeneous polynomial Lyapunov function which proves stability of each mode. We capture this assertion in Remark 1.

Remark 1**.**

[15, Theorem 4.5] If the switched linear system (1) is asymptotically stable under arbitrary switching, then there exists a polynomial Lyapunov function $V(x)$, satisfying (2), which is homogeneous in the entries of $x$.

II-B Quadratic Lyapunov Functions for Switched Systems

Reconsider the system (1). In the special instance that there exists a $V(x)$, satisfying (2), which is quadratic in the entries of $x$, we say that the system (1) is quadratically stable [7]. Such a Lyapunov function will take the form

[TABLE]

where $P\in\mathbb{R}^{n\times n}$ is a symmetric positive definite matrix and

[TABLE]

for all $i\,\in\,\{1,\,\cdots,\,N\}$. Alternatively, one can show that the system (1) is quadratically stable by showing that there exists a symmetric positive definite $Q\in\mathbb{R}^{n\times n}$ that satisfies

[TABLE]

for all $i\,\in\,\{1,\,\cdots,\,N\}$ [7]; in this case,

[TABLE]

is a quadratic Lyapunov function for the system (1).

Quadratic polynomial Lyapunov functions are the simplest substantiation of homogeneous polynomial Lyapunov functions, and thus, the search for a quadratic Lyapunov function for (1) has computational advantages in comparison to other strategies for stability analysis; the search can be reduced to solving a convex feasibility problem involving linear matrix inequalities, and many efficient solvers exist to solve such problems [6, 16]. Recent progress in polynomial optimization systems via sum-of-squares relaxations, however, has shown that more general polynomial Lyapunov functions could be computed as well with added benefits, such as improved system stability margins.

Importantly, if the system (1) is linear time-invariant, i.e. $N=1$, then (1) is asymptotically stable if and only if there exists a $P,Q\in\mathbb{R}^{n\times n}$ satisfying (3) and (4), respectively. This is not true, however, in the case of multiple switched modes; stable switched linear systems exist for which there is no quadratic Lyapunov function certifying the stability of each mode [11, Section 3]. For this reason, we must resort to more complex tools to prove stability in the general setting of (1).

II-C The Lyapunov Differential Equation

We next present the time-variant switched Lyapunov differential equation:

[TABLE]

where $X(t)\in\mathbb{R}^{n\times n}$ and $A(t)$ retains its definition from (1). In this work, we primarily use the Lyapunov differential equation (5) as a stability analysis tool for the initial switched system (1). As is shown in the following proposition, (5) is stable if and only if (1) is stable; moreover, stability guarantees on the Lyapunov differential equation propagate down to stability guarantees on the initial system.

Proposition 1**.**

The switched Lyapunov differential equation (5) is stable if and only if the system (1) is stable.

Proof.

\bigcirc$$\Rightarrow Assume the system (5) is stable, and let $X=xx^{T}$. Then

[TABLE]

Therefore $X=xx^{T}$ converges to zero, which implies $x$ converges to zero along trajectories of (1).

\bigcirc$$\Leftarrow Assume the system (1) is stable, and define $X(t)\in\mathbb{R}^{n\times n}$ with initial condition $X(0)=X_{0}$. Any matrix can be written as the sum of diads; therefore, there exist $p_{1,0},\,\cdots,\,p_{N,0},q_{1,0},\,\cdots,\,q_{N,0}\in\mathbb{R}^{n}$ such that

[TABLE]

Next, consider the $2N$ trajectories that satisfy

[TABLE]

[TABLE]

where $j\in\{1,\,\cdots,\,N\}$, and note that if (1) is stable then (6) and (7) converge to zero.

Taking $X=\sum_{j=1}^{N}p_{j}(t)q_{j}(t)^{T}$ then yields

[TABLE]

Therefore, $X=\sum_{j=1}^{N}p_{j}(t)q_{j}(t)^{T}$ is a (unique) solution to the differential equation (5) with initial condition $X_{0}$, and $p_{j}(t)$ and $q_{j}(t)$ are stable for all $j\in\{1,\,\cdots,\,N\}$. Therefore, $X(t)$ also converges along trajectories of (5). ∎

III Establishing a Hierarchy of Lyapunov Differential Equations

In this section we build on (5) to create a hierarchy of Lyapunov differential equations for the system (1). As was the case in Proposition 1, each system in the hierarchy is shown to have equivalent stability properties.

III-A Notation

Let $A\otimes B\in\mathbb{R}^{np\times mq}$ denote the Kronecker product of $A\in\mathbb{R}^{n\times m}$ and $B\in\mathbb{R}^{p\times q}$. Let $\otimes^{k}x\in\mathbb{R}^{n^{k}}$ denote the $k^{\text{th}}$ Kronecker power of $x\in\mathbb{R}^{n}$, which is defined recursively by

[TABLE]

Let $W^{+}\in\mathbb{R}^{m\times n}$ denote the Moore-Penrose inverse of $W\in\mathbb{R}^{n\times m}$, and let $I_{n}\in\mathbb{R}^{n\times n}$ denote the $n\times n$ identity matrix.

III-B Identifying Meta-Lyapunov Functions

We first rewrite (5) as

[TABLE]

by taking $\vec{X}$ to be the vectorization of $X$, i.e. $\vec{X}=\mbox{vec}(X)\in\mathbb{R}^{n^{2}}$. In this case, $\mathcal{A}(t)\in\mathbb{R}^{n^{2}\times n^{2}}$ evolves nondeterministically in the set $\mathcal{A}(t)\in\{\mathcal{A}_{1},\cdots,\mathcal{A}_{N}\}$, where $\mathcal{A}_{i}$ is defined by

[TABLE]

for $i\in\{1,\,\cdots,\,N\}$.

For convenience, we refer to (8), which is also linear time-variant, as the meta-system relative to system (1). Applying concepts of quadratic stability to meta-systems, the system (8) is stable if there exists a positive definite $P\in\mathbb{R}^{n^{2}\times n^{2}}$ such that

[TABLE]

for all $i\in\{1,\,\cdots,\,N\}$. These constraints correspond to the existence of a Lyapunov function $\mathcal{V}(\vec{X})=\vec{X}^{T}P\vec{X}$ for (8), which is quadratic in the entries of $\vec{X}$. In what follows, we refer to $\mathcal{V}(\vec{X})$ as a meta-Lyapunov function for the system (1), and we formalize the search for such a meta-Lyapunov function as the main inquiry of the section.

Problem 1**.**

Given a system (1), which is known to be stable, find a positive definite matrix $P\in\mathbb{R}^{n^{2}\times n^{2}}$ that satisfies (9).

As was shown in Proposition 1, if the system (1) is stable, then there must be a Lyapunov function $\mathcal{V}(\vec{X})$ that certifies the stability of the meta-system (8); this Lyapunov function, however, need not be quadratic. In what follows, we show that in the special instance that (1) is quadratically stable, there must exist a quadratic Lyapunov function certifying the stability of the meta-system (8), and moreover, there must be a $P\in\mathbb{R}^{n^{2}\times n^{2}}$ that solves Problem 1. We capture this assertion in the following theorem.

Theorem 1**.**

If the system (1) is quadratically stable, then the system (8) is also quadratically stable.

Proof.

Assume there exists of a quadratic Lyapunov function $V(x)=x^{T}Qx$ for (1), and pick $P=Q\otimes Q$. Indeed $P$ is positive definite in the instance $Q$ is positive definite. Moreover,

[TABLE]

for all nonzero $\vec{X}\in\mathbb{R}^{n^{2}}$.

We next show that $\mathcal{V}$ decreases along the trajectories of (8). From (8) we have

[TABLE]

Further, we calculate

[TABLE]

and

[TABLE]

Grouping terms then yields

[TABLE]

We now check that $\dot{\mathcal{V}}(\vec{X})$ is negative for all nonzero $\vec{X}\in\mathbb{R}^{n^{2}}$. To that end, note that

[TABLE]

and

[TABLE]

Since $A(t)^{T}Q+QA(t)$ is negative semidefinite, so is $Q^{1/2}X^{T}(A^{T}Q+QA)XQ^{1/2}$, and its trace is negative. Thus $\dot{\mathcal{V}}(\vec{X})$ is negative for all nonzero $\vec{X}\in\mathbb{R}^{n^{2}}$, and moreover, $\mathcal{V}(\vec{X})=\vec{X}^{T}(Q\otimes Q)\vec{X}$ is a quadratic Lyapunov function certifying the stability of the meta system (8). Additionally, this result confirms that $P=Q\otimes Q$ solves problem 1.

∎

It is of course possible to repeat the process again and certify stability at a deeper level; for instance, one may form the Lyapunov differential equation corresponding to (8),

[TABLE]

$\xi\in\mathbb{R}^{n^{4}}$ and then show that

[TABLE]

is a quadratic Lyapunov function for the new meta-system (10). Pursuing the process further, it is possible to construct a “hierarchy” of Lyapunov differential equations whose state space dimensions are $n^{2^{c}}$, where $c$ is an integer greater than or equal to $1$. In the following section, we complete this hierarchy to include Lyapunov differential equations whose state space dimensions grow as $n^{2c}$.

III-C A Linear Hierarchy of Polynomial Lyapunov Functions

We next develop a hierarchy of dynamical systems whose state space dimensions grow as integer exponents of $n$, the dimension of the state space of (1). This hierarchy complements the hierarchy of systems discussed above.

Theorem 2**.**

System (1) is stable if there exists $c\in\mathbb{N}_{\geq 1}$ and $P_{c}\in\mathbb{R}^{n^{c}\times n^{c}}$ positive definite such that

[TABLE]

for all $i\in\{1,\,\cdots,\,N\}$, where

[TABLE]

Proof.

Taking $\vec{X}=\otimes^{c}x(t)\in\mathbb{R}^{n^{c}}$, we find

[TABLE]

where $\mathcal{A}_{c}$ is given by (11), and the stability of system (13) implies that of System (1). Therefore System (1) is stable if there exists a positive definite $P_{c}\in\mathbb{R}^{n^{c}\times n^{c}}$ such that (11) holds. ∎

Theorem 2 shows that the existence of a $P_{c}\in\mathbb{R}^{n^{c}\times n^{c}}$ satisfying (11) for some integer $c\geq 1$ certifies the stability of (1); such a $P_{c}$ identifies

[TABLE]

as a polynomial Lyapunov function for (1), which is homogeneous in the entries of $x$ and of order $2c$. Importantly, the degree of $V_{c}(x)$ grows linearly with $c$.

III-D Reducing the Dimensionality of the Meta-System

The benefits of searching for meta-Lyapunov functions for (1) using the methods presented thus far are namely structural; (11)-(12) provide an intuitive procedure for generating high-order homogeneous polynomial Lyapunov function for (1) and moreover, this procedure does not require any heavy machinery to implement. In contrast, there are few computational advantages to this approach, at present. This is due in part to internal redundancy built into the Lyapunov constraints given by (11). We demonstrate this assertion through the following example.

Example 1**.**

Consider, for example, the system (1) evolving in $\mathbb{R}^{2}$. In this case, $x=[x_{1},\,x_{2}]^{T}\in\mathbb{R}^{2}$, and $\vec{X}:=x\otimes x\in\mathbb{R}^{4}$ is given by

[TABLE]

When beginning at an initial condition $\vec{X}(0)=x(0)\otimes x(0)$ and evolving along trajectories of the meta-system

[TABLE]

we find that the second and third entries of $\vec{X}$ remain equal to one another, regardless of the switching policy. This is due to the construction of $(I_{2}\otimes A(t)+A(t)\otimes I_{2})$.

The methods presented thus far address the problem of searching for a meta-Lyapunov function $\mathcal{V}(\vec{X})=\vec{X}^{T}P\vec{X}$ for the system (15); the specific choice of $P\in\mathbb{R}^{4\times 4}$ will then correspond to a homogeneous polynomial Lyapunov function $V_{2}(x)=(\otimes^{2}x)^{T}P(\otimes^{2}x)$ for the system (1). Here, it is apparent that the constraints on $P$, given by (11), contain internal redundancy; note, for instance, that one must compute the $10$ unique entries of $P\in\mathbb{R}^{4\times 4}$ in order to find $\mathcal{V}(\vec{X})$, whereas, the resulting Lyapunov function $V_{2}(x)$ will only be defined by $5$ unique monomials.

Now, consider a vector containing the second-order monomials of $x$, this time with no redundancy. Specifically, consider $y(x)=[x_{1}^{2},\,x_{1}x_{2},\,x_{2}^{2}]\in\mathbb{R}^{3}$, and note that $\vec{X}=Wy(x)$ where

[TABLE]

Using (13) as a basis, the dynamics of $y(x)$ can be captured in closed form:

[TABLE]

Therefore, one can now formulate the search for a fourth-order homogeneous polynomial Lyapunov function for (1), as the search for a quadratic Lyapunov function $\overline{V}(y)=y^{T}\overline{P}y$ that certifies the stability of $y$. In this case, the resulting Lyapunov function will have the same number of terms, i.e. 5 distinct monomials, however this search will only require the identification of the $6$ unique entries of $\overline{P}\in\mathbb{R}^{3\times 3}$. ∎

As shown in the previous example, the constraints given by (11) are redundant; that is, a quadratic Lyapunov function that certifies the stability of $\vec{X}$, as in (13), will individually certify the stability of each of the meta-system’s states, whereas, a reduced order meta-Lyapunov function that stabilizes a subset of meta-system’s states may be sufficient.

For this reason, we present a new formulation of the constraints (11)-(12) that contains no redundancy. We begin with the following definition.

Definition 2** ($\mathcal{A}_{c}$-Invariant Subspaces).**

A subspace $S\subset\mathbb{R}^{n}$ is said to be $\mathcal{A}_{c}$-invariant for (13) if for every vector $v\in S$ and every matrix $\mathcal{A}_{c,i}$ with $i\in\{1,\,\cdots,\,N\}$ we have $\mathcal{A}_{c,\,i}v\in S$.

Note that $\mathcal{A}_{c}(t)$ as in (13) will have an inherent invariant subspace, resulting from its construction. We therefore remove this redundancy by analysing a reduced order meta-system, whose states correspond to unique monomials of the initial switched system (1). While the initial meta-Lyapunov conditions (11) are defined by $n^{2c}$ constraints per switched mode, our new formulation only requires $M(n,\,c)^{2}$ such constraints, where $M(n,\,c)$ denotes the number of monomials of order $c\in\mathbb{N}_{\geq 1}$ in the entries of $x\in\mathbb{R}^{n}$ and is given by

[TABLE]

This result is encapsulated in the following theorem.

Theorem 3**.**

Let $y_{c}(x)\in\mathbb{R}^{M(n,\,c)}$ denote a vector containing the monomials of $x$ of order $c$, which we define in conjunction with a matrix $W_{c}\in\mathbb{R}^{n\times M(n,\,c)}$:

[TABLE]

Additionally, let $\overline{P}_{c}\in R^{M(n,\,c)\times M(n,\,c)}$ be symmetric positive definite. If

[TABLE]

for all $i\in\{1,\,\cdots,\ N\}$ where

[TABLE]

then $V_{c}(x)=y_{c}(x)^{T}\,\overline{P}\,y_{c}(x)$ is a homogeneous polynomial Lyapunov function for (1) of order $2c$.

The proof of this result comes from the fact that $A_{c}(t)$ has an inherent invariant subspace, resulting from its construction. As trajectories of

[TABLE]

are known to begin in this subspace, we can encode the search for meta-Lyapunov functions for the system (1) as a search for quadratic Lyapunov functions for the reduced order system

[TABLE]

where $B_{c}(t)\in\{B_{c,1},\,\cdots,\,B_{c,N}\}$ and $y(x)$ is given by (16).

Importantly, Theorem 3 allows the system designer to select $y_{c}(x)$ with whatever ordering properties they like; that is, we do not assume an order to the monomials that are stored in $y_{c}(x)$. However, each ordering will induce a unique $W_{c}$, and thus the resulting Lyapunov conditions will always be the same, regardless of the chosen ordering. Moreover, the constraints given by (17) are equivalent to the constraints given by (11), now with reduced dimensionality.

In the specific case where $n=2$, there is an intuitive ordering to the monomials of $x$; under this assumed ordering, the matrix $W_{c}$, given by (16), can be captured in closed form.

Proposition 2**.**

Consider the system (1) and let $n=2$. In this case we have

[TABLE]

Additionally, for a positive integer $k\in\mathbb{N}_{\geq 0}$, let $0_{k}\in\mathbb{R}^{k}$ denote a vector populated with zeros.

If $y_{c}(x)\in\mathbb{R}^{c+1}$ conforms to the ordering

[TABLE]

then we have $\otimes^{c}x=W_{c}\,y_{c}(x)$, where for an integer $k\in\mathbb{N}_{\geq 1}$ we define $W_{k}$ recursively by

[TABLE]

In the case when $n>2$, it is generally difficult to order the $c^{\text{th}}$ order monomials of $x$ in an intuitive way. For this reason, we do not expand Proposition 2 to account for the case where $n>2$, nor do we suggest a canonical ordering for the entries of $y_{c}(x)$. However, $W_{c}$ can always be solved for using (16) once $y_{c}(x)$ has been chosen.

IV Relation to Homogeneous Polynomial Lyapunov Functions

Traditionally, the search for a polynomial Lyapunov functions systems of the form (1) is encoded as the search for a sum-of-squares polynomial $V(x)$, satisfying (2).

Definition 3**.**

A polynomial $p(x)$ is a sum-of-squares in $x$ if there exist polynomials $g_{1},\cdots,g_{r}$ such that

[TABLE]

The search for a sum-of-squares polynomial $V(x)$, satisfying (2), is known to be a convex optimization problem, computable by solving a semidefinite program [13]. Many efficient solvers exist to handle such problems [6, 16].

We next show that the existence of quadratic Lyapunov functions for the hierarchy of dynamical systems (13) guarantees the existence of a homogeneous sum-of-squares polynomial Lyapunov functions for (1), and vice versa. In this sense, all homogeneous sum-of-squares polynomial Lyapunov functions can be thought of as quadratic Lyapunov functions for a related hierarchy of differential equations. Moreover, one can encode the search for high-order sum-of-squares polynomial Lyapunov functions, which certify the stability of (1), as a search for quadratic Lyapunov functions for the related system (13). Calculating sum-of-squares polynomial Lyapunov functions in this way can be used to reduce the amount of machinery required to certify the stability of general switched linear systems of the form (1).

Theorem 4**.**

There exists a $\overline{P}_{c}\in\mathbb{R}^{M(n,\,c)\times M(n,\,c)}$ satisfying (17) for some positive integer $c\in\mathbb{N}_{\geq 1}$, if and only if there exists a homogeneous sum-of-squares polynomial Lyapunov function $V_{c}(x)$ of degree $2c$ for the system (1).

Proof.

A sum-of-squares polynomial that is homogeneous in the entries $x$ and of order $2c$ will take the form $p(x)=y_{c}(x)^{T}Zy_{c}(x)$, where $Z\in\mathbb{R}^{M(n,c)\times M(n,c)}$ is symmetric, and $y_{c}(x)$ and $M(n,c)$ retain their definitions from Theorem 3. From Theorem 3, we have that if $\overline{P}_{c}\in\mathbb{R}^{M(n,\,c)\times M(n,\,c)}$ satisfies (11) for some positive integer $c\in\mathbb{N}_{\geq 1}$, then we have that $V_{c}(x)=y_{c}(x)^{T}\overline{P}y_{c}(x)$ is a homogeneous polynomial Lyapunov function for (21) and, moreover, $V_{c}(x)$ is a sum-of-squares. To prove the converse, we note that if $p(x)=y_{c}(x)^{T}Zy_{c}(x)$ is a homogeneous sum-of-squares polynomial Lyapunov function for (21) then $Z>0$ and $\dot{p}(x)<0$ for all $x\in\mathbb{R}^{n}$. From the dynamics of $y_{c}(x)$, given as (19), we have $B^{T}_{c,\,i}Z+ZB_{c,i}<0$ for all $i\in\{1,\cdots,N\}$. Therefore $\overline{P}_{c}=Z$ solves (17). ∎

V Numerical Example

In this section, we provide an example case and prove the stability of a switched linear system using a meta-Lyapunov function based approach. An algorithm is provided for generating homogeneous polynomial Lyapunov functions for switched systems, which follows the procedure detailed in Theorem 3; this algorithm is specifically written for implementation with CVX, a convex optimization toolbox made for use with MATLAB [17]. We also provide a comparison to a similar search for homogeneous polynomial Lyapunov functions that was implemented using SOSTOOLS, a sum-of-squares optimization toolbox made for use with MATLAB [18]. Experimental results are provided from MATLAB 2019b, which was run on a 2017 Macbook Pro laptop.

V-A Problem Formulation

We consider the linear time-variant system

[TABLE]

[TABLE]

In the following, we go about showing that (21) is stable. This is done, at first, through the computation of a quadratic Lyapunov function $V_{1}(x)=x^{T}Px$, which satisfies (2), and then through the computation of higher-order homogeneous polynomial Lyapunov functions using the procedure detailed in Theorem 3.

Importantly, if the system (21) begins at an initial position $x_{0}=x(0)$, and there exists a Lyapunov function $V_{c}(x)$ that certifies the stability of (21), then the infinite-time system trajectory is constrained to stay inside

[TABLE]

for all $t\geq 0$. For this reason, we select $V_{c}(x)$ as the minimizers of a suitable objective function, as to shrink the resulting invariant region derived through (22). In what follows, we additionally show that computing higher-order meta-Lyapunov functions allows one to characterise tighter invariant sets by (22), even when the same objective function is used in each computation.

V-B Identifying Meta-Lyapunov Functions

We search for meta-Lyapunov functions for (21) using a semidefinite program. Specifically, when searching for a homogeneous Lyapunov function of order $2c$, we first calculate $B_{c,1}$ and $B_{c,2}$ using equations (12), (18) and (20), and then we search for a symmetric positive-definite matrix $\overline{P}_{c}\in\mathbb{R}^{(c+1)\times(c+1)}$ that satisfies (17). Such a matrix identifies $V_{c}(x)=y_{c}(x)^{T}\,\overline{P}_{c}\,y_{c}(x)$ as a polynomial Lyapunov function for (21), which is homogeneous in the entries of $x$ and of order $2c$. We implement the aforementioned procedure with Algorithm 1, which specifically relies on CVX, a convex optimization toolbox built for use with MATLAB [17, 19]. Algorithm 1 takes as inputs the system parameters $A_{1}$ and $A_{2}$, and a positive integer $c$, and returns a matrix $\overline{P}_{c}$, in the case that one exists, which satisfies (11) at the $c^{\text{th}}$ level.

Note that Algorithm 1 computes the solution to a feasibility problem, rather than an optimization problem; that is, while Algorithm 1 searches for a $\overline{P}_{c}$ that satisfies the meta-Lyapunov constraint (17), this solution is computed without referencing any objective function. Note however, that in the instance that multiple feasible solutions exist, it is preferable to choose $\overline{P}_{c}$ such that the sublevel sets of the resulting homogeneous Lyapunov function $V_{c}(x)=y_{c}(x)^{T}\,\overline{P}_{c}\,y_{c}(x)$ are small; this is due to the fact that $V_{c}(x)$ can be used to find infinite time reachable sets of (21) under arbitrary switching. For this reason, it is desirable to compute $\overline{P}_{c}$ as the solution to an optimisation problem, rather than a feasibility problem.

Little is known, in general about how one can relate the parameters of a polynomial to the volume of its sublevel sets. In our case as well, it is difficult to associate a metric of optimality with the a feasible solution to the meta-Lyapunov constraints (17). Through experimentation, we have generally found that it is preferable to compute numerous solutions using different objective functions, and then compute an invariant region as the intersection of their respective sublevel sets. Specifically we recommend using either using the objective function

which minimises the coefficient on $x_{1}^{2c}$ in the resultant Lyapunov function $V_{c}(x)$, or

which minimises the coefficient on $x_{2}^{2c}$. These objective functions are provided in psuedocode, such that they can easily be inserted in Algorithm 1 at line 11.

V-C Numerical Results and Comparison with SOSTOOLS

We now return to the example system (21), and compute feasible meta-Lyapunov functions with Algorithm 1. Additionally, we compute an over approximation of the infinite time reachable set of (21) when beginning from the initial conditions $x(0)=[1,\,0]^{T}$.

As discussed in the preceding, we compute these invariant sets by implementing Algorithm 1, while attempting to minimize $\overline{P}(1,1)$, i.e. the coefficient on $x_{1}^{2c}$; see Algorithm 1, Line 11. This procedure was computed in MATLAB 2019b using CVX.

In the case of this example, Algorithm 1 was computed for $c\in\{1,2,\cdots,13\}$, thus generating homogeneous polynomial Lyapunov functions for all even orders between 2 and 26. These Lyapunov functions were then used to calculate invariant regions of the state space using (22); see Figure 1. Note that as the order of the meta-Lyapunov function increases, the derived invariant sets shrink in volume. Further, certain higher-order the meta-Lyapunov functions were shown to have non-convex sublevel sets. We provide the number of solver iterations for each experiment, as well as the computations times, in Figure 2.

We next compare the meta-Lyapunov function based method to the more traditional sum-of-squares based approach for calculating invariant regions. Specifically, we search for high-order homogeneous polynomial Lyapunov functions for the system (21) using SOSTOOLS, MATLAB’s sum-of-squares toolbox [18]. We similarly implement SOSTOOLS with the solver SDPT3 and attempt to generate homogeneous polynomial Lyapunov functions for the system (21) while minimizing the coefficient on $x_{1}^{2c}$. As was the case previously, we provide the number of solver iterations for each experiment, as well as the computations times (See Figure 2).

In the experiment, SOSTOOLS was only able to generate homogeneous polynomial Lyapunov functions of order 10 or below; a solver error was returned during each search for more complex Lyapunov function. In contrast, Algorithm 1, when implemented through CVX, was able to generate up to $26^{\text{th}}$-order polynomial Lyapunov functions while minimizing the same objective function. We attribute this discrepancy to the fact that many of the steps required in a traditional sum-of-squares based search optimization are not required by Algorithm 1. For example, Algorithm 1 does not compute the time rate of change of the monomials in $x$ of order $2c$; that is, Algorithm 1 begins with a closed from representation of $\dot{y}_{c}(x)$, which is encoded in the matrices $B_{c,1}$ and $B_{c,2}$. SOSTOOLS must compute $\dot{y}_{c}(x)$ online as a sum of squares of lower-order monomials and, for this reason, one can expect SOSTOOLS to preform less efficiently.

Despite this, in the experiments where SOSTOOLS was able to correctly generate homogeneous Lyapunov functions, we found that SOSTOOLS preformed faster in computation that the meta-Lyapunov search implemented with CVX; however, it did take SDPT3 more solver iterations to generate the solution which minimized the objective function when implemented through SOSTOOLS (Figure 2).

VI Conclusion

This work addresses the problem of searching for homogeneous polynomial Lyapunov functions for stable switched linear systems. An equivalence is shown between polynomial Lyapunov functions for switched linear systems and quadratic Lyapunov functions for a related hierarchy of Lyapunov differential equations. A computationally competitive algorithm is presented for generating high-order homogeneous polynomial Lyapunov functions in this manner.