From 6bb6eb40376eee031b0387ba42edf6b0ed9062ea Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 4 Jul 2025 09:49:24 +0200 Subject: [PATCH 1/2] use new generic zeros from MatrixPencils --- lib/ControlSystemsBase/Project.toml | 2 +- lib/ControlSystemsBase/src/analysis.jl | 210 ++++++++++++------------- 2 files changed, 106 insertions(+), 106 deletions(-) diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml index 30412b1a5..c8238b08e 100644 --- a/lib/ControlSystemsBase/Project.toml +++ b/lib/ControlSystemsBase/Project.toml @@ -36,7 +36,7 @@ ImplicitDifferentiation = "0.7.2" LinearAlgebra = "<0.0.1, 1" MacroTools = "0.5" MatrixEquations = "1, 2.1" -MatrixPencils = "1.6" +MatrixPencils = "1.8.3" Polynomials = "3.0, 4.0" PrecompileTools = "1" Printf = "<0.0.1, 1" diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index b1000b078..81c6c8ad7 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/lib/ControlSystemsBase/src/analysis.jl @@ -276,7 +276,7 @@ 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}(), balance=true, kwargs...) where {T <: BlasFloat, 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, E} if balance A, B, C = balance_statespace(A, B, C) end @@ -289,110 +289,110 @@ 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 +# 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) From e2d204c3d263a8f2ede8264519499fa250892fac Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 24 Jul 2025 15:32:25 +0200 Subject: [PATCH 2/2] update tests --- lib/ControlSystemsBase/test/test_analysis.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index 667247d6d..a250d6665 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/lib/ControlSystemsBase/test/test_analysis.jl @@ -1,6 +1,6 @@ @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 ## # Examples from the Emami-Naeini & Van Dooren Paper @@ -23,7 +23,7 @@ ex_3 = ss(A, B, C, D) @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]) atol=1e-6 # 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; @@ -58,7 +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 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,8 +79,8 @@ C = [0 0 0 1 0 0] 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 tzeros(ex_8) ≈ [-1.0, -1.0] atol=1e-6 +@test tzeros(big(1)ex_8) ≈ [-1.0, -1.0] atol=1e-12 @test ControlSystemsBase.count_integrators(ex_8) == 0 # Example 9 @@ -103,7 +103,7 @@ D = [0 0; 0 0; 0 0] ex_11 = ss(A, B, C, D) -@test tzeros(ex_11) ≈ [4.0, -3.0] +@test tzeros(ex_11) ≈ [4.0, -3.0] atol=1e-5 @test tzeros(big(1)ex_11) ≈ [4.0, -3.0] # Test for multiple zeros, siso tf @@ -405,7 +405,6 @@ end @test length(tzeros(G)) == 3 @test es(tzeros(G)) ≈ es(tzeros(big(1)G)) -end ## large TF poles and zeros