From 427b379904d055865b61103fed8fe34ff6f21388 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 21 Oct 2026 18:22:38 +0000 Subject: [PATCH 1/3] Validate integrator clusters using polynomial coefficients --- lib/ControlSystemsBase/src/analysis.jl | 28 +++++++++-- lib/ControlSystemsBase/test/test_analysis.jl | 50 ++++++++++++++++++++ 2 files changed, 73 insertions(+), 5 deletions(-) diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index f400d38e7..26200e210 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/lib/ControlSystemsBase/src/analysis.jl @@ -44,15 +44,20 @@ 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 mistaking unrelated slow modes for a repeated eigenvalue as the radius grows, 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))) @@ -60,6 +65,19 @@ function count_eigval_multiplicity(p, location, e=nothing; scale=maximum(abs, p, # 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). """ diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index fc0814f46..f02330ea3 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/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 + # An 11-state system with slow modes and a pole at -20 was reported as + # having ten integrators by the radius-only rule (r(10) ≈ 0.862). + 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] From 6d828a3f98e36854e9c1cb5b5b06c19d1bab58ed Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 21 Oct 2026 18:24:08 +0000 Subject: [PATCH 2/3] Refine terminology in integrator validation documentation --- lib/ControlSystemsBase/src/analysis.jl | 2 +- lib/ControlSystemsBase/test/test_analysis.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index 26200e210..3d3616f65 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/lib/ControlSystemsBase/src/analysis.jl @@ -50,7 +50,7 @@ Return the multiplicity `n` of the eigenvalue `location` among the computed eige 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 mistaking unrelated slow modes for a repeated eigenvalue as the radius grows, 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. +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. """ diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index f02330ea3..a4c76751e 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/lib/ControlSystemsBase/test/test_analysis.jl @@ -655,8 +655,8 @@ end @testset "Integrator cluster coefficients" begin count_mult = ControlSystemsBase.count_eigval_multiplicity - # An 11-state system with slow modes and a pole at -20 was reported as - # having ten integrators by the radius-only rule (r(10) ≈ 0.862). + # 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 @@ -751,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)) \ No newline at end of file +@test isunstable(zpk(Float64[], [cis(0.5), cis(0.5), cis(-0.5), cis(-0.5)], 1.0, 1)) From 609d3cb74ae5f73a6f7cfb77abcb8f2e02718a65 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 21 Oct 2026 19:01:42 +0000 Subject: [PATCH 3/3] Bump ControlSystemsBase to 1.25.1 --- lib/ControlSystemsBase/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml index 8ee41be07..2bee2477f 100644 --- a/lib/ControlSystemsBase/Project.toml +++ b/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"