use zero computations form MatrixPencils (#972)

* use zero computations form MatrixPencils

* rm

Author: baggepinnen

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -39,7 +39,6 @@ export  LTISystem,
         place,
         place_knvd,
         # Model Simplification
-        reduce_sys,
         sminreal,
         minreal,
         balreal,

lib/ControlSystemsBase/src/analysis.jl

@@ -237,131 +237,36 @@ dampreport(sys::LTISystem) = dampreport(stdout, sys)
 
 """
     tzeros(sys)
+    tzeros(sys::AbstractStateSpace; extra=Val(false))
 
 Compute the invariant zeros of the system `sys`. If `sys` is a minimal
-realization, these are also the transmission zeros."""
+realization, these are also the transmission zeros.
+
+If `sys` is a state-space system the function has additional keyword arguments, see [`?ControlSystemsBase.MatrixPencils.spzeros`](https://andreasvarga.github.io/MatrixPencils.jl/dev/sklfapps.html#MatrixPencils.spzeros) for more details. If `extra = Val(true)`, the function returns `z, iz, KRInfo` where `z` are the transmission zeros, information on the multiplicities of infinite zeros in `iz` and information on the Kronecker-structure in the KRInfo object. The number of infinite zeros is the sum of the components of iz.
+"""
 function tzeros(sys::TransferFunction)
     if issiso(sys)
         return tzeros(sys.matrix[1,1])
     else
-        return tzeros(ss(sys))
+        return tzeros(ss(sys, minimal=true))
     end
 end
 
-# Implements the algorithm described in:
-# Emami-Naeini, A. and P. Van Dooren, "Computation of Zeros of Linear
-# Multivariable Systems," Automatica, 18 (1982), pp. 415–430.
-#
-# Note that this returns either Vector{ComplexF32} or Vector{Float64}
-tzeros(sys::AbstractStateSpace) = tzeros(sys.A, sys.B, sys.C, sys.D)
+tzeros(sys::AbstractStateSpace; kwargs...) = tzeros(sys.A, sys.B, sys.C, sys.D; kwargs...)
 # Make sure everything is BlasFloat
 function tzeros(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)
     T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
     A2, B2, C2, D2, _ = promote(A,B,C,D, fill(zero(T)/one(T),0,0)) # If Int, we get Float64
     tzeros(A2, B2, C2, D2)
 end
-function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}) where {T <: Union{AbstractFloat,Complex{<:AbstractFloat}}#= For eps(T) =#}
-    # Balance the system
-    A, B, C = balance_statespace(A, B, C)
-
-    # Compute a good tolerance
-    meps = 10*eps(real(T))*norm([A B; C D])
-
-    # Step 1:
-    A_r, B_r, C_r, D_r = reduce_sys(A, B, C, D, meps)
-
-    # Step 2: (conjugate transpose should be avoided since single complex zeros get conjugated)
-    A_rc, B_rc, C_rc, D_rc = reduce_sys(copy(transpose(A_r)), copy(transpose(C_r)), copy(transpose(B_r)), copy(transpose(D_r)), meps)
-    isempty(A) && return complex(T)[]
-
-    # Step 3:
-    # Compress cols of [C D] to [0 Df]
-    mat = [C_rc D_rc]
-    Wr = qr(mat').Q * I
-    W = reverse(Wr, dims=2)
-    mat = mat*W
-    if fastrank(mat', meps) > 0
-        nf = size(A_rc, 1)
-        m = size(D_rc, 2)
-        Af = ([A_rc B_rc] * W)[1:nf, 1:nf]
-        Bf = ([Matrix{T}(I, nf, nf) zeros(nf, m)] * W)[1:nf, 1:nf]
-        zs = eigvalsnosort(Af, Bf)
-        _fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs
-    else
-        zs = complex(T)[]
-    end
-    return zs
-end
 
-
-"""
-    reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
-Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed
-A, B, C, D matrices. These are empty if there are no zeros.
-"""
-function reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
-    T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
-    Cbar, Dbar = C, D
-    if isempty(A)
-        return A, B, C, D
-    end
-    while true
-        # Compress rows of D
-        U = qr(D).Q
-        D = U'D
-        C = U'C
-        sigma = fastrank(D, meps)
-        Cbar = C[1:sigma, :]
-        Dbar = D[1:sigma, :]
-        Ctilde = C[(1 + sigma):end, :]
-        if sigma == size(D, 1)
-            break
-        end
-
-        # Compress columns of Ctilde
-        V = reverse(qr(Ctilde').Q * I, dims=2)
-        Sj = Ctilde*V
-        rho = fastrank(Sj', meps)
-        nu = size(Sj, 2) - rho
-
-        if rho == 0
-            break
-        elseif nu == 0
-            # System has no zeros, return empty matrices
-            A = B = Cbar = Dbar = Matrix{T}(undef, 0,0)
-            break
-        end
-        # Update System
-        n, m = size(B)
-        Vm = [V zeros(T, n, m); zeros(T, m, n) Matrix{T}(I, m, m)] # I(m) is not used for type stability reasons (as of julia v1.7)
-        if sigma > 0
-            M = [A B; Cbar Dbar]
-            Vs = [copy(V') zeros(T, n, sigma) ; zeros(T, sigma, n) Matrix{T}(I, sigma, sigma)]
-        else
-            M = [A B]
-            Vs = copy(V')
-        end
-        sigma, rho, nu
-        M = Vs * M * Vm
-        A = M[1:nu, 1:nu]
-        B = M[1:nu, (nu + rho + 1):end]
-        C = M[(nu + 1):end, 1:nu]
-        D = M[(nu + 1):end,  (nu + rho + 1):end]
-    end
-    return A, B, Cbar, Dbar
-end
-
-# Determine the number of non-zero rows, with meps as a tolerance. For an
-# upper-triangular matrix, this is a good proxy for determining the row-rank.
-function fastrank(A::AbstractMatrix, meps::Real)
-    n, m = size(A)
-    if n*m == 0     return 0    end
-    norms = Vector{real(eltype(A))}(undef, n)
-    for i = 1:n
-        norms[i] = norm(A[i, :])
+function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}; extra::Val{E} = Val{false}(), kwargs...) where {T <: Union{AbstractFloat,Complex{<:AbstractFloat}}, E}
+    (z, iz, KRInfo) = MatrixPencils.spzeros(A, I, B, C, D; kwargs...)
+    if E
+        return (z, iz, KRInfo)
+    else
+        return filter(isfinite, z)
     end
-    mrank = sum(norms .> meps)
-    return mrank
 end
 
 """

lib/ControlSystemsBase/src/plotting.jl

@@ -781,7 +781,7 @@ marginplot
                 end
                 titles[j,i,1,1] *= "["*join([Printf.@sprintf("%3.2g",v) for v in gm],", ")*"] "
                 titles[j,i,1,2] *= "["*join([Printf.@sprintf("%3.2g",v) for v in wgm],", ")*"] "
-                titles[j,i,2,1] *=  "["*join([Printf.@sprintf("%3.2g",v) for v in pm],", ")*"] "
+                titles[j,i,2,1] *=  "["*join([Printf.@sprintf("%3.2g°",v) for v in pm],", ")*"] "
                 titles[j,i,2,2] *=  "["*join([Printf.@sprintf("%3.2g",v) for v in wpm],", ")*"] "
 
                 subplot := min(s2i((plotphase ? (2i-1) : i),j), prod(plotattributes[:layout]))

lib/ControlSystemsBase/src/precompilation.jl

@@ -22,8 +22,10 @@ PrecompileTools.@setup_workload begin
         end
         G = tf(1.0, [1.0, 1])
         ss(G)
+        ss(G, minimal=true)
         G = tf(1.0, [1.0, 1], 1)
         ss(G)
+        ss(G, minimal=true)
 
         # Pdel = P*delay(1.0)
         # pade(Pdel, 2)

lib/ControlSystemsBase/test/test_analysis.jl

@@ -1,6 +1,7 @@
 @test_throws MethodError poles(big(1.0)*ssrand(1,1,1)) # This errors before loading GenericLinearAlgebra
 using GenericLinearAlgebra # Required to compute eigvals of a matrix with exotic element types
 @testset "test_analysis" begin
+es(x) = sort(x, by=LinearAlgebra.eigsortby)
 ## tzeros ##
 # Examples from the Emami-Naeini & Van Dooren Paper
 # Example 3
@@ -19,10 +20,10 @@ D = [1 0;
 
 ex_3 = ss(A, B, C, D)
 @test ControlSystemsBase.count_integrators(ex_3) == 6
-@test tzeros(ex_3) ≈ [0.3411639019140099 + 1.161541399997252im,
+@test es(tzeros(ex_3)) ≈ es([0.3411639019140099 + 1.161541399997252im,
                              0.3411639019140099 - 1.161541399997252im,
                              0.9999999999999999 + 0.0im,
-                             -0.6823278038280199 + 0.0im]
+                             -0.6823278038280199 + 0.0im])
 # Example 4
 A = [-0.129    0.0   0.396e-1  0.25e-1    0.191e-1;
      0.329e-2  0.0  -0.779e-4  0.122e-3  -0.621;
@@ -38,7 +39,7 @@ C = [1 0 0 0 0;
      0 1 0 0 0]
 D = zeros(2, 2)
 ex_4 = ss(A, B, C, D)
-@test tzeros(ex_4) ≈ [-0.06467751189940692,-0.3680512036036696]
+@test es(tzeros(ex_4)) ≈ es([-0.06467751189940692,-0.3680512036036696])
 @test ControlSystemsBase.count_integrators(ex_4) == 1
 
 # Example 5
@@ -56,7 +57,7 @@ B = [0; 0; 1]
 C = [0 -1 0]
 D = [0]
 ex_6 = ss(A, B, C, D)
-@test tzeros(ex_6) == Float64[]
+@test tzeros(ex_6) == [2] # From paper: "Our algorithm will extract the singular part of S(A) and will yield a regular pencil containing the single zero at 2."
 @test ControlSystemsBase.count_integrators(ex_6) == 2
 
 @test ss(A, [0 0 1]', C, D) == ex_6
@@ -194,9 +195,11 @@ 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
-@test length(tzeros(G)) == 1
+@test_broken length(poles(G)) == 1 # The tf poles don't 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]
+@test ss(G, minimal=true).A ≈ [-2]
 
 
 ## MARGIN ##

lib/ControlSystemsBase/test/test_complex.jl

@@ -34,9 +34,14 @@ C_2 = zpk([-1+im], [], 1.0+1im)
 
 s = tf("s");
 @test tzeros(ss(-1, 1, 1, 1.0im)) ≈ [-1.0 + im] rtol=1e-15
-@test tzeros(ss([-1.0-im 1-im; 2 0], [2; 0], [-1+1im -0.5-1.25im], 1)) ≈ [-1-2im, 2-im]
+@test tzeros(ss([-1.0-im 1-im; 2 0], [2; 0], [-1+1im -0.5-1.25im], 1)) ≈ [2-im, -1-2im]
 
-@test tzeros(ss((s-2.0-1.5im)^3/(s+1+im)/(s+2)^3)) ≈ fill(2.0 + 1.5im, 3) rtol=1e-4
-@test tzeros(ss((s-2.0-1.5im)*(s-3.0)/(s+1+im)/(s+2)^2)) ≈ [2.0 + 1.5im, 3.0] rtol=1e-14
+sys = ss((s-2.0-1.5im)^3/(s+1+im)/(s+2)^3)
+@test tzeros(sys) ≈ fill(2.0 + 1.5im, 3) rtol=1e-4
+@test tzeros(ss((s-2.0-1.5im)*(s-3.0)/(s+1+im)/(s+2)^2)) ≈ [3.0, 2.0 + 1.5im] rtol=1e-14
+
+z,iv,info = tzeros(sys, extra=Val(true))
+@test iv == [1]
 
 end
+