From 3b92a274d220e4f2a26c86d213bf45de7b1d869d Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 7 Oct 2026 14:35:27 +0000 Subject: [PATCH] Make the tolerance of count_eigval_multiplicity independent of the time unit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The default radius for multiplicity i was eps(ρ)^(1/i), which is proportional to ρ^(1/i) for i ≥ 2, so the number of integrators depended on the time unit of the model. A backward-stable eigenvalue algorithm perturbs an eigenvalue of multiplicity i by an amount of order ϵ^(1/i) times the scale of the matrix. The default radius is now (100ϵ)^(1/i) * maximum(abs, p), which coincides with the tolerances for simple and double imaginary-axis eigenvalues in the Bruinsma–Steinbuch algorithm. When no multiplicity is found, the returned tolerance is the radius for multiplicity one, since isunstable and sisomargin use it as the threshold on the real part of unstable poles. In integrator_excess, the zeros are counted with the scale of the poles and the zeros. Co-Authored-By: Claude Opus 5.5 (1M context) --- lib/ControlSystemsBase/src/analysis.jl | 33 +++++++++++--------- lib/ControlSystemsBase/test/test_analysis.jl | 33 +++++++++++++++++++- 2 files changed, 51 insertions(+), 15 deletions(-) diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index d37fa60a2..4f99e5206 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/lib/ControlSystemsBase/src/analysis.jl @@ -43,24 +43,33 @@ function poles(sys::TransferFunction{<:TimeEvolution,SisoZpk{T,TR}}) where {T, T end -function count_eigval_multiplicity(p, location, e=eps(maximum(abs, p, init=0.0))) # The init is to handle poor type inference with exotic number types +""" + n, tol = count_eigval_multiplicity(p, location, e = nothing; scale = maximum(abs, p)) + +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)`. + +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. +""" +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 + T = eltype(p) + ϵ = T <: Number && isconcretetype(T) ? eps(real(float(T))) : eps(float(real(typeof(scale)))) + radius(i) = e === nothing ? (100ϵ)^(1/i)*scale : e^(1/i) n = length(p) - tol = zero(e) - n == 0 && return (0, tol) + 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 = e^(1/i) - if count(p->abs(p-location) < tol, p) == i + tol = radius(i) + if count(p->abs(p-location) <= tol, p) == i return (i, tol) end end - (0, tol) + (0, radius(1)) 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 `eps(maximum(abs, p))^(1/n)` around 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 `(100ϵ)^(1/n) * maximum(abs, p)` around the origin (1 in discrete time) equals `n`, see `ControlSystemsBase.count_eigval_multiplicity`. See also [`integrator_excess`](@ref). """ @@ -78,12 +87,7 @@ Count the number of integrators in the system by finding the difference between For discrete-time systems, the origin ``s = 0`` is replaced by the point ``z = 1``. """ function integrator_excess(P::LTISystem) - p = poles(P) - z = tzeros(P) - location = iscontinuous(P) ? 0 : 1 - np, tolp = count_eigval_multiplicity(p, location) - nz, tolz = count_eigval_multiplicity(z, location) - np - nz + integrator_excess_with_tol(P)[1] end function integrator_excess_with_tol(P::LTISystem) @@ -91,7 +95,8 @@ function integrator_excess_with_tol(P::LTISystem) z = tzeros(P) location = iscontinuous(P) ? 0 : 1 np, tolp = count_eigval_multiplicity(p, location) - nz, tolz = count_eigval_multiplicity(z, location) + # The zeros are computed from the system pencil, whose scale is at least that of the poles. A single zero in the origin does not provide a scale of its own. + nz, tolz = count_eigval_multiplicity(z, location; scale=max(maximum(abs, p, init=0.0), maximum(abs, z, init=0.0))) np - nz, p, z, tolp, tolz end diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index feb03d9e8..038a9c34a 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/lib/ControlSystemsBase/test/test_analysis.jl @@ -304,6 +304,36 @@ nintbig = ControlSystemsBase.count_integrators(big(1.0)*pade(OL, 2)) @test ControlSystemsBase.count_integrators(pade(OL, 4)) == nintbig @test ControlSystemsBase.count_integrators(pade(OL, 10)) == nintbig +# The integrator count must not depend on the time unit of the system +let + # The orthogonal similarity transformation prevents the eigenvalues in the origin from being computed exactly + Q = Matrix(qr([1.0 2 3 4; 2 -1 0 1; 3 0 2 -2; 4 1 -2 1]).Q) + G = ss(tf(1, [1, 3, 2, 0, 0])) + Gq = ss(Q'*G.A*Q, Q'*G.B, G.C*Q, G.D) + for c in (1e-4, 1.0, 1e4) + @test ControlSystemsBase.count_integrators(ss(c*Gq.A, c*Gq.B, Gq.C, Gq.D)) == 2 + end + @test ControlSystemsBase.count_integrators(ss(Float32.(Gq.A), Float32.(Gq.B), Float32.(Gq.C), Float32.(Gq.D))) == 2 + + # Slow dynamics without integrators + P = zpk(Float64[], -collect(1.0:10), 2*prod(1.0:10)) + Pslow = zpk(Float64[], -1e-3 .* collect(1.0:10), 2*prod(1e-3 .* (1.0:10))) + @test ControlSystemsBase.count_integrators(Pslow) == 0 + @test margin(Pslow).pm[] ≈ margin(P).pm[] + + Q3 = Matrix(qr([1.0 2 3; 2 -1 0; 3 0 2]).Q) + Gd = c2d(ss(tf(1, [1, 3, 2, 0])), 0.1) + @test ControlSystemsBase.count_integrators(ss(Q3'*Gd.A*Q3, Q3'*Gd.B, Gd.C*Q3, Gd.D, 0.1)) == 1 +end +@test ControlSystemsBase.count_integrators(DemoSystems.double_mass_model()) == 1 # The pole in the origin is not computed exactly +let + # A single zero in the origin that is not computed exactly + θ = 0.3 + Q = [cos(θ) sin(θ); -sin(θ) cos(θ)] + G = ss(tf([1, 0], [1, 2, 1])) + @test ControlSystemsBase.integrator_excess(ss(Q'*G.A*Q, Q'*G.B, G.C*Q, G.D)) == -1 +end + ## G = let # From https://juliacomputing.github.io/JuliaSimControl.jl/v0.10/examples/mtk_disturbance_modeling/#Controller-reduction marginplottestA = [-0.8778418152658126 -0.5578161193611026 0.6099186362718944 -0.739076573488832 0.0 -0.7985033210421671 -0.6437387074246465 -1.09990118296936 -0.9415890954623862 -0.6105710677920562 -0.40609522036480955 -0.02167670496721598 -0.0922028349547116 -0.3659595949709856 0.0 0.0 0.0 0.0; -0.5578161193611026 -0.8220564196754826 -0.4879223680823921 0.8803386749291978 0.0 -0.5305709989906725 -0.42546859154084343 -0.7310995753661961 -0.6252826290802145 -0.4179592857573079 -0.274638775913742 -0.00903355930082124 -0.05630074104219519 -0.24851140114082376 0.0 0.0 0.0 0.0; -12.772894974281625 7.422755287471556 -8.893245512259899 -0.5983738877209197 4.4676642940721933e-7 -2.901639071283097 -2.3556062284249624 -3.998192111860942 -3.428202120843976 -2.1327661674537324 -1.4414918123884637 -0.11684112979233667 -0.37069980312317624 -1.2936533657321307 0.0 0.0 0.0 0.0; 9.260923426511168 -10.119661325070803 1.5927557243537434 -6.10662600432525 0.0 -2.4254235921784115 -1.9671309526934084 -3.3286180099265987 -2.8499659817055276 -1.7302731969932739 -1.1837682861777055 -0.11782983264763462 -0.3265440770990644 -1.0546003416785925 0.0 0.0 0.0 0.0; -1.6061109633450728 0.3963654501432292 -8.331365272058026 -5.264727157230805 -0.001 -7.627995401504769 -6.199926026836587 -10.51167941872737 -9.015447054009817 -5.5777954205960585 -3.778335116508615 -0.3202631508161873 -0.9869003787621273 -3.388660360774488 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -3.548140511799198 -2.6482424961850004 -7.102082038003319 -6.208194955680467 0.0 0.0 0.0 0.0 0.0 -3.548140511799198 -2.6482424961850004 -8.102082038003319 -6.208194955680467; 0.0 0.0 0.0 0.0 0.0 -2.648242496185001 -2.0513640329992824 -5.597632716308486 -3.3299392371108265 0.0 0.0 0.0 0.0 0.0 -2.648242496185001 -2.0513640329992824 -5.597632716308486 -4.3299392371108265; 0.0 0.0 0.0 0.0 0.0 -18.102082038003317 4.4023672836915155 -24.873852981120933 -13.163260306077984 2.3828136105535203 2.0893223444460514 2.283219179019077 2.191129612074663 0.9999995532335706 -8.102082038003319 -5.5976327163084845 -21.873852981120933 -16.163260306077984; 0.0 0.0 0.0 0.0 0.0 3.7918050443195312 -14.329939237110827 -13.16326030607799 -15.408656206873678 0.0 0.0 0.0 0.0 0.0 -6.208194955680469 -4.3299392371108265 -16.16326030607799 -12.408656206873678; 0.0 0.0 0.0 0.0 0.0 -2.350321265186566 -1.6840654718182213 -5.594310418261385 -4.101569094377303 -1.488412883057869 -0.9639113397259123 0.5882419313046784 -0.8312794084435436 -0.3659595949709856 -1.5518179441443989 -1.040326764393575 -4.494409235292025 -3.159979998914917; 0.0 0.0 0.0 0.0 0.0 -2.2279953968524775 -1.5507658631630432 -5.630031465101207 -4.120845772910192 -0.9757754051184105 -1.0966951955892246 -0.4969559273832133 0.8240379338870026 -0.24851140114082376 -1.6974243978618049 -1.1252972716221998 -4.898931889735011 -3.4955631438299775; 0.0 0.0 0.0 0.0 0.0 -3.8534645754285197 -3.020606376036092 -6.73680504859119 -5.278229500140075 -14.905661141735358 5.981263475083093 -9.010086642052235 -0.969073690844096 -1.2936529189657013 -0.9518255041454229 -0.6650001476111296 -2.7386129367302487 -1.850027379296099; 0.0 0.0 0.0 0.0 0.0 -4.441954983917959 -3.3488668822504266 -9.320785111635601 -6.79570619414841 7.530650229517894 -11.303429611248509 1.4749258917061088 -6.433170081424314 -1.0546003416785925 -2.016531391739547 -1.381735929557018 -5.9921671017090015 -3.945740212442883; 0.0 0.0 0.0 0.0 0.0 -24.18765275168101 -17.451102917064617 -57.03910030189987 -42.6620899621136 -7.183906383941132 -3.3819696663653858 -8.651628422874213 -6.251627535992933 -3.3896603607744877 -16.55965735017624 -11.251176890228031 -46.5274208831725 -33.64664290810378; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.0 10.0 -3.0 3.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 10.0 -10.0 3.0 -3.0] @@ -455,4 +485,5 @@ using ControlSystemsBase: isunstable @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 \ No newline at end of file +@test !isunstable(zpk([1], [im+1e-8im, im, -im-1e-8im, -im], 1)) # Almost repeated pole on imaginary axis not in origin +@test isunstable(zpk(Float64[], [0.01; -collect(1.0:9)], 1.0)) # Slightly unstable pole in a system with a large state dimension \ No newline at end of file