From 88f31874af81cd49f1a48dd775eb0168fddab625 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 25 Feb 2025 08:46:22 +0100 Subject: [PATCH] enable tzeros for BigFloat --- lib/ControlSystemsBase/Project.toml | 4 +- .../src/ControlSystemsBase.jl | 2 +- lib/ControlSystemsBase/src/analysis.jl | 115 +++++++++++++++++- .../src/types/conversion.jl | 8 +- lib/ControlSystemsBase/test/test_analysis.jl | 19 ++- 5 files changed, 138 insertions(+), 10 deletions(-) diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml index b00b788ef..8f377288f 100644 --- a/lib/ControlSystemsBase/Project.toml +++ b/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"] diff --git a/lib/ControlSystemsBase/src/ControlSystemsBase.jl b/lib/ControlSystemsBase/src/ControlSystemsBase.jl index ea4851fc7..52b7fdcfb 100644 --- a/lib/ControlSystemsBase/src/ControlSystemsBase.jl +++ b/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) diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index e9e787f21..38a177fac 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/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) diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl index b3777659a..af0d5bc30 100644 --- a/lib/ControlSystemsBase/src/types/conversion.jl +++ b/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 diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index aadf8acee..e29250e76 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/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 \ No newline at end of file