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
— Typestruct 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/γ
.
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
— Functiongripenberg(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
Sum-of-Squares approach
The following method computes upper bounds using Sum Of Squares Programming
SwitchOnSafety.soslyap
— Functionsoslyap(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
SwitchOnSafety.soslyapbs
— Functionsoslyapbs(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.
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
— Functionsosbuildsequence(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.
SwitchOnSafety.findsmp
— Functionfindsmp(seq::HybridSystems.DiscreteSwitchingSequence)
Extract the cycle of highest growth rate in the sequence seq
.
Alternatively, sosextractcycle
can be used to find cycles of high growth rate.
SwitchOnSafety.sosextractcycle
— Functionsosextractcycle(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.
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
— FunctionComputes 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.