From cb6fd45de5a8cbe5f9aa9d52daeccac90c7dfb50 Mon Sep 17 00:00:00 2001 From: DSVarga Date: Thu, 20 Jun 2024 13:17:10 +0200 Subject: [PATCH] Enhanced definition of DescriptorStateSpace object --- Project.toml | 2 +- ReleaseNotes.md | 4 ++ src/DescriptorSystems.jl | 15 +++--- src/analysis.jl | 14 +++--- src/dss.jl | 8 +++- src/dsutils.jl | 4 +- src/types/DescriptorStateSpace.jl | 77 +++++++++++++++++++++++++------ 7 files changed, 94 insertions(+), 30 deletions(-) diff --git a/Project.toml b/Project.toml index 7cea2f3..c342fbf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DescriptorSystems" uuid = "a81e2ce2-54d1-11eb-2c75-db236b00f339" authors = ["Andreas Varga "] -version = "1.4.1" +version = "1.4.2" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/ReleaseNotes.md b/ReleaseNotes.md index e6d511d..8304001 100644 --- a/ReleaseNotes.md +++ b/ReleaseNotes.md @@ -1,5 +1,9 @@ # Release Notes +## Version 1.4.2 + +Enhanced definition of the `DescriptorStateSpace` object. + ## Version 1.4.1 Patch release with enhanced functions for the computation of system norms. diff --git a/src/DescriptorSystems.jl b/src/DescriptorSystems.jl index 3eb2ec8..fe0ef9d 100644 --- a/src/DescriptorSystems.jl +++ b/src/DescriptorSystems.jl @@ -1,16 +1,19 @@ module DescriptorSystems +import Base: *, +, -, /, \, ^, (==), convert, eltype, hcat, hvcat, inv, isapprox, iszero, + lastindex, length, ndims, nothing, one, print, promote_op, require_one_based_indexing, + show, size, vcat, zero +import LinearAlgebra: BlasComplex, BlasFloat, BlasReal, adjoint, copy_oftype, hcat, + normalize, opnorm, rdiv!, transpose +import MatrixPencils: isregular, rmeval +import Polynomials: + AbstractPolynomial, AbstractRationalFunction, degree, isconstant, poles, pqs, variable + using LinearAlgebra using MatrixEquations using MatrixPencils using Polynomials using Random - -import LinearAlgebra: BlasFloat, BlasReal, BlasComplex, copy_oftype, transpose, adjoint, opnorm, normalize, rdiv!, hcat -import Base: +, -, *, /, \, (==), (!=), ^, isapprox, iszero, convert, promote_op, size, length, ndims, - hcat, vcat, hvcat, inv, show, lastindex, require_one_based_indexing, print, show, one, zero, eltype -import MatrixPencils: isregular, rmeval -import Polynomials: AbstractRationalFunction, AbstractPolynomial, poles, isconstant, variable, degree, pqs isdefined(Polynomials,:order) && (import Polynomials: order) export DescriptorStateSpace, AbstractDescriptorStateSpace, dss, dssdata, rdss, rss, iszero, order diff --git a/src/analysis.jl b/src/analysis.jl index 87115a4..b4cdb21 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -64,7 +64,7 @@ function gzero(sys::DescriptorStateSpace;fast = false, prescale = gbalqual(sys) # pzeros([SYS.A SYS.B; SYS.C SYS.D], [SYS.E zeros(SYS.nx,SYS.nu); zeros(SYS.ny,SYS.nx+SYS.nu)]; fast = fast, atol1 = atol1, # atol2 = atol2, rtol = rtol )[1] SYS = prescale ? gprescale(sys)[1] : sys - return spzeros(dssdata(SYS)...; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)[1] + return MatrixPencils.spzeros(dssdata(SYS)...; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)[1] end """ val = gpole(sys; fast = false, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) @@ -206,7 +206,7 @@ function gzeroinfo(sys::DescriptorStateSpace{T}; prescale = gbalqual(sys) > 1000 rtol::Real = sys.nx*eps(real(float(one(T))))*iszero(min(atol1,atol2)), offset::Real = sqrt(eps(float(real(T)))) ) where T SYS = prescale ? gprescale(sys)[1] : sys - val, miz, krinfo = spzeros(dssdata(SYS)...; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol) + val, miz, krinfo = MatrixPencils.spzeros(dssdata(SYS)...; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol) nfsz, nfsbz, nfuz = eigvals_info(val[isfinite.(val)], smarg, SYS.Ts != 0, offset) niev = sum(krinfo.id) nisev = niev == 0 ? 0 : krinfo.id[1] @@ -966,11 +966,11 @@ function norminfd(a, e, b, c, d, ft0, Ts, tol) # Compute lower estimate GMIN as max. gain over the selected frequencies for i = 1:length(testz) - z = testz[i] - bct = copy(bc) - desc ? ldiv!(UpperHessenberg(ac-z*ec),bct) : ldiv!(ac,bct,shift = -z) - gw = opnorm(dc-cc*bct) - gw > gmin && (gmin = gw; compl ? fpeak = angle(z) : fpeak = abs(angle(z))) + z = testz[i] + bct = copy(bc) + desc ? ldiv!(UpperHessenberg(ac-z*ec),bct) : ldiv!(ac,bct,shift = -z) + gw = opnorm(dc-cc*bct) + gw > gmin && (gmin = gw; compl ? fpeak = angle(z) : fpeak = abs(angle(z))) end gmin == 0 && (return TR(0), TR(0)) diff --git a/src/dss.jl b/src/dss.jl index f08ae99..30b9fe1 100644 --- a/src/dss.jl +++ b/src/dss.jl @@ -69,6 +69,10 @@ function dss(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray return DescriptorStateSpace{T}(to_matrix(T,A), I, to_matrix(T,B), to_matrix(T,C,p <=1), typeof(D) <: Number && iszero(D) ? zeros(T,p,m) : p <=1 ? to_matrix(T,D,m > p) : to_matrix(T,D), Ts) end +# function dss(A::SparseMatrixCSC, B::SparseMatrixCSC, C::SparseMatrixCSC, D::SparseMatrixCSC; Ts::Real = 0) +# T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D)) +# return DescriptorStateSpace{T}(A, I, B, C, D, Ts) +# end function dss(D::AbstractNumOrArray; Ts::Real = 0) D == [] && (T = Int64; return DescriptorStateSpace{T}(zeros(T,0,0),I,zeros(T,0,0),zeros(T,0,0),zeros(T,0,0),Ts)) p, m = size(D,1),size(D,2) @@ -127,8 +131,8 @@ for the nonzero elements of `F`, `G` and `H`. The default relative tolerance is of `A`, and `ϵ` is the machine epsilon of the element type of `A`. The keyword argument `atol` can be used to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`. """ -function dss(A::Matrix, E::Union{Matrix,UniformScaling}, B::Matrix, F::Union{Matrix,Missing}, - C::Matrix, G::Union{Matrix,Missing}, D::Matrix, H::Union{Matrix,Missing}; +function dss(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::AbstractMatrix, F::Union{AbstractMatrix,Missing}, + C::AbstractMatrix, G::Union{AbstractMatrix,Missing}, D::AbstractMatrix, H::Union{AbstractMatrix,Missing}; Ts::Real = 0, compacted::Bool = false, atol::Real = zero(real(eltype(A))), atol1::Real = atol, atol2::Real = atol, atol3::Real = atol, rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2,atol3))) diff --git a/src/dsutils.jl b/src/dsutils.jl index 3b45d94..1a68811 100644 --- a/src/dsutils.jl +++ b/src/dsutils.jl @@ -1,7 +1,9 @@ # The function to_matrix is a fork from ControlSystems.jl # covers the case when a vector represents the row of a matrix to_matrix(T, A::AbstractVector, wide::Bool = false) = wide ? Matrix{T}(reshape(A, 1, length(A))) : Matrix{T}(reshape(A, length(A), 1)) -to_matrix(T, A::AbstractMatrix, wide::Bool = false) = Matrix{T}(A) # Fallback +#to_matrix(T, A::AbstractMatrix, wide::Bool = false) = Matrix{T}(A) # Fallback +to_matrix(T, A::AbstractMatrix, wide::Bool = false) = T.(A) # Fallback +# to_matrix(T, A::SparseMatrixCSC, wide::Bool = false) = T.(A) # Fallback for sparse matrices to_matrix(T, A::Number, wide::Bool = true) = fill(T(A), 1, 1) # Handle Adjoint Matrices to_matrix(T, A::Adjoint{R, MT}, wide::Bool = false) where {R<:Number, MT<:AbstractMatrix} = to_matrix(T, MT(A)) diff --git a/src/types/DescriptorStateSpace.jl b/src/types/DescriptorStateSpace.jl index 9650e1f..9d99a85 100644 --- a/src/types/DescriptorStateSpace.jl +++ b/src/types/DescriptorStateSpace.jl @@ -1,7 +1,8 @@ #const AbstractNumOrArray = Union{AbstractVecOrMat{T},Number} where {T <: Number} const AbstractNumOrArray = Union{AbstractVecOrMat,Number} #const ETYPE{T} = Union{Matrix{T},UniformScaling} -const ETYPE{T} = Union{Matrix{T},UniformScaling{Bool}} +#const ETYPE{T} = Union{Matrix{T},UniformScaling{Bool}} +const ETYPE{T} = Union{AbstractMatrix{T},UniformScaling{Bool}} promote_Etype(T0::Type, ::Type{UniformScaling{Bool}}, ::Type{UniformScaling{Bool}}) = UniformScaling{Bool} promote_Etype(T0::Type, T1::Type, ::Type{UniformScaling{Bool}}) = Matrix{promote_type(T0,eltype(T1))} promote_Etype(T0::Type, ::Type{UniformScaling{Bool}}, T1::Type) = Matrix{promote_type(T0,eltype(T1))} @@ -33,24 +34,73 @@ defined by the 4-tuple `SYS = (A-λE,B,C,D)`, then: The dimensions `nx`, `ny` and `nu` can be obtained as `SYS.nx`, `SYS.ny` and `SYS.nu`, respectively. """ -struct DescriptorStateSpace{T, ET <: ETYPE{T}} <: AbstractDescriptorStateSpace -#struct DescriptorStateSpace{T, ET <: Union{Matrix{T},UniformScaling}} <: AbstractDescriptorStateSpace - A::Matrix{T} +# struct DescriptorStateSpace{T, ET <: ETYPE{T}} <: AbstractDescriptorStateSpace +# #struct DescriptorStateSpace{T, ET <: Union{Matrix{T},UniformScaling}} <: AbstractDescriptorStateSpace +# A::Matrix{T} +# E::ET +# B::Matrix{T} +# C::Matrix{T} +# D::Matrix{T} +# Ts::Float64 +# #function DescriptorStateSpace{T}(A::Matrix{T}, E::Union{Matrix{T},UniformScaling}, +# function DescriptorStateSpace{T}(A::Matrix{T}, E::ETYPE{T}, +# B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, Ts::Real) where {T} +# dss_validation(A, E, B, C, D, Ts) +# new{T, typeof(E)}(A, E, B, C, D, Float64(Ts)) +# end +# end +# struct DescriptorStateSpace{T, ET <: ETYPE{T}} <: AbstractDescriptorStateSpace +# #struct DescriptorStateSpace{T, ET <: Union{Matrix{T},UniformScaling}} <: AbstractDescriptorStateSpace +# A::AbstractMatrix{T} +# E::ET +# B::AbstractMatrix{T} +# C::AbstractMatrix{T} +# D::AbstractMatrix{T} +# Ts::Float64 +# #function DescriptorStateSpace{T}(A::Matrix{T}, E::Union{Matrix{T},UniformScaling}, +# function DescriptorStateSpace{T}(A::AbstractMatrix{T}, E::ETYPE{T}, +# B::AbstractMatrix{T}, C::AbstractMatrix{T}, +# D::AbstractMatrix{T}, Ts::Real) where {T} +# dss_validation(A, E, B, C, D, Ts) +# new{T, typeof(E)}(A, E, B, C, D, Float64(Ts)) +# end +# end +struct DescriptorStateSpace{T, ET <: ETYPE{T}, MT<:AbstractMatrix{T}} <: AbstractDescriptorStateSpace + A::MT E::ET - B::Matrix{T} - C::Matrix{T} - D::Matrix{T} + B::MT + C::MT + D::MT Ts::Float64 - #function DescriptorStateSpace{T}(A::Matrix{T}, E::Union{Matrix{T},UniformScaling}, - function DescriptorStateSpace{T}(A::Matrix{T}, E::ETYPE{T}, - B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, Ts::Real) where {T} + function DescriptorStateSpace{T}(A::MT, E::ETYPE{T}, + B::MT, C::MT, D::MT, Ts::Real) where {T, MT<:AbstractMatrix{T}} dss_validation(A, E, B, C, D, Ts) - new{T, typeof(E)}(A, E, B, C, D, Float64(Ts)) + new{T, typeof(E), MT}(A, E, B, C, D, Float64(Ts)) end end + #function dss_validation(A::Matrix{T}, E::Union{Matrix{T},UniformScaling}, -function dss_validation(A::Matrix{T}, E::ETYPE{T}, - B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, Ts::Real) where T +# function dss_validation(A::Matrix{T}, E::ETYPE{T}, +# B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, Ts::Real) where T +# nx = LinearAlgebra.checksquare(A) +# (ny, nu) = size(D) + +# # Validate dimensions +# if typeof(E) <: Matrix +# nx == LinearAlgebra.checksquare(E) || error("A and E must have the same size") +# end +# size(B, 1) == nx || error("B must have the same row size as A") +# size(C, 2) == nx || error("C must have the same column size as A") +# nu == size(B, 2) || error("D must have the same column size as B") +# ny == size(C, 1) || error("D must have the same row size as C") + +# # Validate sampling time +# Ts >= 0 || Ts == -1 || error("Ts must be either a positive number, 0 +# (continuous system), or -1 (unspecified)") +# return nx, nu, ny +# end +function dss_validation(A::AbstractMatrix{T}, E::ETYPE{T}, + B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}, Ts::Real) where T nx = LinearAlgebra.checksquare(A) (ny, nu) = size(D) @@ -68,6 +118,7 @@ function dss_validation(A::Matrix{T}, E::ETYPE{T}, (continuous system), or -1 (unspecified)") return nx, nu, ny end + function Base.getproperty(sys::DescriptorStateSpace, d::Symbol) if d === :nx return size(getfield(sys, :A), 1)