Use algorithms for poles and zeros of rational matrices from MatrixPencils.jl

lib/ControlSystemsBase/src/analysis.jl

@@ -11,7 +11,10 @@ poles(sys::AbstractStateSpace) = eigvalsnosort(sys.A)
 
 # Seems to have a lot of rounding problems if we run the full thing with sisorational,
 # converting to zpk before works better in the cases I have tested.
-poles(sys::TransferFunction) = poles(zpk(sys))
+function poles(sys::TransferFunction{<:TimeEvolution,SisoRational{T}}; kwargs...) where {T}
+    n,d = numpoly(sys), denpoly(sys)
+    sort!(filter!(isfinite, MatrixPencils.rmpoles(n, d; kwargs...)[1]), by=abs)
+end
 
 function poles(sys::TransferFunction{<:TimeEvolution,SisoZpk{T,TR}}) where {T, TR}
     # With right TR, this code works for any SisoTf
@@ -248,11 +251,20 @@ If `sys` is a state-space system the function has additional keyword arguments,
 
 To compute zeros of a system with non-BLAS floats, such as `BigFloat`, install and load the package `GenericSchur.jl` before calling `tzeros`.
 """
-function tzeros(sys::TransferFunction)
+function tzeros(sys::TransferFunction{<:TimeEvolution,SisoRational{T}}; kwargs...) where {T}
+    if issiso(sys)
+        return tzeros(sys.matrix[1,1])
+    else
+        n,d = numpoly(sys), denpoly(sys)
+        filter!(isfinite, MatrixPencils.rmzeros(n, d; kwargs...)[1]) # This uses rm2ls -> rm2lspm internally
+    end
+end
+
+function tzeros(sys::TransferFunction{<:TimeEvolution,SisoZpk{T,TR}}) where {T, TR}
     if issiso(sys)
         return tzeros(sys.matrix[1,1])
     else
-        return tzeros(ss(sys, minimal=true))
+        return tzeros(ss(sys))  # Convert to ss since rmzeros does this anyways, so no reason to pass through tf{SisoRational}
     end
 end
 

lib/ControlSystemsBase/test/test_analysis.jl

@@ -112,8 +112,10 @@ sys = s*(s + 1)*(s^2 + 1)*(s - 3)/((s + 1)*(s + 4)*(s - 4))
 @test tzeros(sys) ≈ [3.0, -1.0, im, -im, 0.0]
 
 ## POLE ##
-@test poles(sys) ≈ [4.0, -4.0, -1.0]
-@test poles([sys sys]) ≈ [4.0, -4.0, -1.0] # Issue #81
+@test es(poles(sys)) ≈ es([4.0, -4.0, -1.0])
+@test es(poles(sys, atol=1e-6)) ≈ es([4.0, -4.0]) # This cancels the -1 pole/zero
+@test es(poles([sys sys], atol=1e-6)) ≈ es([4.0, -4.0]) # Issue #81
+@test poles(zpk([sys sys])) ≈ [4.0, -4.0, -1.0] # This does not cancel the -1 pole/zero
 @test poles(ex_11) ≈ eigvals(ex_11.A)
 @test poles([2/(s+1) 3/(s+2); 1/(s+1) 1/(s+1)]) ≈ [-1, -1, -2]
 
@@ -198,7 +200,7 @@ sys = s*(s + 1)*(s^2 + 1)*(s - 3)/((s + 1)*(s + 4)*(s - 4))
 
 # Example 5.5 from https://www.control.lth.se/fileadmin/control/Education/EngineeringProgram/FRTN10/2019/e05_both.pdf
 G = [1/(s+2) -1/(s+2); 1/(s+2) (s+1)/(s+2)]
-@test_broken length(poles(G)) == 1 # The tf poles don't understand the cancellations
+@test length(poles(G)) == 1 # The tf poles do understand the cancellations
 @test length(poles(ss(G, minimal=true))) == 1 # The ss version with minimal realization does
 @test length(tzeros(G)) == 0 # tzeros converts to minimal ss relalization
 @test minreal(ss(G)).A ≈ [-2]
@@ -398,4 +400,16 @@ end
 @test length(tzeros(G)) == 3
 @test es(tzeros(G)) ≈ es(tzeros(big(1)G))
 
-end
+end
+
+
+## large TF poles and zeros
+G = ssrand(2,3,4)
+Gtf = tf(G)
+
+pss  = poles(G)
+zss  = tzeros(G)
+ptf  = poles(Gtf)
+ztf  = tzeros(Gtf)
+pzpk = poles(zpk(G))
+zzpk = tzeros(zpk(G))

lib/ControlSystemsBase/test/test_discrete.jl

@@ -140,7 +140,8 @@ p = poles(Gcl)
 # Test that all desired poles are in the closed-loop system
 @test norm(minimum(abs.((poles(tf(Bm,Am)) .- sort(p, by=imag)')), dims=2)) < 1e-6
 # Test that the observer poles are in the closed-loop system
-@test norm(minimum(abs.((poles(tf(1,Ao)) .- sort(p, by=imag)')), dims=2)) < 1e-6
+# @test norm(minimum(abs.((poles(tf(1,Ao)) .- sort(p, by=imag)')), dims=2)) < 1e-6
+@test sort(poles(tf(1,Ao)), by=LinearAlgebra.eigsortby)[2:end] ≈ sort(p, by=LinearAlgebra.eigsortby) # One pole is (correctly) cancelled in Gcl, causing the test above to fail
 
 
 pd = c2d_poly2poly(A, 0.1)

lib/ControlSystemsBase/test/test_pid_design.jl

@@ -52,7 +52,7 @@ Tf = 0.01
 
 # Different damping
 Ctf = pid(1,1,1, Tf=0.1, d = 1)
-@test all(p->imag(p) == 0, poles(Ctf))
+@test all(p->isapprox(imag(p), 0, atol=1e-6), poles(Ctf))
 Css = pid(1,1,1, Tf=0.1, d = 1, state_space=true)
 @test all(p->imag(p) == 0, poles(Css))
 @test tf(Css) ≈ Ctf