{'fraction_non_alphanumeric': 0.04325668623613829, 'fraction_numerical': 0.0136986301369863, 'mean_word_length': 4.290983606557377, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Numerical Method to Compute Stability Margins of Switching Linear Systems

Stability margins for linear time-varying (LTV) and switched-linear systems are traditionally computed via quadratic Lyapunov functions, and these functions certify the stability of the system under study. In this work, we show how the more general class of homogeneous polynomial Lyapunov functions is used to compute stability margins with reduced conservatism, and we show how these Lyapunov functions aid in the search for periodic trajectories for marginally stable LTV systems. Our work is premised on the recent observation that the search for a homogeneous polynomial Lyapunov function for some LTV systems is easily encoded as the search for a quadratic Lyapunov function for a related LTV system, and our main contribution is an intuitive algorithm for generating upper and lower bounds on the system's stability margin. We show also how the worst-case switching scheme - which draws an LTV system closest to a periodic orbit - is generated. Three numerical examples are provided to aid the reader and demonstrate the theoretical contributions of the work.


I. INTRODUCTION
Switched-linear and linear time-varying (LTV) systems are two important and widely studied classes of systems [1]. They are useful abstractions for studying systems with uncertain parameters, and they can be used to represent mode-switched systems when the switching scheme may be unknown. In all cases, computing stability margins is an important exercise, with practical benefits for safety critical systems [2]. In the following, we study both switched-linear and LTV systems in an equivalent manner and propose methods to compute upper and lower bounds on a system's stability margin.
We study the stability of LTV systems whose dynamics are influenced by an uncertain system parameter. We characterize the system's robustness with respect to this parameter by computing an upper and a lower bound such that the system is always stable when the parameter is less than the lower bound and such that a diverging trajectory exists when the parameter is greater than the upper bound. As will be shown, these bounds hint at the existence of a specific parameter value that induces periodic orbits for LTV systems but does not induce diverging trajectories. For this reason, we approach the problem by studying methods for finding periodic trajectories in marginally LTV systems.
Stability margins are generally impossible to compute exactly for all but linear time-invariant systems, though much work has been done on estimating the margin through the computation of upper and lower bounds [3] [4]. For timevarying systems, while quadratic stability has been studied in depth since the late 1940s by Lur'e, Yakubovitch, and others [5], exact approaches to the computation of stability margins are more recent and include [6] and [7]. However, most computationally tractable analyses are generally conservative and typically give results which may significantly under-or over-approximate the true stability margin.
In order to gain an understanding of the conservatism of these types of bounding methods, we compute a lessconservative bound by searching for a trajectory which becomes unstable when the system is allowed to operate near its stability margin. Such a trajectory is referred to as a counterexample and is a powerful tool for numerically approximating the true stability margin. Generating a counter-example can be challenging. In many cases, simulations will show that a time-varying system operating near or even beyond its stability margin is always stable. Work by Parrilo and Jadbabaie investigates computation of worst-case switching sequences for discrete time systems [8]. Our work takes inspiration from this approach when computing counterexamples for continuous LTV systems.
In the case where the system matrix remains Hurwitz for all time, the examples studied indicate that instability occurs through a switching sequence which leads the system through a Hopf bifurcation as the uncertain parameter is adjusted. A theoretical result which confirms this observation is not achieved, but the observation nonetheless provides the intuition behind the algorithm for computing a counterexample and therefore an upper bound on the stability margin. In classical linear systems analysis, a system which depends on a single parameter moves from stable to unstable when, as the parameter is adjusted, the poles on the root locus cross through the origin or through the imaginary axis. In the latter case, instability is achieved via a Hopf bifurcation. For the LTV system currently under study, the system matrix can remain Hurwitz so that all poles stay in the left-half plane for all time, but the time-dependent variation of the uncertain parameter can induce a limit cycle.
Our approach to bounding stability margins is aided by the recently presented hierarchy of polynomial Lyapunov functions for switched-linear systems [9]. In that work, the search for homogeneous polynomial Lyapunov functions for LTV systems is recast as a search for quadratic Lyapunov functions for a related hierarchy of time-varying Lyapunov differential equations; thus, this is an elegant and equivalent alternative to sum-of-squares programming techniques for computing Lyapunov functions [10]. Intuitive procedures following from these observations have been shown to provide higher-fidelity performance analyses [11] [12]. The Lyapunov functions returned by these analyses also allow us to improve stability margin analysis by providing valuable information about how to compute counter-examples. A numerical example in this paper shows how a Lyapunov function for an LTV system can be used to produce a trajectory which approaches the bound computed in [11].
This paper is structured as follows. In Section II, we review mathematical preliminaries and the Kronecker productbased formulation of the LTV system hierarchy. Section III defines the system under study and defines the problems addressed in this work. In Section IV, the hierarchy of quadratic Lyapunov functions is used to compute a lower bound on the system stability margin. Sections V and VI use these Lyapunov functions to develop the algorithm for producing a counter-example and an upper bound on the stability margin. Numerical examples are presented in section VII.

II. NOTATION AND PRELIMINARIES
We denote by S n ++ ⊂ R n×n the set of symmetric positive definite n × n matrices. We denote by I n the n × n identity matrix, and we denote by 0 n ∈ R n the zero vector in R n . Given M ∈ R n×m and integer i ≥ 1, we denote by ⊗ i M ∈ R n i ×m i the i th -Kronecker power of M , as defined recursively by We study, in the following, switched linear systems, as iṅ with state x ∈ R n and where A(t) ∈ S ⊆ R n×n evolves for all time inside a finite set of switched modes We are chiefly interested in studying the stability properties of (2).
Definition 1 (Stability). The system (2) is stable if for all initial conditions x(0) ∈ R n , we have that lim t→∞ x(t) = 0 n for all A(t) satisfying (3).
The system (2) is unstable when (2) is not stable.
Definition 2 (Quadratic Stability). For P ∈ S n ++ , the function V (x) = x T P x is a quadratic Lyapunov function for (2) if When (2) is linear time-invariant-equivalently, when k = 1-the system (2) is stable if and only if there exists a quadratic Lyapunov function for (2). This is, however, not true in the general setting of (2); indeed, there exists stable switched systems for which no quadratic Lyapunov function exists [1]. Nonetheless, the existence of such a quadratic Lyapunov function for (2) guarantees the stability of the system.
It was recently shown in [9] how a related hierarchy of switched linear systems is used to assess the stability of (2). Consider the following infinite hierarchy of switched linear systems: is the state of the i th -level system H i and i ≥ 1. Each system H i in the hierarchy is switched linear, as in (2).
The hierarchy (5) is best understood by looking at the first and second level systems H 1 and H 2 . The system H 1 is equivalent to the switched linear system (2), and the system H 2 is nothing but the vectorized version of the Lyapunov differential equatioṅ where A(t) maintains its definition from (2) and where the state of (6) is a matrix X ∈ R n×n . Moreover, it is shown in [9] that if x(t) is a solution to (2), then ξ i (t) = (⊗ i x(t)) is a solution to H i . It follows, therefore, that if the i th -level H i is stable, then the system (2) is stable as well.

Proposition 1 ( [9]
). If the i th -level system H i is quadratically stable, i.e. if there exists a P ∈ S n i ++ such that then the system (2) is stable and is a homogeneous polynomial Lyapunov function for (2) of order 2i.
Note that the system (2) can be stable even when H i is not quadratically stable for a given i ≥ 1. Nonetheless, the stability of (2) implies that there must exist an i ≥ 1 such that H i is quadratically stable [9].
III. PROBLEM STATEMENT In the following, we consider linear time-varying (LTV) systemsẋ for A, A 0 ∈ R n×n and time-varying scalar ∆(t) ≥ 0. We refer to A as the nominal dynamics, and the term ∆(t)A 0 describes a time-varying perturbation. We assume that A is Hurwitz so that (9) is stable when ∆(t) ≡ 0. (2) is a generalization of (9); in particular, if S in (3) is taken as

Remark 1. The switched linear system
then the switched system (2) is stable if and only if the linear time varying system (9) is stable for all ∆(t) that satisfy The stability margin for (9) is the uniqueδ ≥ 0 such that (9) is stable with respect to all δ ∈ [0,δ) and such that (9) is not stable with respect tô δ.
The system (9) is guaranteed to be stable with respect to all δ ∈ [0,δ), and for this reason, we begin by studying methods for under-approximating the stability margin for (9). Problem 1. Under-approximate the stability marginδ for (9) by finding the largest possible δ for which a Lyapunov function can be computed which certifies that (9) is stable with respsect to δ. We denote this lower bound onδ as δ.
Numerous methods exist for computing δ [5], but these methods generally provide conservative estimates of the stability marginδ. In Section IV, we present an intuitive procedure for solving Problem 1 whereby the stability margin is under-approximated using a homogeneous polynomial Lyapunov function for (9). These Lyapunov functions are computed using the hierarchy of switched systems (5), and we show how conservatism is mitigated in our approach when high-order Lyapunov functions are considered.
Next, consider the case that A + δA 0 remains Hurwitz and assume that, for δ =δ + ε, there exists an ε > 0 and a switching function ∆(t) which produces a diverging trajectory. Based on the examples studied, we observe that, when ε = 0, there exists a ∆(t) that can induce a periodic trajectory but not a diverging trajectory. Furthermore, for any ε > 0 there exists a ∆(t) which can produce a diverging trajectory. This observation has not been proven; however, it inspires the development of a useful algorithm for bounding the stability margin from above and finding a periodic trajectory.
Definition 4. A switching signal ∆(t) ∈ [0, δ] which can cause a trajectory of (9) to become periodic when δ =δ is called a worst-case switching function and is denoted as ∆ w (t; δ). The system matrix A(t) produced by the worstcase switching function is denoted by Remark 2. It is not proven that ∆ w (t; δ) always exists, but in practice it has been approximated for every example LTV system studied. It is also not necessarily unique, but that is not important since any ∆ w (t; δ) will produce a periodic trajectory when δ =δ, and it is the approximation of thisδ which is of primary interest.
Even for δ <δ, a worst-case switching function ∆ w (t; δ) is very useful for system analysis as it can be used to slow a trajectory's path to the origin. See Example 3 in Section VII.
While ∆ w (t; δ) can not be computed exactly, it can be closely approximated using a Lyapunov function for (9), and a higher-order Lyapunov function will produce a better approximation. We observe that the transition from stability to instability occurs through a switching sequence which results in a periodic orbit about the surface of a level set of a Lyapunov function.
Problem 2. Given A, A 0 , δ and a Lyapunov function parameterized by P from the solution to Problem 1, approximate ∆ w (t; δ) and A w (t; δ). These approximations are denoted as ∆ w (t; δ, P ) and A w (t; δ, P ). Often, the order of the Lyapunov function used to generate the switching function is included as a subscript, as in ∆ w (t; δ, P 2i ).
Once an approximate worst-case switching function is found, the smallest value of δ for which ∆ w (t; δ, P ) produces a periodic trajectory can be computed.
Problem 3. Find the smallest value of δ such that the system (9) with A w (t; δ, P ) produces a periodic trajectory. In this case, δ approximates the lowest upper bound on the stability marginδ. This upper bound is denoted as δ.

IV. UNDER-APPROXIMATING THE STABILITY MARGIN
We begin by addressing Problem 1. The system (9) is stable if and only if there exists a homogeneous sum-ofsquares polynomial Lyapunov function certifying its stability [13] [14]. Thus, we present an iterative approach for underapproximatingδ, and this approach relies on the hierarchy of switched systems (5).
Proposition 2 is a direct result of Proposition 1 and the discussion provided in Remark 1. Further, Proposition 2 shows how the uncertain system (9) is certified stable with respect to a candidate parameter δ; one can certify that δ solves Problem 1 by computing a quadratic Lyapunov function for the i th -level system H i . Certifying a candidate parameter δ in this way may lead to conservative estimates of the stability marginδ. However, we show next how the homogeneous polynomial Lyapunov function, which is derived in the certification of δ, can be used to compute alternate parameters and, in this way, reduce conservatism.
Proposition 3. For a given i ≥ 1 and δ ≥ 0, assume that there exists a P ∈ S n i ++ satisfying the constraint (7) for S = {A, A + δA 0 }. Define by then δ P solves Problem 1 and δ P ≥ δ.
Proof. Choose integer i ≥ 1 and choose δ ≥ 0 satisfying the constraint (7) at the i th level for some P ∈ S n i ++ . Then (9) is stable with respect to δ as a result of Proposition 2. From this, we know that A is Hurwitz; this is due to the fact that when (9) is stable with respect to δ, (9) is stable with respect to 0. Moreover, the optimization problem (11) is always feasible as δ satisfies the constraint (12).
The constraint (12) is equivalent to the constraint (1), where here P ∈ S n i ++ is fixed, and δ * is a free variable. Therefore, if there exists a δ * ≥ δ satisfying (12) for P then the system (9) is safe with respect to δ * , i.e. δ * solves Problem 1, and δ * ≥ δ. This completes the proof.
We next apply the procedures detailed in Propositions 2 and 3 to form an iterative algorithm for under-approximatinĝ δ with reduced conservatism. This procedure is given in Algorithm 1.

Theorem 1.
If A is Hurwitz, then the output of Algorithm 1 solves Problem 1, for all ε ≥ 0 and all integers i max ≥ 1.

V. COMPUTING THE WORST-CASE SWITCHING FUNCTION
The Lyapunov function returned by Algorithm 1 proves asymptotic stability of (9) with ∆(t) ∈ [0, δ]. We should be able to identifyδ with a procedure such as Algorithm 1 since a homogeneous polynomial Lyapunov function is a necessary condition for the stabiity of (9). However, computing a Lyapunov function of sufficiently high order via semidefinite programming techniques may become computationally intractable for higher dimensional systems or as δ →δ. Therefore, we lay the groundwork for a numerical procedure to produce an upper bound forδ in this section.
The Lyapunov function returned by Algorithm 1 includes valuable information about how to shape ∆(t) in order to produce a periodic or diverging trajectory. Specifically, choosing a function ∆(t) ∈ [0, δ] which maximizes the timederivative of the Lyapunov function returned by Algorithm 1 will push the system as close as possible to instability for the given δ, especially as the order of the Lyapunov function is increased. This intuition is formalized in Proposition 4, and we now use it to approximate the worst-case switching function ∆ w (t; δ). δ ← δ + ε 5: for i ∈ {1, · · · , i max } do 6: Compute P i ∈ S n i ++ satisfying (7) for δ 8: if P i found then 9: δ ← δ Compute δ Pi satisfying (11) for P i

15:
if δ Pi found then 16: δ ← δ Pi 17: end function Proposition 4. Let (8) be a Lyapunov function, parameterized by P , for the system (9). Then the function ∆ w (t; δ, P ) which maximizesV (x(t)) is given by where the indicator function I(t) is given by for A = A i and A 0 = A i 0 is solved by (14). The optimization problem (16) shows that the ∆(t) selected is the value in the interval [0, δ] which maximizeṡ V (x(t)) for all t ≥ 0. The state-dependence of I(t) suggests that ∆ w (t; δ, P ) is not unique and generally differs for trajectories beginning at different initial conditions x 0 . Therefore, Algorithm 2 uses a simulation given an initial condition x 0 Algorithm 2 Solve Problem 2 input : A, A 0 from (9), desired margin δ. P from Algorithm 1. Simulation parameters: initial condition x 0 and time horizon t f . output: T, Σ from (17) to describe ∆ w (t; δ, P ).

Remark 3.
Since ∆ w (t; δ, P ) is a piecewise-constant function, it is most easily utilized in a numerical procedure when expressed as a pair of finite sets T = {t 0 , . . . , t f } and Σ = {σ 0 , . . . , σ f −1 } such that The worst-case switching function (14) contains all of the necessary information needed to compute an upper bound on the stability margin of (9).

VI. BOUNDING THE STABILITY MARGIN FROM ABOVE AND PRODUCING A PERIODIC TRAJECTORY
We now seek to approximate the lowest upper bound on the stability margin of (9) with an iterative procedure which, starting with δ = δ, increments δ until A w (t; δ, P ) can be shown to produce a periodic trajectory. Such a trajectory with an initial condition x 0 and time horizon t f can become periodic, but it is not necessarily periodic from x 0 . Therefore, the switching sequence described by T and Σ in (17) is searched to see if there exist indices k ≥ j ≥ 0 such that a Algorithm 3 Solve Problem 3 input : A, A 0 from (9), δ from Alg. 1.
P from Algorithm 1. Simulation parameters: initial condition x 0 and time horizon t f . output: δ, an upper bound onδ.
A d , discrete transition matrix certifying marginal stability.
Proposition 5. A switching function described by (17) with sets T and Σ produces a limit cycle for (9) if there exists some k ≥ j ≥ 0 such that the discrete transition matrix has an eigenvalue of magnitude equal to 1.
Algorithm 3 seeks to find the parameter δ which induces a limit cycle. In doing so, it calls Algorithm 2 to compute ∆ w (t; δ, P ) and searches the switching sequence for a discrete transition matrix A d that has an eigenvalue of magnitude 1. Such an A d certifies that ∆ w (t; δ, P ) produces a periodic trajectory.
Remark 4. The δ found by Algorithm 3 depends on the parameter P returned by Algorithm 1. A higher-order Lyapunov function will cause a smaller value of δ to be found. Therefore, δ is an upper bound onδ, and in order to find the lowest possible δ, a high-order Lyapunov function should be used with Algorithm 3, if possible.

VII. NUMERICAL EXAMPLES AND APPLICATIONS
The practical importance of the techniques presented in this article is highlighted in the following examples. Example 1. [Second Order System] We first study the system (9) with n = 2 and This system resembles a spring-mass damper with uncertainty in the spring constant. Algorithm 1 is used to compute a lower bound on the stability margin. The algorithm returns δ = 2.15 with a 14 th order polynomial Lyapunov function. By increasing i max to compute a 28 th order Lyapunov function, the algorithm computes δ = 2.16. Beyond that, only slight increases in δ can be obtained for significantly higherorder Lyapunov functions such that computation quickly becomes intractable. The next step is to find a δ for which we can produce a periodic trajectory.
Using the 14 th order Lyapunov function, Algorithm 3 provides an upper bound δ on the system's stability margin and provides the subset of the switching sequence from Algorithm 2 which certifies periodicity of the system (or near periodicity-the numerical method is approximate) with A w (t; δ, P ). The value δ = 2.21 was computed, and two switching sequences were returned: one using simulation data from 8] T , one using A w (t; δ, P 14 ) described by (T 1 , Σ 1 ) and the other using A w (t; δ, P 14 ) described by (T 2 , Σ 2 ). Both switching sequences produce a limit cycle but at different locations. This demonstrates that a function ∆ w (t; δ, P ) is not unique but that any such function can produce a periodic trajectory for the same value of δ.

Example 2. [Fourth Order
System] For higher-order systems, limit cycles cannot be visualized on the phase plane; instead we search for periodic orbits via analysis of the indicator function (15) and Algorithm 3. To demonstrate this, we study (9) with 4 th order dynamics inspired by the linearized, non-dimensional lateral dynamics of a fixed-wing aircraft with the following parameters [15]: Algorithm 1 computes a maximum value of δ = 0.24 for a 6th-order polynomial Lyapunov function. It may be possible to find a larger δ, but increasing i max to attain a higher-order Lyapunov function would make computation intractable without additionally employing an algorithm to reduce the redundancies produced by the use of the Kronecker product.
To get a better picture of system robustness, we run Algorithm 3 using the 6th-order Lyapunov function, and compute δ = 0.27. The switching sequence ∆(t; δ, P 6 ) which produces a periodic trajectory is returned, described by (T, Σ): Switching Sequence using x 0 = [1 1 1 1] T T 1 {0, 0.027, 0.060} Σ 1 {0.27, 0} Figure 2 plots the indicator function (15) for the cases where ∆ w (t; δ, P 6 ) and ∆ w (t; δ, P 6 ) are used. The former generates a limit cycle starting from about t = 0.332, while the latter produces values of I(t) which slowly decay over time. The initial erratic jumps in I(t) show that the trajectory using δ does not begin on a limit cycle but rather takes a few tenths of a second to reach the periodic orbit. Of practical interest to an analysis of uncertain systems are worst-case impulse responses. In [11], the authors use the system hierarchy (5) to easily compute exponential bounds on an impulse response. The methods in this paper can be used to find a worst-case impulse response which approaches these bounds. We add a control input and define an output for the system (9) and studẏ The impulse response h(t) of (21) is described parametrically by a solution ϕ(t) tȯ Since the impulse response of (21) is closely related to the unforced response, the algorithms presented above can be used to compute various worst-case and marginally stable trajectories in these scenarios as well.
Example 3 (Worst-Case Impulse Response). Consider the system studied in Example 1 with the additional parameters B = [0 1] T and C = [1 0]. We use Theorem 4 in [11], with parameter α = 0.11, to produce an exponential bound on the impulse response. Using Algorithm 2, a worst-case switching function ∆ w (t; δ = 1, P 14 ) is computed for (22). Figure  3 shows how this switching function produces an impulse response which approaches the bound much more closely than the system with no switching.

VIII. CONCLUSION
In this work, we show how the more general class of homogeneous polynomial Lyapunov functions is used to compute stability margins with reduced conservatism, and we show how these Lyapunov functions aid in the search for periodic trajectories for marginally stable LTV systems. Our work is premised on the recent observation that the search for a homogeneous polynomial Lyapunov function for some LTV systems is easily encoded as the search for a quadratic Lyapunov function for a related LTV system, and our main contribution is an intuitive algorithm for generating upper and lower bounds on the system's stability margin. We show also how the worst-case switching scheme-which draws an LTV system closest to a periodic orbit-is generated.