add isunstable (#1039)

* add `isunstable`

* bump version

* update docstring

Author: baggepinnen

lib/ControlSystemsBase/Project.toml

@@ -2,7 +2,7 @@ name = "ControlSystemsBase"
 uuid = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "1.20.1"
+version = "1.20.2"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -39,7 +39,7 @@ Hungarian = "0.7.0"
 ImplicitDifferentiation = "0.7.2"
 LinearAlgebra = "<0.0.1, 1"
 MacroTools = "0.5"
-Makie = "0.24"
+Makie = "0.22, 0.23, 0.24"
 MatrixEquations = "1, 2.1"
 MatrixPencils = "1.8.3"
 Polynomials = "3.0, 4.0"

lib/ControlSystemsBase/src/types/Lti.jl

@@ -107,10 +107,34 @@ common_timeevol(systems::LTISystem...) = common_timeevol(timeevol(sys) for sys i
 """
     isstable(sys)
 
-Returns `true` if `sys` is stable, else returns `false`."""
+Returns `true` if `sys` is stable, else returns `false`.
+
+Marginally stable systems are considered unstable by this function, see [`isunstable`](@ref) for a function that returns true only for exponentially unstable systems, that is, `!isunstable(sys)` implies that `sys` is either stable or marginally stable.
+"""
 isstable(sys::LTISystem{Continuous}) = all(real.(poles(sys)) .< 0)
 isstable(sys::LTISystem{<:Discrete}) = all(abs.(poles(sys)) .< 1)
 
+
+"""
+    isunstable(sys)
+
+Returns `true` if `sys` is exponentially unstable, else returns `false`.
+Marginally stable systems (systems with a simple poles on the imaginary axis) are considered stable by this function, see [`isstable`](@ref) for a function that returns true only for exponentially stable systems.
+"""
+function isunstable(sys::LTISystem)
+    inte, p, z, tolp, tolz = integrator_excess_with_tol(sys)
+    if inte > 1
+        return true
+    end
+    # Go through all poles on the imaginary axis and check if they are duplicated
+    for pi in p
+        abs(real(pi)) > sqrt(sqrt(eps(abs(pi)))) && continue # The pole is too far away to be on the imaginary axis
+        # check if there are multiple poles with this imaginary part
+        count_eigval_multiplicity(p, complex(0.0, imag(pi)))[1] > 1 && return true
+    end
+    return iscontinuous(sys) ? any(real.(p) .> tolp) : any(abs.(p) .> tolp)
+end
+
 # Fallback since LTISystem not AbstractArray
 Base.size(sys::LTISystem, i::Integer) = size(sys)[i]
 

lib/ControlSystemsBase/test/test_analysis.jl

@@ -428,4 +428,14 @@ zss  = tzeros(G)
 ptf  = poles(Gtf)
 ztf  = tzeros(Gtf)
 pzpk = poles(zpk(G))
-zzpk = tzeros(zpk(G))
+zzpk = tzeros(zpk(G))
+
+
+
+## Test isunstable
+using ControlSystemsBase: isunstable
+@test !isunstable(tf(1, [1,1]))
+@test isunstable(tf(1, [1,-1]))
+@test isunstable(tf(1, [1,0,0])) # Double integrator
+@test isunstable(zpk([1], [im, im, -im, -im], 1)) # Repeated pole on imaginary axis not in origin
+@test !isunstable(zpk([1], [im+1e-8im, im, -im-1e-8im, -im], 1)) # Almost repeated pole on imaginary axis not in origin