usability improvements for multiplication with scalars

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)

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)

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}

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

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