Joint Spectral Radius (JSR)

The Joint Spectral Radius (JSR) of a discrete-time switched system is the minimal value of γ such that the system scaled by γ is asymptotically stable.

SwitchOnSafety.ScaledHybridSystem — Type

struct ScaledHybridSystem{T, H <: Union{DiscreteSwitchedLinearSystem, ConstrainedDiscreteSwitchedLinearSystem}} <: HybridSystems.AbstractHybridSystem
    system::H
    γ::T
end

Discrete-time system where each reset map is scaled by γ, that is, the reset map x ↦ Ax of system is replaced by x ↦ Ax/γ.

source

The JSR is NP-hard to compute but several methods exist to approximate this quantity.

Gripenberg algorithm, a Branch-and-Bound approach

The following approachs searches through all the different switching sequences, pruning some using the running estimate of the lower bound.

SwitchOnSafety.gripenberg — Function

gripenberg(s::AbstractDiscreteSwitchedSystem; δ=1e-2,
                max_eval = 10000, max_ρ_eval = max_eval,
                max_norm_eval = max_eval, max_length = 50,
                matrix_norm = A -> opnorm(A, 2), verbose = 1)

Gripenberg algorithm [G96] for computing an upper bound ub and a lower bound lb to the Joint Spectral Radius such that ub - lb ≤ 1e-2.

[G96] Gripenberg, G. Computing the joint spectral radius. Linear Algebra and its Applications, Elsevier, 1996, 234, 43-60

source

Sum-of-Squares approach

The following method computes upper bounds using Sum Of Squares Programming

SwitchOnSafety.soslyap — Function

soslyap(s::AbstractSwitchedSystem, d; optimizer_constructor=nothing)

Find Sum-of-Squares Lyapunov functions; i.e. solves [(5), PJ08] or gives moment matrices certifying the infeasibility of the problem. Use ScaledHybridSystem to use a different growth rate than 1.

[PJ08] P. Parrilo and A. Jadbabaie. Approximation of the joint spectral radius using sum of squares. Linear Algebra and its Applications, Elsevier, 2008, 428, 2385-2402

source

SwitchOnSafety.soslyapbs — Function

soslyapbs(s::AbstractSwitchedSystem, d::Integer,
          soslb, dual,
          sosub, primal;
          verbose=0, tol=1e-5, step=0.5, scaling=quickub(s),
          ranktols=tol, disttols=tol, kws...)

Find the smallest γ such that soslyap is feasible.

source

The infeasibility certificates computed in the binary search carried out by soslyapbs can be used to produce cycles of high growth rate using findsmp on the sequence produced by sosbuildsequence.

SwitchOnSafety.sosbuildsequence — Function

sosbuildsequence(s::AbstractSwitchedSystem, d::Integer;
                 v_0=:Random, p_0=:Random, l::Integer=1,
                 Δt::Float64=1., niter::Integer=42,
                 kws...)

Compute the truncation of length l of the high growth rate infinite sequence produced by the algorithm introduced in [LJP17]. The trajectory ends at mode v_0 (or a random one if v_0 is :Random) and is built backward as explained in [LJP17]. The measures used to guide the construction are the infeasibility certificates of highest growth rate computed by soslyap with polynomials of degree 2d. The starting polynomial is either p_0, or a random strictly sum-of-squares polynomial if p_0 is :Random or the primal solution of soslyap certifying the best upper bound for mode v_0.

[LJP17] B. Legat, R. M. Jungers, and P. A. Parrilo.

Certifying unstability of Switched Systems using Sum of Squares Programming, arXiv preprint arXiv:1710.01814, 2017.

source

SwitchOnSafety.findsmp — Function

findsmp(seq::HybridSystems.DiscreteSwitchingSequence)

Extract the cycle of highest growth rate in the sequence seq.

source

Alternatively, sosextractcycle can be used to find cycles of high growth rate.

SwitchOnSafety.sosextractcycle — Function

sosextractcycle(s::AbstractDiscreteSwitchedSystem, dual, d::Integer;
                ranktols=1e-5, disttols=1e-5)

Extract cycles of high growth rate from atomic occupation measures given by the infeasibility certificates of highest growth rate computed by soslyap. The method is detailed in [LJP17].

[LJP17] B. Legat, R. M. Jungers, and P. A. Parrilo.

Certifying unstability of Switched Systems using Sum of Squares Programming, arXiv preprint arXiv:1710.01814, 2017.

source

Polytopic approach

The following method can verify numerically that a cycle is an s.m.p. by computing invariant polytopes. When it is not an s.m.p., it find a candidate of higher growth rate.

SwitchOnSafety.invariant_polytopes — Function

Computes invariant balanced polytopes for system s using the algorithm of [GP13] with the spectral maximizing product candidate (s.m.p.) smp. The choice between Complex Balanced Polytope, Real Balanced Polytope or Conitope using invariance of the positive orthant is done automatically. The solver used is optimizer_constructor. This solver should support second order cone if the Complex Balanced Polytope is used (e.g. if the leading eigenvalue has an imaginary part). The maximum depth considered is max_length, a point x is considered to be in the interior of the polytope if the maximum λ such that λ * x is in the polytope is larger than 1 + tol.

New s.m.p. candidates are considered with tolerance new_candidate_tol. A maximum of max_cycles times the current s.m.p. candidate is prepended to this new candidate to try to find a candidate with growth rate higher than the current one, we also stop adding cycles if the length reaches max_smp_length. If gready, then cycles continue to be prepended even if the it has a higher growth rate than the current s.m.p. until the growth rate decreases when prepending a new cycle.

If verbose is 1, the number of leaves is printed when the depth is a multiple of log_step_length.

[GP13] N. Guglielmi and V. Protasov. Exact computation of joint spectral characteristics of linear operators. Foundations of Computational Mathematics 13.1, 2013, 37-97.

source