diff --git a/src/analysis.jl b/src/analysis.jl index 559d87e1a..688454371 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -366,6 +366,9 @@ https://arxiv.org/pdf/1805.10312.pdf """ function relative_gain_array(A::AbstractMatrix; tol = 1e-15) m, n = size(A) + if m == n && LinearAlgebra.det(A) != 0 + return A .* inv(transpose(A)) + end L = zeros(m, n) M = ones(m, n) S = sign.(A) diff --git a/src/types/StateSpace.jl b/src/types/StateSpace.jl index 219f3288f..913ab7575 100644 --- a/src/types/StateSpace.jl +++ b/src/types/StateSpace.jl @@ -140,7 +140,13 @@ struct HeteroStateSpace{TE <: TimeEvolution, AT<:AbstractMatrix,BT<:AbstractVecO # Explicit constructor function HeteroStateSpace{TE,AT,BT,CT,DT}(A, B, C, D, timeevol) where {TE,AT,BT,CT,DT} state_space_validation(A,B,C,D) - new{TE,AT,BT,CT,DT}(AT(A), BT(B), CT(C), DT(D), TE(timeevol)) + new{TE,AT,BT,CT,DT}( + A isa AT ? A : AT(A), + B isa BT ? B : BT(B), + C isa CT ? C : CT(C), + D isa DT ? D : DT(D), + TE(timeevol) + ) end # Base constructor function HeteroStateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, timeevol::TimeEvolution) @@ -218,6 +224,7 @@ end ## ADDITION ## Base.zero(sys::AbstractStateSpace) = ss(zero(sys.D), sys.timeevol) +Base.zero(::Type{StateSpace{Continuous, F}}) where {F} = ss([zero(F)], Continuous()) # Cannot make a zero of discrete system since sample time is not stored in type. function +(s1::StateSpace{TE,T}, s2::StateSpace{TE,T}) where {TE,T} #Ensure systems have same dimensions @@ -266,7 +273,13 @@ end function *(sys1::StateSpace{TE,T}, sys2::StateSpace{TE,T}) where {TE,T} #Check dimension alignment #Note: sys1*sys2 = y <- sys1 <- sys2 <- u - if sys1.nu != sys2.ny + if xor(issiso(sys1), issiso(sys2)) + if issiso(sys1) + sys1 = append(fill(sys1, sys2.ny)...) + else + sys2 = append(fill(sys2, sys1.nu)...) + end + elseif sys1.nu != sys2.ny error("sys1*sys2: sys1 must have same number of inputs as sys2 has outputs") end timeevol = common_timeevol(sys1,sys2) @@ -282,7 +295,13 @@ end function *(sys1::HeteroStateSpace, sys2::HeteroStateSpace) #Check dimension alignment #Note: sys1*sys2 = y <- sys1 <- sys2 <- u - if sys1.nu != sys2.ny + if xor(issiso(sys1), issiso(sys2)) + if issiso(sys1) + sys1 = append(fill(sys1, sys2.ny)...) + else + sys2 = append(fill(sys2, sys1.nu)...) + end + elseif sys1.nu != sys2.ny error("sys1*sys2: sys1 must have same number of inputs as sys2 has outputs") end timeevol = common_timeevol(sys1,sys2) diff --git a/src/types/promotion.jl b/src/types/promotion.jl index 421cc0cec..bb4bf7d60 100644 --- a/src/types/promotion.jl +++ b/src/types/promotion.jl @@ -67,5 +67,10 @@ Base.promote_rule(::Type{TransferFunction{TE, SisoRational{T1}}}, ::Type{MT}) wh Base.promote_rule(::Type{StateSpace{TE, T1}}, ::Type{MT}) where {TE, T1, MT<:AbstractMatrix} = StateSpace{TE, promote_type(T1,eltype(MT))} +function Base.promote_rule(::Type{HeteroStateSpace{TE, T1, T2, T3, T4}}, ::Type{MT}) where {TE, T1, T2, T3, T4, MT<:AbstractMatrix} + T = promote_type(eltype(T2),eltype(T3),eltype(T4),eltype(MT)) + StateSpace{TE, T} # We can't know how MT will interact with the hss, so we fall back to standard ss +end + Base.promote_rule(::Type{DelayLtiSystem{T1,S}}, ::Type{MT}) where {T1, S, MT<:AbstractMatrix} = DelayLtiSystem{promote_type(T1, eltype(MT)),S} diff --git a/test/test_analysis.jl b/test/test_analysis.jl index 03820c368..84cc072d6 100644 --- a/test/test_analysis.jl +++ b/test/test_analysis.jl @@ -211,25 +211,21 @@ marginplot(P, w) # RGA a = 10 -P = [ - tf([1,-a^2], [1, 0, a^2]) tf([a, a], [1, 0, a^2]) - -tf([a, a], [1, 0, a^2]) tf([1,-a^2], [1, 0, a^2]) - ] -P = minreal(ss(P)) +P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) w = exp10.(LinRange(-1, 2, 1000)) rgaplot(P, w) rgaplot([P, 2P], w) R = relative_gain_array(P, w) -@test maximum(abs, R) > 1e14 +@test maximum(abs, R) > 50 # Inf/NaN at w=10 @test minimum(abs, R) ≈ 1e-2 atol=1e-4 -R = relative_gain_array(P, 10) +R = relative_gain_array(P, 9.99) @test size(R) == size(P) -@test R[1,1] ≈ R[2,2] -@test R[2,1] ≈ R[1,2] -@test R[1,1] ≈ -R[1,2] +@test R[1,1] ≈ R[2,2] rtol=0.01 +@test R[2,1] ≈ R[1,2] rtol=0.01 +@test R[1,1] ≈ -R[1,2] rtol=0.01 # tests from https://arxiv.org/pdf/1805.10312.pdf diff --git a/test/test_statespace.jl b/test/test_statespace.jl index d9caf7019..f759d238b 100644 --- a/test/test_statespace.jl +++ b/test/test_statespace.jl @@ -61,7 +61,6 @@ @test D_111 - D_211 == SS([-0.5 0 0; 0 0.2 -0.8; 0 -0.8 0.07],[2; 1; 2], [3 -1 -0],[0], 0.005) - # Multiplication @test C_111 * C_221 == SS([-5 2 0; 0 -5 -3; 0 2 -9], [0 0; 1 0; 0 2],[3 0 0],[0 0]) @@ -70,6 +69,9 @@ @test 4*C_222 == SS([-5 -3; 2 -9],[1 0; 0 2],[4 0; 0 4],[0 0; 0 0]) @test D_111 * D_221 == SS([-0.5 2 0; 0 0.2 -0.8; 0 -0.8 0.07], [0 0; 1 0; 0 2],[3 0 0],[0 0],0.005) + @test C_111 * I(2) == I(2) * C_111 == SS(diagm([a_1; a_1]), 2*I(2), 3*I(2), 0*I(2)) + @test minreal(C_111*C_222_d - C_222_d*C_111, atol=1e-3) == ss(0*I(2)) # scalar times MIMO + @test C_111*C_222 == ss([-5 0 2 0; 0 -5 0 2; 0 0 -5 -3; 0 0 2 -9], [0 0; 0 0; 1 0; 0 2], [3 0 0 0; 0 3 0 0], 0) # Division @test 1/C_222_d == SS([-6 -3; 2 -11],[1 0; 0 2],[-1 0; -0 -1],[1 -0; 0 1]) @@ -121,7 +123,6 @@ # Errors @test_throws ErrorException C_111 + C_222 # Dimension mismatch @test_throws ErrorException C_111 - C_222 # Dimension mismatch - @test_throws ErrorException C_111 * C_222 # Dimension mismatch @test_throws ErrorException D_111 + C_111 # Sampling time mismatch @test_throws ErrorException D_111 - C_111 # Sampling time mismatch @test_throws ErrorException D_111 * C_111 # Sampling time mismatch