Skip to content

Commit

Permalink
Enhanced definition of DescriptorStateSpace object
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasvarga
andreasvarga committed Jun 20, 2024
1 parent e326f38 commit cb6fd45
Show file tree
Hide file tree
Showing 7 changed files with 94 additions and 30 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "DescriptorSystems"
uuid = "a81e2ce2-54d1-11eb-2c75-db236b00f339"
authors = ["Andreas Varga <[email protected]>"]
version = "1.4.1"
version = "1.4.2"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
4 changes: 4 additions & 0 deletions ReleaseNotes.md
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand Down
15 changes: 9 additions & 6 deletions src/DescriptorSystems.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down
14 changes: 7 additions & 7 deletions src/analysis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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*ϵ)
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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))

Expand Down
8 changes: 6 additions & 2 deletions src/dss.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)))
Expand Down
4 changes: 3 additions & 1 deletion src/dsutils.jl
View file
Original file line number Diff line number Diff line change
@@ -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))
Expand Down
77 changes: 64 additions & 13 deletions src/types/DescriptorStateSpace.jl
View file
Original file line number Diff line number Diff line change
@@ -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))}
Expand Down Expand Up @@ -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

This comment has been minimized.

Sign in to view

Copy link
@Keno

Keno Aug 1, 2024

As written, this forces all the matrices to be of the same type. I have a case where they are not:

  LoadError: MethodError: no method matching (DescriptorStateSpace{Float64, ET, MT} where {ET<:Union{LinearAlgebra.UniformScaling{Bool}, AbstractMatrix{Float64}}, MT<:AbstractMatrix{Float64}})(::Matrix{Float64}, ::SparseArrays.SparseMatrixCSC{Float64, Int64}, ::Matrix{Float64}, ::LinearAlgebra.Diagonal{Float64, Vector{Float64}}, ::Matrix{Float64}, ::Int64)
  The type `DescriptorStateSpace{Float64, ET, MT} where {ET<:Union{LinearAlgebra.UniformScaling{Bool}, AbstractMatrix{Float64}}, MT<:AbstractMatrix{Float64}}` exists, but no method is defined for this combination of argument types when trying to construct it.

  Closest candidates are:
    (DescriptorStateSpace{T, ET, MT} where {ET<:Union{LinearAlgebra.UniformScaling{Bool}, AbstractMatrix{T}}, MT<:AbstractMatrix{T}})(::MT, ::Union{LinearAlgebra.UniformScaling{Bool}, AbstractMatrix{T}}, ::MT, ::MT, ::MT, ::Real) where {T, MT<:AbstractMatrix{T}}
     @ DescriptorSystems ~/.julia/packages/DescriptorSystems/J8hbA/src/types/DescriptorStateSpace.jl:75

  Stacktrace:
    [1] dss(A::Matrix{Float64}, E::SparseArrays.SparseMatrixCSC{Float64, Int64}, B::Vector{Float64}, C::LinearAlgebra.Diagonal{Bool, Vector{Bool}}, D::Int64; Ts::Int64, check_reg::Bool, atol::Float64, atol1::Float64, atol2::Float64, rtol::Float64)
      @ DescriptorSystems ~/.julia/packages/DescriptorSystems/J8hbA/src/dss.jl:62
    [2] dss(A::Matrix{Float64}, E::SparseArrays.SparseMatrixCSC{Float64, Int64}, B::Vector{Float64}, C::LinearAlgebra.Diagonal{Bool, Vector{Bool}}, D::Int64)
      @ DescriptorSystems ~/.julia/packages/DescriptorSystems/J8hbA/src/dss.jl:52

This worked (although suboptimally) before this commit, because they were all coerced to Matrix{Float64}. Now it gives above error. I think this needs either more type parameters or additional methods that decide which matrix type is controlling and coerces the others.
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)

Expand All @@ -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)
Expand Down

3 comments on commit cb6fd45

@andreasvarga
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/109443

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.4.2 -m "<description of version>" cb6fd45de5a8cbe5f9aa9d52daeccac90c7dfb50
git push origin v1.4.2

@andreasvarga
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thank you bringing this issue to my attention. I could ensure to automatically convert the basic structured matrices (e.g., Diagonal, UpperTriangular, LowerTriangular, UpperHessenberg) to full matrices. Doing this for sparse matrices, would however imply to include the SparseArrays package among the required packages. Also this conversion could lead to inefficiency. I would like to avoid this, because there is practically no support (yet) for sparse matrices in the DescriptorSystems. The required conversion can be done explicitly before calling the function dss, so the responsability for inefficiency remains on the user side.

I realize now that besides nice opportunities, the new definition of the descriptor system object brings also some unexpected troubles. I would like to minimize the impact of this change as much as possible.

Please sign in to comment.