From 5357668db0f52f06388c31ea898bc7b70a622fde Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 14 Feb 2025 13:03:43 +0100 Subject: [PATCH] use zero computations form MatrixPencils (#972) * use zero computations form MatrixPencils * rm --- .../src/ControlSystemsBase.jl | 1 - lib/ControlSystemsBase/src/analysis.jl | 121 ++---------------- lib/ControlSystemsBase/src/plotting.jl | 2 +- lib/ControlSystemsBase/src/precompilation.jl | 2 + lib/ControlSystemsBase/test/test_analysis.jl | 15 ++- lib/ControlSystemsBase/test/test_complex.jl | 11 +- 6 files changed, 33 insertions(+), 119 deletions(-) diff --git a/lib/ControlSystemsBase/src/ControlSystemsBase.jl b/lib/ControlSystemsBase/src/ControlSystemsBase.jl index 090352853..ea4851fc7 100644 --- a/lib/ControlSystemsBase/src/ControlSystemsBase.jl +++ b/lib/ControlSystemsBase/src/ControlSystemsBase.jl @@ -39,7 +39,6 @@ export LTISystem, place, place_knvd, # Model Simplification - reduce_sys, sminreal, minreal, balreal, diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index bbf35e762..e9e787f21 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/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 """ diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl index c3f0ddd65..df3934e47 100644 --- a/lib/ControlSystemsBase/src/plotting.jl +++ b/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])) diff --git a/lib/ControlSystemsBase/src/precompilation.jl b/lib/ControlSystemsBase/src/precompilation.jl index 5aa8ad6d0..9bf7fdc83 100644 --- a/lib/ControlSystemsBase/src/precompilation.jl +++ b/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) diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index 507767147..aadf8acee 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/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 ## diff --git a/lib/ControlSystemsBase/test/test_complex.jl b/lib/ControlSystemsBase/test/test_complex.jl index 0a44493ac..e9e709096 100644 --- a/lib/ControlSystemsBase/test/test_complex.jl +++ b/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 +