Validate integrator clusters using polynomial coefficients

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.25.0"
+version = "1.25.1"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

lib/ControlSystemsBase/src/analysis.jl

@@ -44,22 +44,40 @@ end
 
 
 """
-    n, tol = count_eigval_multiplicity(p, location, e = nothing; scale = maximum(abs, p))
+    n, tol = count_eigval_multiplicity(p, location, e = nothing; scale = maximum(abs, p), coefftol = nothing)
 
-Return the multiplicity `n` of the eigenvalue `location` among the computed eigenvalues `p`, together with the tolerance `tol` used to obtain it. The multiplicity is the first value of `i` for which exactly `i` elements of `p` lie within the radius `r(i)` around `location`. If no such `i` exists, `n = 0` and `tol = r(1)`.
+Return the multiplicity `n` of the eigenvalue `location` among the computed eigenvalues `p`, together with the tolerance `tol` used to obtain it. A candidate multiplicity `i` requires exactly `i` elements of `p` within the radius `r(i)` around `location`. Return the first candidate that also passes the coefficient check below. If none passes, `n = 0` and `tol = r(1)`.
 
-A backward-stable eigenvalue algorithm perturbs an eigenvalue of multiplicity `i` by an amount of order `ϵ^(1/i)` relative to the scale of the problem. The default radius is therefore `r(i) = (100ϵ)^(1/i) * scale`, which coincides with the tolerances for simple and double imaginary-axis eigenvalues in the Bruinsma–Steinbuch algorithm for the ``H_∞`` norm (`100ϵ` and `10√ϵ`). The radius is proportional to `scale`, so that the result does not depend on the time unit of the system. If `e` is provided, the radius `r(i) = e^(1/i)` is used instead.
+A backward-stable eigenvalue algorithm perturbs an eigenvalue of multiplicity `i` by an amount of order `ϵ^(1/i)` relative to the scale of the problem. The default radius is therefore `r(i) = (100ϵ)^(1/i) * scale`, which coincides with the tolerances for simple and double imaginary-axis eigenvalues in the Bruinsma–Steinbuch algorithm for the ``H_∞`` norm (`100ϵ` and `10√ϵ`). The radius is proportional to `scale`, so that the result does not depend on the time unit of the system.
+
+To avoid classifying unrelated slow modes as a repeated eigenvalue as the radius increases, form the monic polynomial with roots `(p[j] - location)/scale` for each candidate cluster. All coefficients other than the leading one must have magnitude at most `i*coefftol`, where `coefftol` defaults to `100ϵ`. This checks that the normalized cluster polynomial is close to `z^i`, and is invariant under a change of time unit for continuous-time integrators. This remains a heuristic based on computed eigenvalues, rather than a matrix backward-error bound.
+
+If `e` is provided, retain the explicit-radius rule `r(i) = e^(1/i)` without the coefficient check.
 """
-function count_eigval_multiplicity(p, location, e=nothing; scale=maximum(abs, p, init=0.0)) # The init is to handle poor type inference with exotic number types
+function count_eigval_multiplicity(p, location, e=nothing; scale=maximum(abs, p, init=0.0), coefftol=nothing) # The init is to handle poor type inference with exotic number types
     T = eltype(p)
     ϵ = T <: Number && isconcretetype(T) ? eps(real(float(T))) : eps(float(real(typeof(scale))))
+    η = coefftol === nothing ? 100ϵ : coefftol
     radius(i) = e === nothing ? (100ϵ)^(1/i)*scale : e^(1/i)
     n = length(p)
     n == 0 && return (0, zero(radius(1)))
     for i = 1:n
         # if we count i poles within the circle assuming i integrators, we return i
         tol = radius(i)
         if count(p->abs(p-location) <= tol, p) == i
+            if e === nothing && !iszero(scale)
+                q = [(pj-location)/scale for pj in p if abs(pj-location) <= tol]
+                # Ascending coefficients of prod(z - qj), built in place. Check all
+                # coefficients: a zero centroid alone also accepts symmetric slow modes.
+                coefficients = fill(zero(first(q)), i+1)
+                coefficients[end] = one(first(q))
+                for (j, qj) in enumerate(q)
+                    for k in i-j+1:i
+                        coefficients[k] -= qj*coefficients[k+1]
+                    end
+                end
+                maximum(abs, @view(coefficients[1:i])) <= i*η || continue
+            end
             return (i, tol)
         end
     end
@@ -69,7 +87,7 @@ end
 """
     count_integrators(P)
 
-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 `(100ϵ)^(1/n) * maximum(abs, p)` around the origin (1 in discrete time) equals `n`, see `ControlSystemsBase.count_eigval_multiplicity`.
+Count the number of poles in the origin (1 in discrete time), allowing for numerical splitting of repeated poles. Candidate clusters within a radius `(100ϵ)^(1/n) * maximum(abs, p)` must also have a normalized polynomial close to `z^n`, see `ControlSystemsBase.count_eigval_multiplicity`.
 
 See also [`integrator_excess`](@ref).
 """

lib/ControlSystemsBase/test/test_analysis.jl

@@ -652,6 +652,56 @@ P = let
 end
 @test ControlSystemsBase.count_integrators(P) == 2
 
