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

L yapunov Differ ential Equation Hierar chy and Polynomial L yapunov Functions for Switched Linear Systems Matthe w Abate, Corbin Klett, Samuel Coogan, and Eric Feron Abstract — This work studies the problem of sear ching for homogeneous polynomial L yapuno v functions for stable switched linear systems. Specifically , we show an equiv alence between polynomial L yapunov functions f or systems of this class and quadratic L yapunov functions for a r elated hierarch y of L yapunov differential equations. This creates an intuitive procedur e f or checking the stability properties of switched linear systems, and a computationally competitive algorithm is presented f or generating high-order homogeneous polynomial L yapunov functions in this manner . Additionally , we pro vide a comparison between polynomial L yapunov functions generated with our proposed approach and polynomial L yapunov func- tions generated with a more traditional sum-of-squares based approach. I . I N T RO D U C T I O N Switched dynamical system models appear throughout the field of control theory , and the structure of such models has been widely e xplored 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 e xample, 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]–[7], as can random noise [8]. Additionally , switched lin- ear systems can be used as an ov er -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 L yapunov 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 This material is based upon work supported by the United States Govern- ment under Air Force Of fice of Scientific Research grant number F A9550- 19-1-0015. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the Government. M. Abate is with the School of Mechanical Engineering and the School of Electrical and Computer Engineering, Georgia Institute of T echnology , Atlanta, 30332, USA: Matt.Abate@GaTech.edu . C. Klett and E. Feron are with the School of Aerospace En- gineering, Georgia Institute of T echnology , Atlanta, 30332, USA: Corbin@GaTech.edu and Feron@GaTech.edu . S. Coogan is with the School of Electrical and Computer Engineering and the School of Civil and Environmental Engineering, Georgia Institute of T echnology , Atlanta, 30332, USA: Sam.Coogan@GaTech.edu . and all switched modes. The simplest class of polynomial L yapunov function is the class of quadratic L yapunov func- tions and, as such, the search for quadratic L yapunov func- tions has computational adv antages in comparison to other methods of stability analysis. Numerous works, including [7], [11] explore the guarantees attainable when solely searching for quadratic L yapunov functions, ho wever , recent progress in sum-of-squares based techniques hav e shown that higher- order polynomial L yapunov 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 L yapunov function as a conv ex 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 L ya- punov functions has the ability to help characterize in variant regions of the state space with complex geometries; this is not possible when computing quadratic L yapunov functions. This work provides an algorithm for constructing ho- mogeneous polynomial L yapunov functions for switched linear systems. The aforementioned algorithm searches for polynomial L yapunov functions through a con ve x feasibility problem, ho we ver , the structure of our algorithm dif fers significantly from traditional sum-of-squares formulations. Specifically , we encode the search for polynomial L yapunov functions as a search for quadratic L yapunov functions for a related hierarchy of L yapunov differential equations. This creates an intuiti ve procedure for checking the stabil- ity properties of switched linear systems and enables ne w applications as well [8]. Moreover , we show that ev ery homogeneous sum-of-squares polynomial L yapunov function for a gi ven initial system can be transformed to a quadratic polynomial L yapunov function for a system in the related hierarchy; this procedure can also be conducted in the re verse order , allowing one to generate sum-of-squares polynomial L yapunov functions for an initial system through the identi- fication of a quadratic polynomial L yapunov function for a related system. This paper is or ganized in the follo wing way . W e re view common analysis tools for assessing the stability of switched linear systems in Section II. Specifically , we introduce a time-varying L yapunov differential equation, which we de- fine in reference to an initial switched linear system. Using the time-varying L yapunov differential equation as an initial case, we then form a hierarchy of L yapunov dif ferential equations in Section III. Quadratic L yapunov functions for differential equations in this hierarchy are sho wn to cor - respond to homogeneous polynomial L yapunov functions for the initial switched system later in the same section. Section IV explores the relation between quadratic L yapunov functions for the aforementioned hierarchy of L yapunov dif- ferential equations and homogeneous sum-of-squares poly- nomial L yapunov functions for the initial switched linear system. Finally , we provide an algorithm, formulated as a con vex optimization problem, for computing high-order homogeneous polynomial L yapunov functions for switched linear systems in Section V ; this algorithm is presented in conjunction with a numerical e xample. I I . S T A B I L I T Y A N D S W I T C H E D L I N E A R S Y S T E M S A. Pr eliminaries Consider the linear time-variant system ˙ x = A ( t ) x, (1) where x ( t ) ∈ R n denotes the system state, and A ( t ) ∈ R n × n ev olves nondeterministically inside a finite set of switched linear modes A ( t ) ∈ { A 1 , · · · , A N } . W e assume that each of the switched modes ˙ x = A i x , with i ∈ { 1 , · · · , N } , con- ver ges asymptotically to the origin for all initial conditions x (0) ∈ 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, in volving the search for a common polynomial L yapunov function that stabilizes each switched mode (Definition 1). Definition 1. A common L yapunov function for the system (1) is a mapping V : R n → R such that V ( x ) > 0 ˙ V ( x ) = h∇ V , A i x i < 0 ∀ x 6 = 0 ∀ i ∈ { 1 , · · · , N } . (2) It is well kno wn that the system (1) is stable if and only if there e xists a L yapunov function V ( x ) which satisfies (2). Moreo ver , the authors of [14] sho w that (1) is stable if and only if there exists a common homogeneous polynomial L yapunov function which prov es stability of each mode. W e 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 L yapunov function V ( x ) , satisfying (2), which is homogeneous in the entries of x . 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 L yapunov function will take the form V ( x ) = x T P x where P ∈ R n × n is a symmetric positiv e definite matrix and A T i P + P A i < 0 (3) for all i ∈ { 1 , · · · , N } . Alternatively , one can show that the system (1) is quadratically stable by sho wing that there exists a symmetric positi ve definite Q ∈ R n × n that satisfies A i Q + QA T i < 0 (4) for all i ∈ { 1 , · · · , N } [7]; in this case, V ( x ) = x T Q − 1 x is a quadratic L yapunov function for the system (1). Quadratic polynomial L yapunov functions are the simplest substantiation of homogeneous polynomial L yapunov func- tions, and thus, the search for a quadratic L yapunov function for (1) has computational advantages in comparison to other strategies for stability analysis; the search can be reduced to solving a con ve x feasibility problem in volving linear matrix inequalities, and many efficient solvers exist to solve such problems [6], [16]. Recent progress in polynomial optimiza- tion systems via sum-of-squares relaxations, howe ver , has shown that more general polynomial L yapunov functions could be computed as well with added benefits, such as improv ed system stability margins. Importantly , if the system (1) is linear time-in variant, i.e. N = 1 , then (1) is asymptotically stable if and only if there exists a P , Q ∈ R n × n satisfying (3) and (4), respecti vely . This is not true, howe ver , in the case of multiple switched modes; stable switched linear systems exist for which there is no quadratic L yapunov function certifying the stability of each mode [11, Section 3]. For this reason, we must resort to more comple x tools to pro ve stability in the general setting of (1). C. The Lyapuno v Differ ential Equation W e next present the time-v ariant switched L yapunov dif- ferential equation: ˙ X = A ( t ) X + X A ( t ) T , (5) where X ( t ) ∈ R n × n and A ( t ) retains its definition from (1). In this work, we primarily use the L yapunov differential equation (5) as a stability analysis tool for the initial switched system (1). As is sho wn in the follo wing proposition, (5) is stable if and only if (1) is stable; moreov er, stability guarantees on the L yapunov dif ferential equation propagate down to stability guarantees on the initial system. Proposition 1. The switched L yapunov differ ential equation (5) is stable if and only if the system (1) is stable . Pr oof.  ⇒ Assume the system (5) is stable, and let X = xx T . Then ˙ X = ˙ xx T + x ˙ x T = A ( t ) xx T + xx T A ( t ) T = A ( t ) X + X A ( t ) T . Therefore X = xx T con verges to zero, which implies x con verges to zero along trajectories of (1).  ⇐ Assume the system (1) is stable, and define X ( t ) ∈ R n × n with initial condition X (0) = X 0 . Any matrix can be written as the sum of diads; therefore, there exist p 1 , 0 , · · · , p N , 0 , q 1 , 0 , · · · , q N , 0 ∈ R n such that X 0 = N X j =1 p j, 0 q T j, 0 . Next, consider the 2 N trajectories that satisfy d dt p j ( t ) = A ( t ) p j ( t ) , p j (0) = p j, 0 , (6) d dt q j ( t ) = A ( t ) q j ( t ) , q j (0) = q j, 0 . (7) where j ∈ { 1 , · · · , N } , and note that if (1) is stable then (6) and (7) conv erge to zero. T aking X = P N j =1 p j ( t ) q j ( t ) T then yields ˙ X ( t ) = N X j =1  ˙ p j q T j + p j ˙ q T j  = A ( t )  N X j =1 p j q T j  +  N X j =1 p j q T j  A ( t ) T = A ( t ) X ( t ) + X ( t ) A ( t ) T . Therefore, X = P N j =1 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 ∈ { 1 , · · · , N } . Therefore, X ( t ) also con ver ges along trajectories of (5). I I I . E S T A B L I S H I N G A H I E R A R C H Y O F L Y A P U N O V D I FF E R E N T I A L E Q U A T I O N S In this section we b uild on (5) to create a hierarchy of L yapunov dif ferential equations for the system (1). As was the case in Proposition 1, each system in the hierarchy is shown to ha ve equiv alent stability properties. A. Notation Let A ⊗ B ∈ R np × mq denote the Kronecker product of A ∈ R n × m and B ∈ R p × q . Let ⊗ k x ∈ R n k denote the k th Kronecker power of x ∈ R n , which is defined recursiv ely by ⊗ 1 x = x ∈ R n , ⊗ k x = x ⊗ ( ⊗ k − 1 x ) ∈ R n k , k ≥ 2 . Let W + ∈ R m × n denote the Moore-Penrose inv erse of W ∈ R n × m , and let I n ∈ R n × n denote the n × n identity matrix. B. Identifying Meta-Lyapuno v Functions W e first re write (5) as ˙ ~ X = A ( t ) ~ X (8) by taking ~ X to be the vectorization of X , i.e. ~ X = vec ( X ) ∈ R n 2 . In this case, A ( t ) ∈ R n 2 × n 2 ev olves nondeterministi- cally in the set A ( t ) ∈ {A 1 , · · · , A N } , where A i is defined by A i := I n ⊗ A i + A i ⊗ I n for i ∈ { 1 , · · · , N } . For con venience, we refer to (8), which is also linear time- variant, as the meta-system relati ve to system (1). Applying concepts of quadratic stability to meta-systems, the system (8) is stable if there e xists a positi ve definite P ∈ R n 2 × n 2 such that A T i P + P A i < 0 , (9) for all i ∈ { 1 , · · · , N } . These constraints correspond to the existence of a L yapunov function V ( ~ X ) = ~ X T P ~ X for (8), which is quadratic in the entries of ~ X . In what follo ws, we refer to V ( ~ X ) as a meta-L yapunov function for the system (1), and we formalize the search for such a meta-L yapunov function as the main inquiry of the section. Problem 1. Gi ven a system (1), which is known to be stable, find a positi ve definite matrix P ∈ R n 2 × n 2 that satisfies (9). As was shown in Proposition 1, if the system (1) is stable, then there must be a L yapunov function V ( ~ X ) that certifies the stability of the meta-system (8); this L yapunov function, howe ver , need not be quadratic. In what follo ws, we sho w that in the special instance that (1) is quadratically stable, there must exist a quadratic L yapunov function certifying the stability of the meta-system (8), and moreov er , there must be a P ∈ R n 2 × n 2 that solves Problem 1. W e capture this assertion in the following theorem. Theorem 1. If the system (1) is quadratically stable, then the system (8) is also quadratically stable. Pr oof. Assume there exists of a quadratic L yapunov function V ( x ) = x T Qx for (1), and pick P = Q ⊗ Q . Indeed P is positiv e definite in the instance Q is positi ve definite. Moreov er , V ( ~ X ) := ~ X T P ~ X > 0 for all nonzero ~ X ∈ R n 2 . W e next show that V decreases along the trajectories of (8). From (8) we ha ve ˙ V ( ~ X ) = ~ X T  A ( t ) T P + P A ( t )  ~ X . Further , we calculate A ( t ) T P = ( I n ⊗ A ( t ) T + A ( t ) T ⊗ I n ) Q ⊗ Q = Q ⊗ ( A ( t ) T Q ) + ( A ( t ) T Q ) ⊗ Q and P A ( t ) = Q ⊗ Q ( I n ⊗ A ( t ) + A ( t ) ⊗ I n ) = Q ⊗ ( QA ( t )) + ( QA ( t )) ⊗ Q. Grouping terms then yields A ( t ) T P + P A ( t ) = Q ⊗ ( A ( t ) T Q + QA ( t )) + · · · + ( A ( t ) T Q + QA ( t )) ⊗ Q. W e now check that ˙ V ( ~ X ) is negati ve for all nonzero ~ X ∈ R n 2 . T o that end, note that ~ X T  Q ⊗ ( A ( t ) T Q + QA ( t ))  ~ X = · · · = ~ X T vec  ( A ( t ) T Q + QA ( t )) X Q  · · · = trace  X T ( A ( t ) T Q + QA ( t )) X Q  · · · = trace  Q 1 / 2 X T ( A ( t ) T Q + QA ( t )) X Q 1 / 2  , and ~ X T  ( A ( t ) T Q + QA ( t )) ⊗ Q  ~ X = · · · = ~ X T vec  QX ( A ( t ) T Q + QA ( t ))  · · · = trace  X T QX ( A ( t ) T Q + QA ( t ))  · · · = trace  Q 1 / 2 X T ( A ( t ) T Q + QA ( t )) X Q 1 / 2  Since A ( t ) T Q + QA ( t ) is negati ve semidefinite, so is Q 1 / 2 X T ( A T Q + QA ) X Q 1 / 2 , and its trace is negati ve. Thus ˙ V ( ~ X ) is negati ve for all nonzero ~ X ∈ R n 2 , and moreover , V ( ~ X ) = ~ X T ( Q ⊗ Q ) ~ X is a quadratic L yapunov function certifying the stability of the meta system (8). Additionally , this result confirms that P = Q ⊗ Q solv es 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 L yapunov differential equation corresponding to (8), d dt ξ = ( I ⊗ A ( t ) + A ( t ) ⊗ I ) ξ , (10) ξ ∈ R n 4 and then sho w that V ( ξ ) = ξ T ( Q ⊗ Q ⊗ Q ⊗ Q ) ξ is a quadratic L yapunov function for the new meta-system (10). Pursuing the process further, it is possible to construct a “hierarchy” of L yapunov differential equations whose state space dimensions are n 2 c , where c is an integer greater than or equal to 1 . In the follo wing section, we complete this hierarchy to include L yapunov dif ferential equations whose state space dimensions grow as n 2 c . C. A Linear Hierar chy of P olynomial Lyapuno v Functions W e next de velop 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 hierarch y of systems discussed abo ve. Theorem 2. System (1) is stable if ther e e xists c ∈ N ≥ 1 and P c ∈ R n c × n c positive definite suc h that A T c,i P c + P c A c,i < 0 (11) for all i ∈ { 1 , · · · , N } , wher e A c,i := c − 1 X j =0 I n j ⊗ A i ⊗ I n c − 1 − j . (12) Pr oof. T aking ~ X = ⊗ c x ( t ) ∈ R n c , we find ˙ ~ X = A c ( t ) ~ X (13) where A c is gi ven by (11), and the stability of system (13) implies that of System (1). Therefore System (1) is stable if there exists a positiv e definite P c ∈ R n c × n c such that (11) holds. Theorem 2 sho ws that the existence of a P c ∈ R n c × n c satisfying (11) for some integer c ≥ 1 certifies the stability of (1); such a P c identifies V c ( x ) = ( ⊗ c x ( t ) T ) P c ( ⊗ c x ( t )) (14) as a polynomial L yapunov function for (1), which is homo- geneous in the entries of x and of order 2 c . Importantly , the degree of V c ( x ) grows linearly with c . D. Reducing the Dimensionality of the Meta-System The benefits of searching for meta-L yapunov functions for (1) using the methods presented thus far are namely struc- tural; (11)-(12) provide an intuitiv e procedure for generating high-order homogeneous polynomial L yapunov function for (1) and moreover , this procedure does not require any heavy machinery to implement. In contrast, there are fe w computational adv antages to this approach, at present. This is due in part to internal redundancy built into the L yapunov constraints giv en by (11). W e demonstrate this assertion through the follo wing e xample. Example 1. Consider , for example, the system (1) ev olving in R 2 . In this case, x = [ x 1 , x 2 ] T ∈ R 2 , and ~ X := x ⊗ x ∈ R 4 is giv en by ~ X =  x 2 1 x 1 x 2 x 1 x 2 x 2 2  T . (15) When beginning at an initial condition ~ X (0) = x (0) ⊗ x (0) and ev olving along trajectories of the meta-system ˙ ~ X = ( I 2 ⊗ A ( t ) + A ( t ) ⊗ I 2 ) ~ X , we find that the second and third entries of ~ X remain equal to one another , regardless of the switching policy . This is due to the construction of ( I 2 ⊗ A ( t ) + A ( t ) ⊗ I 2 ) . The methods presented thus far address the problem of searching for a meta-L yapunov function V ( ~ X ) = ~ X T P ~ X for the system (15); the specific choice of P ∈ R 4 × 4 will then correspond to a homogeneous polynomial L yapunov function V 2 ( x ) = ( ⊗ 2 x ) T P ( ⊗ 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 ∈ R 4 × 4 in order to find V ( ~ X ) , whereas, the resulting L yapunov function V 2 ( x ) will only be defined by 5 unique monomials. Now , consider a vector containing the second-order mono- mials of x , this time with no redundancy . Specifically , consider y ( x ) = [ x 2 1 , x 1 x 2 , x 2 2 ] ∈ R 3 , and note that ~ X = W y ( x ) where W =     1 0 0 0 1 0 0 1 0 0 0 1     . Using (13) as a basis, the dynamics of y ( x ) can be captured in closed form: ˙ y = W + A c ( t ) W y . Therefore, one can no w formulate the search for a fourth- order homogeneous polynomial L yapunov function for (1), as the search for a quadratic L yapunov function V ( y ) = y T P y that certifies the stability of y . In this case, the resulting L yapunov function will hav e the same number of terms, i.e. 5 distinct monomials, howe ver this search will only require the identification of the 6 unique entries of P ∈ R 3 × 3 . As shown in the pre vious example, the constraints giv en by (11) are redundant; that is, a quadratic L yapunov function that certifies the stability of ~ X , as in (13), will indi vidually certify the stability of each of the meta-system’ s states, whereas, a reduced order meta-L yapunov function that sta- bilizes a subset of meta-system’ s states may be sufficient. For this reason, we present a ne w formulation of the constraints (11)-(12) that contains no redundanc y . W e begin with the follo wing definition. Definition 2 ( A c -In variant Subspaces) . A subspace S ⊂ R n is said to be A c -in variant for (13) if for every vector v ∈ S and e very matrix A c,i with i ∈ { 1 , · · · , N } we hav e A c, i v ∈ S . Note that A c ( t ) as in (13) will hav e an inherent in vari- ant subspace, resulting from its construction. W e therefore remov e 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- L yapunov conditions (11) are defined by n 2 c 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 ∈ N ≥ 1 in the entries of x ∈ R n and is giv en by M ( n, c ) =  c + n − 1 n − 1  . This result is encapsulated in the follo wing theorem. Theorem 3. Let y c ( x ) ∈ R M ( n, c ) denote a vector containing the monomials of x of or der c , whic h we define in conjunction with a matrix W c ∈ R n × M ( n, c ) : ⊗ c x = W c y c ( x ) . (16) Additionally , let P c ∈ R M ( n, c ) × M ( n, c ) be symmetric posi- tive definite. If B T c, i P c + P c B c,i < 0 (17) for all i ∈ { 1 , · · · , N } where B c, i := W + c A c, i W c , (18) then V c ( x ) = y c ( x ) T P y c ( x ) is a homogeneous polynomial L yapunov function for (1) of or der 2 c . The proof of this result comes from the fact that A c ( t ) has an inherent in variant subspace, resulting from its construc- tion. As trajectories of d dt ( ⊗ c x ) = A c ( t )( ⊗ c x ) are kno wn to begin in this subspace, we can encode the search for meta-L yapunov functions for the system (1) as a search for quadratic L yapunov functions for the reduced order system ˙ y c = B c ( t ) y c (19) where B c ( t ) ∈ { B c, 1 , · · · , 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 ) . Howe ver , each ordering will induce a unique W c , and thus the resulting L yapunov conditions will always be the same, regardless of the chosen ordering. Moreov er , the constraints giv en by (17) are equiv alent to the constraints giv en 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 , gi ven by (16), can be captured in closed form. Proposition 2. Consider the system (1) and let n = 2 . In this case we have M ( n, c ) = c + 1 . Additionally , for a positive integ er k ∈ N ≥ 0 , let 0 k ∈ R k denote a vector populated with zeros. If y c ( x ) ∈ R c +1 conforms to the ordering y c ( x ) =  x c 1 x c − 1 1 x 2 · · · x c 2  T , then we have ⊗ c x = W c y c ( x ) , wher e for an inte ger k ∈ N ≥ 1 we define W k r ecursively by W 1 = I 2 W k =  W k − 1 0 2 k − 1 0 2 k − 1 W k − 1  k ≥ 2 . (20) In the case when n > 2 , it is generally difficult to order the c th order monomials of x in an intuiti ve 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 ) . Howe ver , W c can always be solved for using (16) once y c ( x ) has been chosen. I V . R E L A T I O N T O H O M O G E N E O U S P O LY N O M I A L L Y A P U N OV F U N C T I O N S T raditionally , the search for a polynomial L yapunov func- tions 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 , · · · , g r such that p ( x ) = r X i =1 g i ( x ) 2 . The search for a sum-of-squares polynomial V ( x ) , satis- fying (2), is kno wn to be a conv ex optimization problem, computable by solving a semidefinite program [13]. Many efficient solvers exist to handle such problems [6], [16]. W e next show that the existence of quadratic L yapunov functions for the hierarch y of dynamical systems (13) guarantees the existence of a homogeneous sum-of-squares polynomial L yapunov functions for (1), and vice versa. In this sense, all homogeneous sum-of-squares polynomial L yapunov functions can be thought of as quadratic L yapunov functions for a related hierarchy of differential equations. Moreov er , one can encode the search for high-order sum- of-squares polynomial L yapunov functions, which certify the stability of (1), as a search for quadratic L yapunov functions for the related system (13). Calculating sum-of- squares polynomial L yapunov 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. Ther e exists a P c ∈ R M ( n, c ) × M ( n, c ) satisfying (17) for some positive integ er c ∈ N ≥ 1 , if and only if ther e exists a homogeneous sum-of-squares polynomial L yapunov function V c ( x ) of de gr ee 2 c for the system (1) . Pr oof. A sum-of-squares polynomial that is homogeneous in the entries x and of order 2 c will take the form p ( x ) = y c ( x ) T Z y c ( x ) , where Z ∈ R M ( n,c ) × M ( n,c ) is symmetric, and y c ( x ) and M ( n, c ) retain their definitions from Theorem 3. From Theorem 3, we hav e that if P c ∈ R M ( n, c ) × M ( n, c ) satisfies (11) for some positi ve inte ger c ∈ N ≥ 1 , then we ha ve that V c ( x ) = y c ( x ) T P y c ( x ) is a homogeneous polynomial L yapunov function for (21) and, moreover , V c ( x ) is a sum- of-squares. T o prove the conv erse, we note that if p ( x ) = y c ( x ) T Z y c ( x ) is a homogeneous sum-of-squares polynomial L yapunov function for (21) then Z > 0 and ˙ p ( x ) < 0 for all x ∈ R n . From the dynamics of y c ( x ) , given as (19), we hav e B T c, i Z + Z B c,i < 0 for all i ∈ { 1 , · · · , N } . Therefore P c = Z solves (17). V . N U M E R I C A L E X A M P L E In this section, we provide an example case and prove the stability of a switched linear system using a meta- L yapunov function based approach. An algorithm is pro vided for generating homogeneous polynomial L yapunov functions for switched systems, which follows the procedure detailed in Theorem 3; this algorithm is specifically written for imple- mentation with CVX, a con vex optimization toolbox made for use with MA TLAB [17]. W e also provide a comparison to a similar search for homogeneous polynomial L yapunov functions that was implemented using SOST OOLS, a sum- of-squares optimization toolbox made for use with MA TLAB [18]. Experimental results are provided from MA TLAB 2019b, which w as run on a 2017 Macbook Pro laptop. A. Pr oblem F ormulation W e consider the linear time-variant system ˙ x = A ( t ) x A ( t ) ∈ { A 1 , A 2 } (21) A 1 =  − . 5 . 5 − . 5 − . 5  A 2 =  − 2 . 5 2 . 5 − 2 . 5 1 . 5  . In the follo wing, we go about showing that (21) is stable. This is done, at first, through the computation of a quadratic L yapunov function V 1 ( x ) = x T P x , which satisfies (2), and then through the computation of higher-order homogeneous polynomial L yapunov 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 L yapunov function V c ( x ) that certifies the stability of (21), then the infinite-time system trajectory is constrained to stay inside x ( t ) ∈ { x ∈ R n | V c ( x ) ≤ V c ( x 0 ) } (22) for all t ≥ 0 . For this reason, we select V c ( x ) as the minimizers of a suitable objective function, as to shrink the resulting in variant re gion deri ved through (22). In what follows, we additionally show that computing higher -order meta-L yapunov functions allows one to characterise tighter in variant sets by (22), even when the same objectiv e function is used in each computation. B. Identifying Meta-Lyapuno v Functions W e search for meta-L yapunov functions for (21) using a semidefinite program. Specifically , when searching for a homogeneous L yapunov function of order 2 c , we first calculate B c, 1 and B c, 2 using equations (12), (18) and (20), and then we search for a symmetric positi ve-definite matrix P c ∈ R ( c +1) × ( c +1) that satisfies (17). Such a matrix iden- tifies V c ( x ) = y c ( x ) T P c y c ( x ) as a polynomial L yapunov function for (21), which is homogeneous in the entries of x and of order 2 c . W e implement the aforementioned procedure with Algorithm 1, which specifically relies on CVX, a con vex optimization toolbox built for use with MA TLAB [17], [19]. Algorithm 1 takes as inputs the system parameters A 1 and A 2 , and a positiv e integer c , and returns a matrix P c , in the case that one exists, which satisfies (11) at the c th lev el. Note that Algorithm 1 computes the solution to a feasi- bility problem, rather than an optimization problem; that is, while Algorithm 1 searches for a P c that satisfies the meta- L yapunov constraint (17), this solution is computed without referencing any objective function. Note howe ver , that in the instance that multiple feasible solutions exist, it is preferable to choose P c such that the suble vel sets of the resulting homogeneous L yapunov function V c ( x ) = y c ( x ) T 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 P c as the solution to an optimisation problem, rather than a feasibility problem. Little is known, in general about ho w one can relate the parameters of a polynomial to the volume of its sublev el sets. In our case as well, it is difficult to associate a metric of optimality with the a feasible solution to the meta-L yapunov constraints (17). Through e xperimentation, we have generally found that it is preferable to compute numerous solutions using dif ferent objecti ve functions, and then compute an in variant re gion as the intersection of their respecti ve sublevel Algorithm 1 Computing Meta-L yapunov Functions input : A 1 , A 2 ∈ R 2 × 2 from (1). c ∈ N ≥ 1 . output : P c ∈ R ( c +1) × ( c +1) satisfying (17). 1: function M E T A L Y A P U N O V ( A 1 , A 2 , c ) 2: Initialize: Compute A c, 1 and A c, 2 by (12) 3: Compute W c by (20) 4: B c, 1 ← W + c A c, 1 W c 5: B c, 2 ← W + c A c, 2 W c 6: cvx begin sdp 7: variable P c ( c + 1 , c + 1) semidefinite 8: 0 > B T c, 1 P c + P c B c, 1 9: 0 > B T c, 2 P c + P c B c, 2 10: P c > I n 11: %% Possibly Insert Objecti ve Function 12: cvx end 13: if Program feasible then 14: retur n P c 15: else 16: retur n ‘infeasible’ 17: end function sets. Specifically we recommend using either using the objectiv e function 11: minimize P c (1 , 1) which minimises the coefficient on x 2 c 1 in the resultant L yapunov function V c ( x ) , or 11: minimize P c ( c + 1 , c + 1) which minimises the coefficient on x 2 c 2 . These objecti ve functions are provided in psuedocode, such that the y can easily be inserted in Algorithm 1 at line 11. C. Numerical Results and Comparison with SOST OOLS W e no w return to the example system (21), and compute feasible meta-L yapunov functions with Algorithm 1. Addi- tionally , we compute an over approximation of the infinite time reachable set of (21) when be ginning from the initial conditions x (0) = [1 , 0] T . As discussed in the preceding, we compute these in variant sets by implementing Algorithm 1, while attempting to minimize P (1 , 1) , i.e. the coef ficient on x 2 c 1 ; see Algorithm 1, Line 11. This procedure was computed in MA TLAB 2019b using CVX. In the case of this example, Algorithm 1 was computed for c ∈ { 1 , 2 , · · · , 13 } , thus generating homogeneous polynomial L yapunov functions for all ev en orders between 2 and 26. These L yapunov functions were then used to calculate in variant regions of the state space using (22); see Figure 1. Note that as the order of the meta-L yapunov function increases, the deriv ed inv ariant sets shrink in volume. Fur - ther , certain higher-order the meta-L yapunov functions were − 1 0 1 − 2 − 1 0 1 x 1 x 2 Reachable Set 2nd order 10th order 16th order 26th order Initial State Fig. 1: Simulated system response of (21). When starting from x 0 = [1 , 0] T , the system can only reach the re gion shown in light yellow , which was computed via simulation. The equipotentials of high-order homogeneous L yapunov functions are also shown. Specifically , the dark blue, light blue, orange and red regions represent inv ariant sets calcu- lated using quadratic, 10 th -order , 16 th -order and 26 th -order meta-L yapunov functions, respectiv ely . The in variance of these regions is sho wn by (22). shown to ha ve non-con ve x sublev el sets. W e provide the number of solv er iterations for each experiment, as well as the computations times, in Figure 2. W e ne xt compare the meta-L yapunov function based method to the more traditional sum-of-squares based ap- proach for calculating inv ariant regions. Specifically , we search for high-order homogeneous polynomial L yapunov functions for the system (21) using SOSTOOLS, MA TLAB’ s sum-of-squares toolbox [18]. W e similarly implement SOS- TOOLS with the solver SDPT3 and attempt to generate homogeneous polynomial L yapunov functions for the system (21) while minimizing the coefficient on x 2 c 1 . As was the case pre viously , we provide the number of solv er iterations for each experiment, as well as the computations times (See Figure 2). In the experiment, SOSTOOLS was only able to generate homogeneous polynomial L yapunov functions of order 10 or below; a solver error was returned during each search for more complex L yapunov function. In contrast, Algorithm 1, when implemented through CVX, was able to generate up to 26 th -order polynomial L yapunov functions while minimizing the same objectiv e function. W e attribute this discrepancy L yapuno v function V c ( x ) , or 11: minimize P c ( c +1 ,c + 1) which minimises the coef ficient on x 2 c 2 . These objecti v e functions are pro vided in psuedocode, such that the y can easily be inserted in Algorithm 1 at line 11. C. Numerical Results and Comparison with SOStools W e no w return to the e xample system (21), and compute feasible meta-L yapuno v functions with Algorithm 1. Addi- tionally , we compute an o v er approximation of the infinite time reachable set of ( 2 1) , when be ginning from tw o possible initial conditions x ( 0) 2 ⇢ 1 0  ,  0 1  . As discussed in the preceding, we compute these in v ariant sets by implementing Algorithm (1), while attempting to minimize P (1 , 1) , i.e. the c o e f ficient on x 2 c 1 ; see Algorithm 1, Line 11. This procedure w as computed in MA TLAB 2019b using CVX, implemented with the solv er SDPT3 [19]. In the case of this e xample, Algorithm 1 w as computed for c 2 { 1 , 2 , ·· · , 13 } . As such, homogeneous polynomial L ya- puno v functions were generated for all e v en orders between 2 and 26; these polynomials were then used to calculate in v ariant re gions of the state space (see Figure 2). Note that as the order of the meta-L yapuno v function increases, the deri v ed i n v ar iant sets shrink in v olume. Further , certain higher -order the meta-L yapuno v functions were sho wn to ha v e non-con v e x suble v el sets. W e pro vide the number of solv er iterations for each e xperiment, as well as the compu- tations times, in Figure 1. W e ne xt compare the meta-L yapuno v function based method for o v er -approximating reachable sets t o the more traditional sum-of-squares based approximation method. Specifically , we search for high-order homogeneous poly- nomial L yapuno v functions for the system (21) using SOS- T OOLS, MA TLAB ’ s sum-of-squares toolbox [20]. W e simi- larly implement SOST OOLS with the solv er SDPT3 and at- tempt to generate homogeneous polynomial L yapuno v func- tions for the system (21) while minimizing the coef ficient on x 2 c 1 . As w as the case pre viously , we pro vide the number of solv er iterations for each e xperiment, as well as the computations times, in Figure 1. In the e xperiment, SOST OOLS w as only able to generate homogeneous polynomial L yapuno v functions of order 10 or belo w; a solv er error w as returned during each search for more comple x L yapuno v function. In contrast, Algorithm 1, when i mplemented through CVX, w as able to generate up to 26 th -order polynomial L yapuno v functions while minimizing the same objecti v e function. W e attrib ute this discrepanc y to the f act that man y of the steps required in a traditional sum-of-squares based search optimization are not required by Algorithm 1. F or e xample, Algorithm 1 does not compute the time rate of change of the monomials in x of order 2 c ; that is, Algorithm 1 be gins with a closed L yapuno v function V c ( x ) , or 11: minimize P c ( c +1 ,c + 1) which minimises the coef ficient on x 2 c 2 . These objecti v e functions are pro vided in psuedocode, such that the y can easily be inserted in Algorithm 1 at line 11. C. Numerical Results and Comparison with SOStools W e no w return to the e xample system (21), and compute feasible meta-L yapuno v functions with Algorithm 1. Addi- tionally , we compute an o v er approximation of the infinite time reachable sets of tw o possible initial conditions x ( 0) 2 ⇢ 1 0  ,  0 1  . As discussed in the preceding, we compute these in v ariant sets by implementing Algorithm (1), while attempting to minimize P (1 , 1) , i.e. the c o e f ficient on x 2 c 1 ; see Algorithm 1, Line 11. This procedure w as computed in MA TLAB 2019b using CVX, implemented with the solv er SDPT3 [16]. In the case of this e xample, Algorithm 1 w as computed for c 2 { 1 , 2 , ·· · , 13 } . As such, homogeneous polynomial L ya- puno v functions were generated for all e v en orders between 2 and 26; these polynomials were then used to calculate in v ariant re gions of the state space (see Figure 1). Note that as the order of the meta-L yapuno v function increases, the deri v ed i n v ar iant sets shrink in v olume. Further , certain higher -order the meta-L yapuno v functions were sho wn to ha v e non-con v e x suble v el sets. W e pro vide the number of solv er iterations for each e xperiment,as well as the compu- tations times, in T able —. W e ne xt compare the meta-L yapuno v function based method for o v er -approximating reachable sets t o the more traditional sum-of-squares based approximation method. Specifically , we search for high-order homogeneous poly- nomial L yapuno v functions for the system (21) using SOS- T OOLS, MA TLAB ’ s sum-of-squares toolbox [17]. W e simi- larly implement SOST OOLS with the solv er SDPT3 and at- tempt to generate homogeneous polynomial L yapuno v func- tions for the system (21) while minimizing the coef ficient on x 2 c 1 . In the e xperiment, SOST OOLS w as only able to generate homogeneous polynomial L yapuno v functions of order 10 or belo w; a solv er error w as returned during each search for more comple x L yapuno v function. In contrast, Algorithm 1, when im plemented through CVX, w as able to generate up to 26 th -order polynomial L yapuno v functions whi le minimizing the same objecti v e function. W e attrib ute this discrepanc y to the f act that man y of the steps required in a traditional sum-of-squares based search optimization are not required by Algorithm 1. F or e xample, Algorithm 1 does not compute the time rate of change of the monomials in x of order 2 c ; that is, Algorithm 1 be gins with a clos ed from representation of ˙ y c ( x ) , which is encoded in the matrices B c, 1 and B c, 2 . In contrast, SOS tools must compute ˙ y c ( x ) online as a sum of squares of lo wer -order monomials. F or these reason, one can e xpect SOStools to preform less ef ficiently . c MetaL yapuno v SOST OOLS Iterations T ime (s) Iterations T imes (s) 1 9 1.510 11 0.396 2 10 1.580 14 0.416 5 20 1.680 25 0.562 6 21 1.730 33 N/ A 10 36 2.400 33 N/ A 12 40 2.480 31 N/ A 13 54 3.300 31 N/ A 14 36 N/A 31 N/ A VI. C ON CLU SION This w ork addresses the problem of searching for homo- geneous polynomial L yapuno v functions for stable switched linear systems. An equi v alence is sho wn between polyno- mial L yapuno v functions for switched linear systems and quadratic L yapuno v functions for a related hierarch y of L ya- puno v dif ferential equations. A computationally competiti v e algorithm is presented for generating high-order homoge- neous polynomial L yapuno v functions in this manner . R EF E REN CES [1] D. Liberzon, Switc hing in Systems and Contr ol . Systems & Control: F oundations & Applications, Birkh ¨ auser Boston, 2003. [2] E. Feron, Quadr atic stabilizability of switc hed systems via state and output feedbac k . Center for Intelligent Control Systems, 1996. [3] R. Goebel, R. G. Sanfelice, and A. R. T eel, “Hybrid dynamical systems, ” IEEE Contr ol Systems Ma gazine , v ol. 29, pp. 28–93, April 2009. [4] R. W . Brock ett, Hybrid Models for Motion Contr ol Systems , pp. 29–53. Boston, MA: Birkh ¨ auser Boston, 1993. [5] S. Bo yd, L. El Ghaoui, E. Feron, and V . Balakrishnan, Linear Matrix Inequalities in System and Contr ol Theory , v ol. 15 of Studies in Applied Mathematics . Philadelphia, P A: SIAM, June 1994. [6] B. R. Barmish, “Necessary and suf ficient conditions for quadratic stabilizability of an uncertain system, ” J ournal of Optimization Theory and Applications , v ol. 46, pp. 399–408, Aug 1985. [7] Y . Y oon, C. Klett , and E. Feron, “Bounding the stat e co v ariance matrix for a randomly switching linear system with noise, ” arXiv pr eprint , 2019. [8] M. Johansson and A. Rantzer , “Computation of piece wise quadratic lyapuno v functions for h ybrid systems, ” in 1997 Eur opean Contr ol Confer ence (ECC) , pp. 2005–2010, IEEE, 1997. [9] P . A. P arrilo, Structur ed semidefinite pr o gr ams and semialg ebr aic g e- ometry methods in r ob ustness and optimization . PhD thesis, California Institute of T echnology , 2000. [10] P . A. P arrilo, “Semidefinite programming relaxations for semialgebraic problems, ” Mathematical pr o gr amming , v ol. 96, no. 2, pp. 293–320, 2003. [11] P . Mason, U. Boscain, and Y . Chitour , “Common polynomial lyapuno v functions for linear switched systems, ” SIAM J ournal on Contr ol and Optimization , v ol. 45, 04 2004. [12] A. A. Ahmadi and P . A. P arrilo, “Con v erse results on e xistence of sum of squares lyapuno v functions, ” in 2011 50th IEEE Confer ence on Decision and Contr ol and Eur opean Contr ol Confer ence , pp. 6516– 6521, Dec 2011. [13] J. F . Sturm, “Using sedumi 1.02, a m atlab toolbox for optimization o v er symmetric cones, ” Optimization Methods and Softwar e , v ol. 11, no. 1-4, pp. 625–653, 1999. [14] M. Grant and S. Bo yd, “CVX: Matlab softw are for disciplined con v e x programming, v ersion 2.1. ” http://cvxr .com/cvx, Mar . 2014. Fig. 1: Comparing Algorithm 1 to a search for homogeneous L yapuno v functions using SOST OOLs. In each e xperiment, from representation of ˙ y c ( x ) , which is encoded in the matrices B c, 1 and B c, 2 . SOS tools must compute ˙ y c ( x ) online as a sum of squares of lo wer -order monomials and, for this reason, one can e xpect SOST OOLS to preform less ef ficiently . 2c MetaL yapunov SOSTOOLS Iterations T ime (s) Iterations Times (s) 2 9 1.510 11 0.396 4 10 1.580 14 0.416 10 20 1.680 25 0.562 12 21 1.730 33 N/A 20 36 2.400 33 N/A 24 40 2.480 31 N/A 26 54 3.300 31 N/A 28 36 N/A 31 N/A VI. C ON CLU SION This w ork addresses the problem of searching for homo- geneous polynomial L yapuno v functions for stable switched linear systems. An equi v alence is sho wn between polyno- mial L yapuno v functions for switched linear systems and quadratic L yapuno v functions for a related hierarch y of L ya- puno v dif ferential equations. A computationally competiti v e algorithm is presented for generating high-order homoge- neous polynomial L yapuno v functions in this manner . R EF E REN CES [1] D. Liberzon, Switc hing in Systems and Contr ol . Systems & Control: F oundations & Applications, Birkh ¨ auser Boston, 2003. [2] E. Feron, Quadr atic stabilizability of switc hed systems via state and output feedbac k . Center for Intelligent Control Systems, 1996. [3] R. Goebel, R. G. Sanfelice, and A. R. T eel, “Hybrid dynamical systems, ” IEEE Contr ol Systems Ma gazine , v ol. 29, pp. 28–93, April 2009. [4] R. W . Brock ett, Hybrid Models for Motion Contr ol Systems , pp. 29–53. Boston, MA: Birkh ¨ auser Boston, 1993. [5] V . A. Y akubo vich, “The solution of some matrix inequalities encoun- tered in automatic control theory , ” Societ Math Dokl. , v ol. 143, no. 5, pp. 652–656, 1964. [6] S. Bo yd, L. El Ghaoui, E. Feron, and V . Balakrishnan, Linear Matrix Inequalities in System and Contr ol Theory , v ol. 15 of Studies in Applied Mathematics . Philadelphia, P A: SIAM, June 1994. [7] B. R. Barmish, “Necessary and suf ficient conditions for quadratic stabilizability of an uncertain system, ” J ournal of Optimization Theory and Applications , v ol. 46, pp. 399–408, Aug 1985. Fig. 2: Comparing Algorithm 1 to a search for homogeneous L yapunov functions using SOSTOOLS. Algorithm 1 is used to compute homogeneous L yapunov functions of orders 2, 4, 10, 12, 20, 24, and 26 for the system (21). SOST OOLS, how- ev er, is only able to correctly generate L yapunov functions of order 10 and belo w . The computation time and number of solver iterations are provided for each experiment. The symbol N/A is used when the solv er is unable to find a homogeneous L yapunov function of a certain order; in this case, the number of solver iterations which were preformed before failure is also provided. 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 2 c ; that is, Algorithm 1 begins with a closed from representation of ˙ y c ( x ) , which is encoded in the matrices B c, 1 and B c, 2 . SOSTOOLS must compute ˙ y c ( x ) online as a sum of squares of lo wer-order monomials and, for this reason, one can expect SOSTOOLS to preform less ef ficiently . Despite this, in the experiments where SOSTOOLS was able to correctly generate homogeneous L yapunov functions, we found that SOSTOOLS preformed f aster in computation that the meta-L yapunov search implemented with CVX; howe ver , it did take SDPT3 more solver iterations to generate the solution which minimized the objectiv e function when implemented through SOST OOLS (Figure 2). V I . C O N C L U S I O N This w ork addresses the problem of searching for homo- geneous polynomial L yapunov functions for stable switched linear systems. An equi v alence is sho wn between polyno- mial L yapunov functions for switched linear systems and quadratic L yapunov functions for a related hierarchy of L ya- punov differential equations. A computationally competitiv e algorithm is presented for generating high-order homoge- neous polynomial L yapunov functions in this manner . R E F E R E N C E S [1] D. Liberzon, Switching in Systems and Contr ol . Systems & Control: Foundations & Applications, Birkh ¨ auser Boston, 2003. [2] E. Feron, Quadratic stabilizability of switched systems via state and output feedback . Center for Intelligent Control Systems, 1996. [3] R. Goebel, R. G. Sanfelice, and A. R. T eel, “Hybrid dynamical systems, ” IEEE Control Systems Magazine , vol. 29, pp. 28–93, April 2009. [4] R. W . Brockett, Hybrid Models for Motion Contr ol Systems , pp. 29–53. Boston, MA: Birkh ¨ auser Boston, 1993. [5] V . A. Y akubovich, “The solution of some matrix inequalities encoun- tered in automatic control theory , ” Societ Math Dokl. , vol. 143, no. 5, pp. 652–656, 1964. [6] S. Boyd, L. El Ghaoui, E. Feron, and V . Balakrishnan, Linear Matrix Inequalities in System and Contr ol Theory , v ol. 15 of Studies in Applied Mathematics . Philadelphia, P A: SIAM, June 1994. [7] B. R. Barmish, “Necessary and sufficient conditions for quadratic stabilizability of an uncertain system, ” Journal of Optimization Theory and Applications , vol. 46, pp. 399–408, Aug 1985. [8] Y . Y oon, C. Klett, and E. Feron, “Bounding the state covariance matrix for a randomly switching linear system with noise, ” arXiv pr eprint arXiv:1905.09427 , 2019. [9] V . Y akubovich, “The method of matrix inequalities in the stability theory of nonlinear control systems, ” Automation and Remote Control , vol. 26, pp. 577–592, 1965. [10] V . A. Y akubovich, “Frequency conditions for the existence of abso- lutely stable periodic and almost periodic limiting regimes of control systems with many nonstationary elements, ” IF A C W orld Congr ess , 1966. [11] M. Johansson and A. Rantzer, “Computation of piecewise quadratic lyapunov functions for hybrid systems, ” in 1997 European Control Confer ence (ECC) , pp. 2005–2010, IEEE, 1997. [12] P . A. Parrilo, Structur ed semidefinite pr ogr ams and semialgebraic ge- ometry methods in rob ustness and optimization . PhD thesis, California Institute of T echnology , 2000. [13] P . A. Parrilo, “Semidefinite programming relaxations for semialgebraic problems, ” Mathematical pr ogramming , vol. 96, no. 2, pp. 293–320, 2003. [14] P . Mason, U. Boscain, and Y . Chitour, “Common polynomial lyapunov functions for linear switched systems, ” SIAM journal on control and optimization , vol. 45, no. 1, pp. 226–245, 2006. [15] A. A. Ahmadi and P . A. Parrilo, “Conv erse results on existence of sum of squares Lyapuno v functions, ” in 2011 50th IEEE Conference on Decision and Contr ol and Eur opean Contr ol Confer ence , pp. 6516– 6521, Dec 2011. [16] J. F . Sturm, “Using sedumi 1.02, a matlab toolbox for optimization over symmetric cones, ” Optimization Methods and Softwar e , vol. 11, no. 1-4, pp. 625–653, 1999. [17] M. Grant and S. Boyd, “CVX: Matlab software for disciplined conv ex programming, version 2.1, ” Mar . 2014. [18] G. V . S. P . P . S. A. Papachristodoulou, J. Anderson and P . A. Parrilo, SOSTOOLS: Sum of squar es optimization toolbox for MATLAB . [19] M. Grant and S. Boyd, “Graph implementations for nonsmooth conv ex programs, ” in Recent Advances in Learning and Contr ol , Lecture Notes in Control and Information Sciences, pp. 95–110, Springer- V erlag Limited, 2008.