tweak computation of integrator excess

to handle new difficult test case

Author: baggepinnen

lib/ControlSystemsBase/src/analysis.jl

@@ -38,24 +38,22 @@ function poles(sys::TransferFunction{<:TimeEvolution,SisoZpk{T,TR}}) where {T, T
 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 maximum 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

@@ -250,6 +250,33 @@ P = tf(1,[5, 10.25, 6.25, 1])
 w_180, gm, w_c, pm = margin(50P)
 @test pm[] ≈ -35.1 rtol=1e-2
 
+## Tricky case from https://www.reddit.com/r/ControlTheory/comments/1inhxsz/understanding_stability_in_highorder/
+s = tf("s")
+kpu = -10.593216768722073; kiu = -0.00063; t = 1000; tau = 180; a = 1/8.3738067325406132e-5;
+kpd = 15.92190277847431; kid = 0.000790960718241793;
+kpo = -10.39321676872207317; kio = -0.00063;
+kpb = kpd; kib = kid;
+
+C1 = (kpu + kiu/s)*(1/(t*s + 1))
+C2 = (kpu + kiu/s)*(1/(t*s + 1))
+C3 = (kpo + kio/s)*(1/(t*s + 1))
+Cb = (kpb + kib/s)*(1/(t*s + 1))
+OL = (ss(Cb)*ss(C1)*ss(C2)*ss(C3)*exp(-3*tau*s))/((C1 - a*s)*(C2 - a*s)*(C3 - a*s));
+
+wgm, gm, ωϕₘ, ϕₘ = margin(OL; full=true, allMargins=true)
+@test ϕₘ[][] ≈ -320 rtol=1e-2
+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)