tweak computation of integrator excess

lib/ControlSystemsBase/Project.toml

@@ -48,9 +48,10 @@ julia = "1.6"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -59,5 +60,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "GR", "Plots", "SparseArrays", "StaticArrays"]
-
+test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "GenericLinearAlgebra", "GR", "Plots", "SparseArrays", "StaticArrays"]

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -242,6 +242,8 @@
             print(io, "\n$(exc.f) with continuous-time systems, including delay systems and nonlinear systems, require the user to first ")
             printstyled(io, "install and load ControlSystems.jl, or pass the keyword method = :zoh", color=:green, bold=true)
             print(io, " for automatic discretization (applicable to systems without delays or nonlinearities only).")
+        elseif exc.f ∈ (eigvals!, ) && argtypes[1] <: AbstractMatrix{<:Number}
+            printstyled(io, "\nComputing eigenvalues of a matrix with exotic element types may require `using GenericLinearAlgebra`.", color=:green, bold=true)
         end
         plots_id = Base.PkgId(UUID("91a5bcdd-55d7-5caf-9e0b-520d859cae80"), "Plots")
         if exc.f isa Function && nameof(exc.f) === :plot && parentmodule(argtypes[1]) == @__MODULE__() && !haskey(Base.loaded_modules, plots_id)

lib/ControlSystemsBase/src/analysis.jl

@@ -38,24 +38,22 @@
 end
 
 
-function count_eigval_multiplicity(p, location)
-    T = float(real(eltype(p)))
+function count_eigval_multiplicity(p, location, e=eps(maximum(abs, p)))
     n = length(p)
     n == 0 && return 0
-    maximum(1:n) do i
+    for i = 1:n
         # if we count i poles within the circle assuming i integrators, we return i
-        if count(p->abs(p-location) < (i+1)*(eps(T)^(1/i)), p) >= i
-            i
-        else
-            0
+        if count(p->abs(p-location) < (e^(1/i)), p) == i
+            return i
         end
     end
+    0
 end
 
 """
     count_integrators(P)
 
-Count the number of poles in the origin by finding the maximum value of `n` for which the number of poles within a circle of radius `(n+1)*eps(numeric_type(sys))^(1/n)` arount the origin (1 in discrete time) equals `n`.
+Count the number of poles in the origin by finding the first value of `n` for which the number of poles within a circle of radius `eps(maximum(abs, p))^(1/n)` around the origin (1 in discrete time) equals `n`.
 
 See also [`integrator_excess`](@ref).
 """

lib/ControlSystemsBase/test/test_analysis.jl

@@ -1,3 +1,5 @@
+@test_throws MethodError poles(big(1.0)*ssrand(1,1,1)) # This errors before loading GenericLinearAlgebra
+using GenericLinearAlgebra # Required to compute eigvals of a matrix with exotic element types
 @testset "test_analysis" begin
 ## tzeros ##
 # Examples from the Emami-Naeini & Van Dooren Paper
@@ -269,6 +271,15 @@ for wgm in wgm[]
      @test mod(rad2deg(angle(freqresp(OL, wgm)[])), 360)-180 ≈ 0 atol=1e-1
 end
 
+nint = ControlSystemsBase.count_integrators(pade(OL, 2))
+nintbig = ControlSystemsBase.count_integrators(big(1.0)*pade(OL, 2))
+@test nint == nintbig # This one is very tricky and tests the threshold value of the eigval counting
+
+@test ControlSystemsBase.count_integrators(pade(OL, 3)) == nintbig
+@test ControlSystemsBase.count_integrators(pade(OL, 4)) == nintbig
+@test ControlSystemsBase.count_integrators(pade(OL, 10)) == nintbig
+
+
 # RGA
 a = 10
 P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)