+
+@testset "Integrator cluster coefficients" begin
+    count_mult = ControlSystemsBase.count_eigval_multiplicity
+    # The radius-only rule classified the slow modes as integrators when the
+    # candidate radius included every pole except the one with largest magnitude.
+    A = [
+        -0.02314734339987039 -0.23470977613890898 0.40744939133965374 -0.614172004152354 -0.7336985077201985 -0.7188547646417938 -0.6116827006955479 0.28701228553159197 28.932835460401698 -0.8577016118934753 -0.002900248526486404
+        0.04747182843197798 0.028473586849358613 -0.7336202770573863 0.8169604093939422 1.0737298575562346 0.22485825328995757 0.0779364215281102 -0.010519719302367558 0.8399890329652873 -0.026616152383339912 -8.029617097986969e-5
+        0.0 0.5169499360495327 -0.18168990332330218 0.6841609486445348 0.7588442551565007 0.15992175120178942 0.024995361037593256 0.004936045695480118 1.8302561456391446 -0.054932217097018884 -0.00016795573573835654
+        0.0 0.0 -0.013220361979022434 0.01230952085420172 -0.9010774527457568 0.02301808232117536 -0.0404752161566272 -0.03731112146998583 -2.7446448440610656 0.09802006669271178 0.0005426830092505972
+        0.0 0.0 0.0 0.6076623552098359 -0.40125742352027083 0.05519756476123633 -0.10706393385556773 -0.036672569915258635 0.9288661567402832 -0.0337764354485775 -0.0002001539651115842
+        0.0 0.0 0.0 0.0 -0.14787625350432967 -0.06636050823375049 -0.33810141209957845 0.3256454596279698 8.326333586378503 -0.29552019036140387 -0.0015981690223459586
+        0.0 0.0 0.0 0.0 0.0 0.4530941732907558 0.10965212099932263 0.20248460567824594 -10.804442594382103 0.3896892668191028 0.0022606629283560917
+        0.0 0.0 0.0 0.0 0.0 0.0 0.07781306774762253 -0.11421090913742127 -9.709505482383845 0.3621161336134016 0.0025039855060416507
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.13310699861538938 -19.935290888483895 -0.002413679382998811 -1.670177472499906e-5
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.48501499005802 -0.10382400507966154 -0.001266235535267126
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.9999995379331097 -2.2768266629953853e-8
+    ]
+    for c in (1e-4, 1.0, 1e4)
+        P = ss(c*A, ones(11,1), ones(1,11), 0)
+        @test ControlSystemsBase.count_integrators(P) == 0
+        n, tol = count_mult(poles(P), 0)
+        @test n == 0
+        @test tol ≈ 100eps()*maximum(abs, poles(P))
+    end
+
+    @test count_mult(Float64[], 0) == (0, 0.0)
+    @test count_mult(zeros(15), 0) == (15, 0.0)
+    @test count_mult(ones(5), 1)[1] == 5
+    @test count_mult(Number[0.0, -1.0], 0)[1] == 1
+    for T in (Float32, Float64, BigFloat), m in (3, 5, 10, 15)
+        # Splitting of an m-fold root, with a nonzero pole setting the scale.
+        r = eps(T)^(one(T)/m)
+        cluster = r .* cispi.(2 .* T.(0:m-1) ./ m)
+        @test count_mult([cluster; -one(T)], 0)[1] == m
+        @test count_mult([cluster .+ T(0.1); -one(T)], 0)[1] == 0
+    end
+    # The centroid vanishes, but the quadratic coefficient reveals distinct roots.
+    for pair in ([0.01, -0.01], [0.01im, -0.01im])
+        p = [repeat(pair, 5); -20]
+        @test count_mult(p, 0)[1] == 0
+        @test count_mult(p, 0; coefftol=1e-6)[1] == 10
+        # Zeros use the same check when computing integrator excess.
+        G = zpk(p, [-1.0], 1)
+        @test ControlSystemsBase.integrator_excess(G) == 0
+    end
+    # Preserve the existing explicit-radius override.
+    @test count_mult([-0.05, 0.05], 0, 0.01) == (2, 0.1)
+end
+
 ## Difficult test case for zeros
 G = let
      tempA = [-0.6841991610512457 -0.0840213470263692 -0.0004269818661494616 -2.7625001165862086e-18; 2.081491248616774 0.0 0.0 8.404160870560225e-18; 0.0 24.837541148074962 0.12622006230897712 0.0; -1.2211265763794115e-14 -2.778983834717109e8 -1.4122312296634873e6 -4.930380657631326e-32]
@@ -701,4 +751,4 @@ using ControlSystemsBase: isunstable
 @test isunstable(tf(1, [1,2,1], 1)) # Double pole at z = -1
 @test !isunstable(ss(c2d(tf(1, [1,0,1]), 0.1))) # Simple pole pair on the unit circle
 @test isunstable(ss(c2d(tf(1, [1,0,1])^2, 0.1))) # Repeated pole pair on the unit circle
-@test isunstable(zpk(Float64[], [cis(0.5), cis(0.5), cis(-0.5), cis(-0.5)], 1.0, 1))
+@test isunstable(zpk(Float64[], [cis(0.5), cis(0.5), cis(-0.5), cis(-0.5)], 1.0, 1))