Make the tolerance of the integrator count independent of the time unit

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,20 +87,16 @@ 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)
     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)
+    # 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
 

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
+@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