Merge pull request #976 from JuliaControl/bigfloatzeros

enable `tzeros` for `BigFloat`

lib/ControlSystemsBase/Project.toml

@@ -51,7 +51,7 @@ ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
-GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -60,4 +60,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "GenericLinearAlgebra", "GR", "Plots", "SparseArrays", "StaticArrays"]
+test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "GenericSchur", "GR", "Plots", "SparseArrays", "StaticArrays"]

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -242,7 +242,7 @@ function __init__()
             printstyled(io, "install and load ControlSystems.jl, or pass the keyword method = :zoh", color=:green, bold=true)
             print(io, " for automatic discretization (applicable to systems without delays or nonlinearities only).")
         elseif exc.f ∈ (eigvals!, ) && argtypes[1] <: AbstractMatrix{<:Number}
-            printstyled(io, "\nComputing eigenvalues of a matrix with exotic element types may require `using GenericLinearAlgebra`.", color=:green, bold=true)
+            printstyled(io, "\nComputing eigenvalues of a matrix with exotic element types may require `using GenericSchur`.", color=:green, bold=true)
         end
         plots_id = Base.PkgId(UUID("91a5bcdd-55d7-5caf-9e0b-520d859cae80"), "Plots")
         if exc.f isa Function && nameof(exc.f) === :plot && parentmodule(argtypes[1]) == @__MODULE__() && !haskey(Base.loaded_modules, plots_id)

lib/ControlSystemsBase/src/analysis.jl

@@ -4,6 +4,8 @@
 Compute the poles of system `sys`.
 
 Note: Poles with multiplicity `n > 1` may suffer numerical inaccuracies on the order `eps(numeric_type(sys))^(1/n)`, i.e., a double pole in a system with `Float64` coefficients may be computed with an error of about `√(eps(Float64)) ≈ 1.5e-8`.
+
+To compute the poles of a system with non-BLAS floats, such as `BigFloat`, install and load the package `GenericSchur.jl` before calling `poles`.
 """
 poles(sys::AbstractStateSpace) = eigvalsnosort(sys.A)
 
@@ -243,6 +245,8 @@ Compute the invariant zeros of the system `sys`. If `sys` is a minimal
 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.
+
+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)
     if issiso(sys)
@@ -260,7 +264,11 @@ function tzeros(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::Abst
     tzeros(A2, B2, C2, D2)
 end
 
-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}
+function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}; extra::Val{E} = Val{false}(), balance=true, kwargs...) where {T <: BlasFloat, E}
+    if balance
+        A, B, C = balance_statespace(A, B, C)
+    end
+
     (z, iz, KRInfo) = MatrixPencils.spzeros(A, I, B, C, D; kwargs...)
     if E
         return (z, iz, KRInfo)
@@ -269,6 +277,111 @@ function tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}
     end
 end
 
+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}
+    isempty(A) && return complex(T)[]
+
+    # Balance the system
+    A, B, C = balance_statespace(A, B, C; verbose=false)
+
+    # 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)
+
+    # 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]
+        Z = schur(complex.(Af), complex.(Bf)) # Schur instead of eigvals to handle BigFloat
+        zs = Z.values
+        ControlSystemsBase._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, :])
+    end
+    mrank = sum(norms .> meps)
+    return mrank
+end
+
 """
     relative_gain_array(G, w::AbstractVector)
     relative_gain_array(G, w::Number)

lib/ControlSystemsBase/src/types/conversion.jl

@@ -159,11 +159,11 @@ The inverse of `sysb, T = balance_statespace(sys)` is given by `similarity_trans
 
 This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see [`balreal`](@ref)
 """
-function balance_statespace(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, perm::Bool=false)
+function balance_statespace(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, perm::Bool=false; verbose=true)
     try
         return _balance_statespace(A,B,C, perm)
     catch
-        @warn "Unable to balance state-space, returning original system" maxlog=10
+        verbose && @warn "Unable to balance state-space, returning original system" maxlog=10
         return A,B,C,convert(typeof(A), I(size(A, 1)))
     end
 end
@@ -175,8 +175,8 @@ end
 #     balance_statespace(A2, B2, C2, perm)
 # end
 
-function balance_statespace(sys::S, perm::Bool=false) where S <: AbstractStateSpace
-    A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
+function balance_statespace(sys::S, perm::Bool=false; kwargs...) where S <: AbstractStateSpace
+    A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm; kwargs...)
     ss(A,B,C,sys.D,sys.timeevol), T
 end
 

lib/ControlSystemsBase/test/test_analysis.jl

@@ -1,5 +1,5 @@
-@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
+@test_throws MethodError poles(big(1.0)*ssrand(1,1,1)) # This errors before loading GenericSchur
+using GenericSchur # Required to compute eigvals (in tzeros and poles) of a matrix with exotic element types
 @testset "test_analysis" begin
 es(x) = sort(x, by=LinearAlgebra.eigsortby)
 ## tzeros ##
@@ -58,6 +58,7 @@ C = [0 -1 0]
 D = [0]
 ex_6 = ss(A, B, C, D)
 @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_broken tzeros(big(1.0)ex_6) == [2]
 @test ControlSystemsBase.count_integrators(ex_6) == 2
 
 @test ss(A, [0 0 1]', C, D) == ex_6
@@ -79,6 +80,7 @@ D = [0]
 ex_8 = ss(A, B, C, D)
 # TODO : there may be a way to improve the precision of this example.
 @test tzeros(ex_8) ≈ [-1.0, -1.0] atol=1e-7
+@test tzeros(big(1)ex_8) ≈ [-1.0, -1.0] atol=1e-7
 @test ControlSystemsBase.count_integrators(ex_8) == 0
 
 # Example 9
@@ -102,6 +104,7 @@ D = [0 0;
      0 0]
 ex_11 = ss(A, B, C, D)
 @test tzeros(ex_11) ≈ [4.0, -3.0]
+@test tzeros(big(1)ex_11) ≈ [4.0, -3.0]
 
 # Test for multiple zeros, siso tf
 s = tf("s")
@@ -368,4 +371,16 @@ P = let
 end
 @test ControlSystemsBase.count_integrators(P) == 2
 
+## Difficult test case for zeros
+G = let
+     tempA = [-0.6841991610512457 -0.0840213470263692 -0.0004269818661494616 -2.7625001165862086e-18; 2.081491248616774 0.0 0.0 8.404160870560225e-18; 0.0 24.837541148074962 0.12622006230897712 0.0; -1.2211265763794115e-14 -2.778983834717109e8 -1.4122312296634873e6 -4.930380657631326e-32]
+     tempB = [-0.5316255605902501; 2.0811471051085637; -45.068824982602656; 5.042589978197361e8;;]
+     tempC = [0.0 0.0 0.954929658551372 0.0]
+     tempD = [0.0;;]
+     ss(tempA, tempB, tempC, tempD)
+end
+
+@test length(tzeros(G)) == 3
+@test es(tzeros(G)) ≈ es(tzeros(big(1)G))
+
 end