From 03ce6a3c3f68b0dcd86f827e7241f32b5da2cf40 Mon Sep 17 00:00:00 2001 From: DSVarga Date: Sat, 10 Aug 2024 13:21:31 +0200 Subject: [PATCH] Some fixes to handle models with structured and sparse matrices. --- Project.toml | 9 ++++- ReleaseNotes.md | 4 ++ src/connections.jl | 63 ++++++++++++++++++++++++++++++- src/dss.jl | 8 ++-- src/dsutils.jl | 5 +++ src/types/DescriptorStateSpace.jl | 18 ++++++++- test/test_dss.jl | 26 +++++++++++++ 7 files changed, 123 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index d9322a2..69f282f 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.3" +version = "1.4.4" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -14,13 +14,18 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" LinearAlgebra = "1" MatrixEquations = "2.4" MatrixPencils = "1.8" +Measurements = "2.11" Polynomials = "2.0, 3, 4" Random = "1" +SparseArrays = "1" Test = "1" julia = "1.8" [extras] +Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Measurements", "SparseArrays", "Test"] + diff --git a/ReleaseNotes.md b/ReleaseNotes.md index 92f9c55..7c8884e 100644 --- a/ReleaseNotes.md +++ b/ReleaseNotes.md @@ -1,5 +1,9 @@ # Release Notes +## Version 1.4.4 + +Some fixes to handle structured and sparse matrices in the `dss` function. + ## Version 1.4.3 Some fixes of errors originating from the enhanced definition of the `DescriptorStateSpace` object. diff --git a/src/connections.jl b/src/connections.jl index df74eb7..d12523b 100644 --- a/src/connections.jl +++ b/src/connections.jl @@ -67,14 +67,43 @@ function append(A::Union{DescriptorStateSpace,AbstractNumOrArray,UniformScaling} return append(promote_to_systems(Ts, fill(1,length(A)), 1, promote_type(eltype.(A)...), A...)...) end +# function hcat(SYS1::DescriptorStateSpace{T1, TE1}, SYS2::DescriptorStateSpace{T2, TE2}) where {T1<:Number, TE1<:ETYPE, T2<:Number, TE2<:ETYPE} +# #function hcat(SYS1 :: DescriptorStateSpace, SYS2 :: DescriptorStateSpace) +# ny = SYS1.ny +# ny == size(SYS2, 1) || error("The systems must have the same output dimension") +# #T = promote_type(eltype(SYS1), eltype(SYS2)) +# T = promote_type(T1, T2) +# #TE = promote_type(TE1, TE2) +# TE = promote_Etype(T, TE1, TE2) +# Ts = promote_Ts(SYS1.Ts, SYS2.Ts) -function hcat(SYS1::DescriptorStateSpace{T1, TE1}, SYS2::DescriptorStateSpace{T2, TE2}) where {T1<:Number, TE1<:ETYPE, T2<:Number, TE2<:ETYPE} +# A = blockdiag(T.(SYS1.A), T.(SYS2.A)) + +# local E::TE +# if SYS1.E === I && SYS2.E === I +# #if SYS1.E == I && SYS2.E == I +# E = I +# elseif SYS1.E == I || SYS2.E == I +# blockdims = [size(SYS1.A,1), size(SYS2.A,1)] +# E = sblockdiag(blockdims, SYS1.E == I ? SYS1.E : T.(SYS1.E), SYS2.E == I ? SYS2.E : T.(SYS2.E)) +# else +# E = blockdiag(T.(SYS1.E), T.(SYS2.E)) +# end +# B = blockdiag(T.(SYS1.B), T.(SYS2.B)) +# C = [ T.(SYS1.C) T.(SYS2.C)] +# D = [ T.(SYS1.D) T.(SYS2.D)] +# return DescriptorStateSpace{T}(A, E, B, C, D, Ts) +# end + +function hcat(SYS1::DescriptorStateSpace{T1}, SYS2::DescriptorStateSpace{T2}) where {T1<:Number, T2<:Number} #function hcat(SYS1 :: DescriptorStateSpace, SYS2 :: DescriptorStateSpace) ny = SYS1.ny ny == size(SYS2, 1) || error("The systems must have the same output dimension") #T = promote_type(eltype(SYS1), eltype(SYS2)) T = promote_type(T1, T2) #TE = promote_type(TE1, TE2) + TE1 = typeof(SYS1.E) + TE2 = typeof(SYS2.E) TE = promote_Etype(T, TE1, TE2) Ts = promote_Ts(SYS1.Ts, SYS2.Ts) @@ -179,13 +208,43 @@ function Base.vcat(systems::DescriptorStateSpace...) return DescriptorStateSpace{T}(A, E, B, C, D, Ts) end -function vcat(SYS1::DescriptorStateSpace{T1, TE1}, SYS2::DescriptorStateSpace{T2, TE2}) where {T1<:Number, TE1<:ETYPE, T2<:Number, TE2<:ETYPE} +# function vcat(SYS1::DescriptorStateSpace{T1, TE1}, SYS2::DescriptorStateSpace{T2, TE2}) where {T1<:Number, TE1<:ETYPE, T2<:Number, TE2<:ETYPE} +# #function vcat(SYS1 :: DescriptorStateSpace, SYS2 :: DescriptorStateSpace) +# nu = SYS1.nu +# nu == size(SYS2, 2) || error("The systems must have the same input dimension") +# #T = promote_type(eltype(SYS1), eltype(SYS2)) +# T = promote_type(T1, T2) +# #TE = promote_type(TE1, TE2) +# TE = promote_Etype(T, TE1, TE2) +# Ts = promote_Ts(SYS1.Ts, SYS2.Ts) + +# A = blockdiag(T.(SYS1.A), T.(SYS2.A)) + +# local E::TE +# if SYS1.E === I && SYS2.E === I +# #if SYS1.E == I && SYS2.E == I +# E = I +# elseif SYS1.E == I || SYS2.E == I +# blockdims = [size(SYS1.A,1), size(SYS2.A,1)] +# E = sblockdiag(blockdims, SYS1.E == I ? SYS1.E : T.(SYS1.E), SYS2.E == I ? SYS2.E : T.(SYS2.E)) +# else +# E = blockdiag(T.(SYS1.E), T.(SYS2.E)) +# end +# C = blockdiag(T.(SYS1.C), T.(SYS2.C)) +# B = [ T.(SYS1.B); T.(SYS2.B)] +# D = [ T.(SYS1.D); T.(SYS2.D)] +# return DescriptorStateSpace{T}(A, E, B, C, D, Ts) +# end + +function vcat(SYS1::DescriptorStateSpace{T1}, SYS2::DescriptorStateSpace{T2}) where {T1<:Number, T2<:Number} #function vcat(SYS1 :: DescriptorStateSpace, SYS2 :: DescriptorStateSpace) nu = SYS1.nu nu == size(SYS2, 2) || error("The systems must have the same input dimension") #T = promote_type(eltype(SYS1), eltype(SYS2)) T = promote_type(T1, T2) #TE = promote_type(TE1, TE2) + TE1 = typeof(SYS1.E) + TE2= typeof(SYS2.E) TE = promote_Etype(T, TE1, TE2) Ts = promote_Ts(SYS1.Ts, SYS2.Ts) diff --git a/src/dss.jl b/src/dss.jl index 30b9fe1..995a55b 100644 --- a/src/dss.jl +++ b/src/dss.jl @@ -58,7 +58,7 @@ function dss(A::AbstractNumOrArray, E::Union{AbstractNumOrArray,UniformScaling}, check_reg && E != I && !MatrixPencils.isregular(to_matrix(T,A), to_matrix(T,E), atol1 = atol1, atol2 = atol2, rtol = rtol ) && error("The pencil A-λE is not regular") p = typeof(C) <: Union{Vector,Number} ? (size(A,1) <= 1 ? size(C,1) : 1) : size(C,1) - m = typeof(B) <: Union{Vector,Number} ? 1 : size(B,2) + m = typeof(B) <: Union{Vector,Number} ? 1 : size(B,2) return DescriptorStateSpace{T}(to_matrix(T,A), E == I ? I : to_matrix(T,E), 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 @@ -79,18 +79,18 @@ function dss(D::AbstractNumOrArray; Ts::Real = 0) T = eltype(D) return DescriptorStateSpace{T}(zeros(T,0,0),I,zeros(T,0,m),zeros(T,p,0),to_matrix(T,D),Ts) end -function dss(;A::Union{AbstractNumOrArray} = zeros(Bool,0,0), E::Union{AbstractNumOrArray,UniformScaling} = I, B::Union{AbstractNumOrArray,Missing} = missing, +function dss(;A::AbstractNumOrArray = zeros(Bool,0,0), E::Union{AbstractNumOrArray,UniformScaling} = I, B::Union{AbstractNumOrArray,Missing} = missing, C::Union{AbstractNumOrArray,Missing} = missing, D::Union{AbstractNumOrArray,Missing} = missing, Ts::Real = 0, check_reg::Bool = false, atol::Real = zero(real(eltype(A))), atol1::Real = atol, atol2::Real = atol, rtol::Real = (typeof(A) <: Number ? 1 : min(size(A)...))*eps(real(float(one(eltype(A)))))*iszero(min(atol1,atol2))) T = Bool - T = promote_type(T, ismissing(A) ? T : eltype(A), E == I ? Bool : eltype(E), ismissing(B) ? T : eltype(B), + T = promote_type(T, eltype(A), E == I ? Bool : eltype(E), ismissing(B) ? T : eltype(B), ismissing(C) ? T : eltype(C), ismissing(D) ? T : eltype(D)) A = to_matrix(T,A) n = size(A,1) E == I || (E = to_matrix(T,E)) - check_reg && E != I && !isregular(A, E, atol1 = atol1, atol2 = atol2, rtol = rtol ) && + check_reg && E != I && !MatrixPencils.isregular(A, E, atol1 = atol1, atol2 = atol2, rtol = rtol ) && error("The pencil A-λE is not regular") ismissing(D) || (D1 = to_matrix(T,D)) diff --git a/src/dsutils.jl b/src/dsutils.jl index 9c61cf3..1f905f1 100644 --- a/src/dsutils.jl +++ b/src/dsutils.jl @@ -8,6 +8,11 @@ to_matrix(T, A::AbstractMatrix, wide::Bool = false) = AbstractMatrix{T}(A) # Fa 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)) +to_matrix(T, A::Diagonal, wide::Bool = false) = Matrix{T}(A) # Fallback +to_matrix(T, A::UpperTriangular, wide::Bool = false) = Matrix{T}(A) # Fallback +to_matrix(T, A::LowerTriangular, wide::Bool = false) = Matrix{T}(A) # Fallback +to_matrix(T, A::UpperHessenberg, wide::Bool = false) = Matrix{T}(A) # Fallback + # eye(n) = Matrix{Bool}(I, n, n) diff --git a/src/types/DescriptorStateSpace.jl b/src/types/DescriptorStateSpace.jl index 9d99a85..eae2f72 100644 --- a/src/types/DescriptorStateSpace.jl +++ b/src/types/DescriptorStateSpace.jl @@ -65,20 +65,34 @@ The dimensions `nx`, `ny` and `nu` can be obtained as `SYS.nx`, `SYS.ny` and `SY # new{T, typeof(E)}(A, E, B, C, D, Float64(Ts)) # end # end -struct DescriptorStateSpace{T, ET <: ETYPE{T}, MT<:AbstractMatrix{T}} <: AbstractDescriptorStateSpace +# struct DescriptorStateSpace{T, ET <: ETYPE{T}, MT<:AbstractMatrix{T}} <: AbstractDescriptorStateSpace +# A::MT +# E::ET +# B::MT +# C::MT +# D::MT +# Ts::Float64 +# 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), MT}(A, E, B, C, D, Float64(Ts)) +# end +# end +struct DescriptorStateSpace{T, ET <: Union{AbstractMatrix{T},UniformScaling{Bool}}, MT<:AbstractMatrix{T}} <: AbstractDescriptorStateSpace A::MT E::ET B::MT C::MT D::MT Ts::Float64 - function DescriptorStateSpace{T}(A::MT, E::ETYPE{T}, + function DescriptorStateSpace{T}(A::MT, E::Union{MT,UniformScaling{Bool}}, 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), 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 diff --git a/test/test_dss.jl b/test/test_dss.jl index 9799f32..8c6f25a 100644 --- a/test/test_dss.jl +++ b/test/test_dss.jl @@ -3,6 +3,8 @@ module Test_dss using DescriptorSystems using LinearAlgebra using Polynomials +using SparseArrays +using Measurements using Test println("Test_dss") @@ -225,6 +227,30 @@ catch @test false end +# some tests involving the new DescriptorStateSpace structure + +t = rand(3,3) +sys = dss(UpperHessenberg(t), UpperTriangular(t),LowerTriangular(t),Diagonal(t),0) +@test istriu(sys.A,-1) && istriu(sys.E) && istril(sys.B) && isdiag(sys.C) && iszero(sys.D) + +ssys = dss(sparse(sys.A),sparse(sys.E),sparse(sys.B),sparse(sys.C),sparse(sys.D)) +@test istriu(ssys.A,-1) && istriu(ssys.E) && istril(ssys.B) && isdiag(ssys.C) && iszero(ssys.D) + +ρ1 = measurement(0, 0.25); ρ2 = measurement(0, 0.25); +# build uncertain state matrix A(p) +A = [-.8 0 0;0 -.5*(1+ρ1) .6*(1+ρ2); 0 -0.6*(1+ρ2) -0.5*(1+ρ1)]; +B = [1 1;1 0;0 1]; +C = [0 1 1; 1 1 0]; D = zeros(2,2); +# build an uncertain system +try + usys = dss(A,B,C,D) + @test true +catch + @test false +end + + + end