Removed the matrix type parameter from StateSpace, closes #413 (#437)

* Removed the matrix type parameter from StateSpace, closes #413

* Added some tests for conversions and promotion

* Improved TransferFunction -> StateSpace conversion

Author: olof3

docs/src/man/creating_systems.md

@@ -91,7 +91,7 @@ ControlSystems._string_mat_with_headers(a::SizedArray) = ControlSystems._string_
 Notice the different matrix types used
 ```jldoctest HSS
 julia> sys = ss([-5 0; 0 -5],[2; 2],[3 3],[0])
-StateSpace{Continuous,Int64,Array{Int64,2}}
+StateSpace{Continuous,Int64}
 A =
  -5   0
   0  -5

src/types/DelayLtiSystem.jl

@@ -4,7 +4,7 @@
 Represents an LTISystem with internal time-delay. See `?delay` for a convenience constructor.
 """
 struct DelayLtiSystem{T,S<:Real} <: LTISystem
-    P::PartionedStateSpace{StateSpace{Continuous,T,Matrix{T}}}
+    P::PartionedStateSpace{StateSpace{Continuous,T}}
     Tau::Vector{S} # The length of the vector tau implicitly defines the partitionging of P
 
     # function DelayLtiSystem(P::StateSpace{Continuous,T, MT}, Tau::Vector{T})
@@ -33,19 +33,19 @@ function DelayLtiSystem{T,S}(sys::StateSpace, Tau::AbstractVector{S} = Float64[]
         throw(ArgumentError("The delay vector of length $length(Tau) is too long."))
     end
 
-    psys = PartionedStateSpace{StateSpace{Continuous,T,Matrix{T}}}(sys, nu, ny)
+    psys = PartionedStateSpace{StateSpace{Continuous,T}}(sys, nu, ny)
     DelayLtiSystem{T,S}(psys, Tau)
 end
 # For converting DelayLtiSystem{T,S} to different T
-DelayLtiSystem{T}(sys::DelayLtiSystem) where {T} = DelayLtiSystem{T}(PartionedStateSpace{StateSpace{Continuous,T,Matrix{T}}}(sys.P), Float64[])
+DelayLtiSystem{T}(sys::DelayLtiSystem) where {T} = DelayLtiSystem{T}(PartionedStateSpace{StateSpace{Continuous,T}}(sys.P), Float64[])
 DelayLtiSystem{T}(sys::StateSpace) where {T} = DelayLtiSystem{T, Float64}(sys, Float64[])
 
 # From StateSpace, infer type
-DelayLtiSystem(sys::StateSpace{Continuous,T,MT}, Tau::Vector{S}) where {T, MT,S} = DelayLtiSystem{T,S}(sys, Tau)
-DelayLtiSystem(sys::StateSpace{Continuous,T,MT}) where {T, MT} = DelayLtiSystem{T,T}(sys, T[])
+DelayLtiSystem(sys::StateSpace{Continuous,T}, Tau::Vector{S}) where {T,S} = DelayLtiSystem{T,S}(sys, Tau)
+DelayLtiSystem(sys::StateSpace{Continuous,T}) where {T} = DelayLtiSystem{T,T}(sys, T[])
 
 # From TransferFunction, infer type TODO Use proper constructor instead of convert here when defined
-DelayLtiSystem(sys::TransferFunction{TE,S}) where {TE,T,S<:SisoTf{T}} = DelayLtiSystem{T}(convert(StateSpace{Continuous,T, Matrix{T}}, sys))
+DelayLtiSystem(sys::TransferFunction{TE,S}) where {TE,T,S<:SisoTf{T}} = DelayLtiSystem{T}(convert(StateSpace{Continuous,T}, sys))
 
 
 # TODO: Think through these promotions and conversions
@@ -69,7 +69,7 @@ end
 Base.convert(::Type{<:DelayLtiSystem}, sys::TransferFunction)  = DelayLtiSystem(sys)
 # Catch convertsion between T
 Base.convert(::Type{V}, sys::DelayLtiSystem)  where {T, V<:DelayLtiSystem{T}} =
-    sys isa V ? sys : V(StateSpace{Continuous,T,Matrix{T}}(sys.P.P), sys.Tau)
+    sys isa V ? sys : V(StateSpace{Continuous,T}(sys.P.P), sys.Tau)
 
 
 
@@ -109,6 +109,13 @@ function /(anything, sys::DelayLtiSystem)
     /(anything, sys.P.P) # If all delays are zero, invert the inner system
 end
 
+
+# Test equality (of realizations)
+function ==(sys1::DelayLtiSystem, sys2::DelayLtiSystem)
+    all(getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(DelayLtiSystem))
+end
+
+
 ninputs(sys::DelayLtiSystem) = size(sys.P.P, 2) - length(sys.Tau)
 noutputs(sys::DelayLtiSystem) = size(sys.P.P, 1) - length(sys.Tau)
 nstates(sys::DelayLtiSystem) = nstates(sys.P.P)

src/types/PartionedStateSpace.jl

@@ -220,3 +220,8 @@ end
 function Base.convert(::Type{<:PartionedStateSpace}, sys::T) where T<: StateSpace
     PartionedStateSpace(sys,sys.nu,sys.ny)
 end
+
+# Test equality (of realizations)
+function ==(sys1::PartionedStateSpace, sys2::PartionedStateSpace)
+    all(getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(PartionedStateSpace))
+end

src/types/StateSpace.jl

@@ -24,28 +24,28 @@ end
 
 abstract type AbstractStateSpace{TE<:TimeEvolution} <: LTISystem end
 
-struct StateSpace{TE, T, MT<:AbstractMatrix{T}} <: AbstractStateSpace{TE}
-    A::MT
-    B::MT
-    C::MT
-    D::MT
+struct StateSpace{TE, T} <: AbstractStateSpace{TE}
+    A::Matrix{T}
+    B::Matrix{T}
+    C::Matrix{T}
+    D::Matrix{T}
     timeevol::TE
-    function StateSpace{TE, T, MT}(A::MT, B::MT, C::MT, D::MT, timeevol::TE) where {TE<:TimeEvolution, T, MT <: AbstractMatrix{T}}
+    function StateSpace{TE, T}(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, timeevol::TE) where {TE<:TimeEvolution, T}
         state_space_validation(A,B,C,D)
-        new{TE, T, MT}(A, B, C, D, timeevol)
+        new{TE, T}(A, B, C, D, timeevol)
     end
-    function StateSpace(A::MT, B::MT, C::MT, D::MT, timeevol::TE) where {TE<:TimeEvolution, T, MT <: AbstractMatrix{T}}
+    function StateSpace(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, timeevol::TE) where {TE<:TimeEvolution, T}
         state_space_validation(A,B,C,D)
-        new{TE, T, MT}(A, B, C, D, timeevol)
+        new{TE, T}(A, B, C, D, timeevol)
     end
 end
 # Constructor for Discrete system
-function StateSpace(A::MT, B::MT, C::MT, D::MT, Ts::Number) where {T, MT <: AbstractMatrix{T}}
+function StateSpace(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, Ts::Number) where {T}
     state_space_validation(A,B,C,D)
     StateSpace(A, B, C, D, Discrete(Ts))
 end
 # Constructor for Continuous system
-function StateSpace(A::MT, B::MT, C::MT, D::MT) where {T, MT <: AbstractMatrix{T}}
+function StateSpace(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}) where {T}
     state_space_validation(A,B,C,D)
     StateSpace(A, B, C, D, Continuous())
 end
@@ -72,19 +72,19 @@ end
 const AbstractNumOrArray = Union{Number, AbstractVecOrMat}
 
 # Explicit construction with different types
-function StateSpace{TE,T,MT}(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, timeevol::TimeEvolution) where {TE, T, MT <: AbstractMatrix{T}}
+function StateSpace{TE,T}(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, timeevol::TimeEvolution) where {TE, T}
     D = fix_D_matrix(T, B, C, D)
-    return StateSpace{TE, T,Matrix{T}}(MT(to_matrix(T, A)), MT(to_matrix(T, B)), MT(to_matrix(T, C)), MT(D), TE(timeevol))
+    return StateSpace{TE, T}(to_matrix(T, A), to_matrix(T, B), to_matrix(T, C), to_matrix(T, D), TE(timeevol))
 end
 
 # Explicit conversion
-function StateSpace{TE,T,MT}(sys::StateSpace) where {TE, T, MT <: AbstractMatrix{T}}
-    StateSpace{TE,T,MT}(sys.A,sys.B,sys.C,sys.D,sys.timeevol)
+function StateSpace{TE,T}(sys::StateSpace) where {TE, T}
+    StateSpace{TE,T}(sys.A,sys.B,sys.C,sys.D,sys.timeevol)
 end
 
 function StateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, timeevol::TimeEvolution)
     A, B, C, D, T = to_similar_matrices(A,B,C,D)
-    return StateSpace{typeof(timeevol),T,Matrix{T}}(A, B, C, D, timeevol)
+    return StateSpace{typeof(timeevol),T}(A, B, C, D, timeevol)
 end
 # General Discrete constructor
 StateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, Ts::Number) =
@@ -115,13 +115,13 @@ StateSpace(sys::LTISystem) = convert(StateSpace, sys)
 """
     `sys = ss(A, B, C, D [,Ts])`
 
-Create a state-space model `sys::StateSpace{TE, T, MT<:AbstractMatrix{T}}`
-where `MT` is the type of matrixes `A,B,C,D`, `T` the element type and TE is Continuous or Discrete.
+Create a state-space model `sys::StateSpace{TE, T}`
+with matrix element type `T` and TE is `Continuous` or `<:Discrete`.
 
 This is a continuous-time model if `Ts` is omitted.
-Otherwise, this is a discrete-time model with sampling period Ts.
+Otherwise, this is a discrete-time model with sampling period `Ts`.
 
-`sys = ss(D [, Ts])` specifies a static gain matrix D.
+`sys = ss(D [, Ts])` specifies a static gain matrix `D`.
 """
 ss(args...;kwargs...) = StateSpace(args...;kwargs...)
 
@@ -214,7 +214,7 @@ function isapprox(sys1::ST1, sys2::ST2; kwargs...) where {ST1<:AbstractStateSpac
 end
 
 ## ADDITION ##
-function +(s1::StateSpace{TE,T,MT}, s2::StateSpace{TE,T,MT}) where {TE,T, MT}
+function +(s1::StateSpace{TE,T}, s2::StateSpace{TE,T}) where {TE,T}
     #Ensure systems have same dimensions
     if size(s1) != size(s2)
         error("Systems have different shapes.")
@@ -227,7 +227,7 @@ function +(s1::StateSpace{TE,T,MT}, s2::StateSpace{TE,T,MT}) where {TE,T, MT}
     C = [s1.C s2.C;]
     D = [s1.D + s2.D;]
 
-    return StateSpace{TE,T,MT}(A, B, C, D, timeevol)
+    return StateSpace{TE,T}(A, B, C, D, timeevol)
 end
 
 function +(s1::HeteroStateSpace, s2::HeteroStateSpace)
@@ -258,7 +258,7 @@ end
 -(sys::ST) where ST <: AbstractStateSpace = ST(sys.A, sys.B, -sys.C, -sys.D, sys.timeevol)
 
 ## MULTIPLICATION ##
-function *(sys1::StateSpace{TE,T,MT}, sys2::StateSpace{TE,T,MT}) where {TE,T, MT}
+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
@@ -271,7 +271,7 @@ function *(sys1::StateSpace{TE,T,MT}, sys2::StateSpace{TE,T,MT}) where {TE,T, MT
     B = [sys1.B*sys2.D ; sys2.B]
     C = [sys1.C   sys1.D*sys2.C;]
     D = [sys1.D*sys2.D;]
-    return StateSpace{TE,T,MT}(A, B, C, D, timeevol)
+    return StateSpace{TE,T}(A, B, C, D, timeevol)
 end
 
 function *(sys1::HeteroStateSpace, sys2::HeteroStateSpace)

src/types/conversion.jl

@@ -19,37 +19,34 @@
 # Base.convert(::Type{<:TransferFunction{<:SisoRational}}, b::Number) = tf(b)
 # Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, b::Number) = zpk(b)
 #
-Base.convert(::Type{TransferFunction{TE,SisoZpk{T1, TR1}}}, b::AbstractMatrix{T2}) where {TE, T1, TR1, T2<:Number} =
-    zpk(T1.(b), undef_sampletime(TE))
-Base.convert(::Type{TransferFunction{TE,SisoRational{T1}}}, b::AbstractMatrix{T2}) where {TE, T1, T2<:Number} =
-    tf(T1.(b), undef_sampletime(TE))
+Base.convert(::Type{TransferFunction{TE,SisoZpk{T1, TR1}}}, D::AbstractMatrix{T2}) where {TE, T1, TR1, T2<:Number} =
+    zpk(convert(Matrix{T1}, D), undef_sampletime(TE))
+Base.convert(::Type{TransferFunction{TE,SisoRational{T1}}}, D::AbstractMatrix{T2}) where {TE, T1, T2<:Number} =
+    tf(Matrix{T1}(D), undef_sampletime(TE))
 
-function convert(::Type{StateSpace{TE,T,MT}}, D::AbstractMatrix{<:Number}) where {TE,T, MT}
+function convert(::Type{StateSpace{TE,T}}, D::AbstractMatrix{<:Number}) where {TE,T}
     (ny, nu) = size(D)
-    A = MT(fill(zero(T), (0,0)))
-    B = MT(fill(zero(T), (0,nu)))
-    C = MT(fill(zero(T), (ny,0)))
-    D = convert(MT, D)
-    return StateSpace{TE,T,MT}(A,B,C,D,undef_sampletime(TE))
+    A = fill(zero(T), (0,0))
+    B = fill(zero(T), (0,nu))
+    C = fill(zero(T), (ny,0))
+    D = convert(Matrix{T}, D)
+    return StateSpace{TE,T}(A,B,C,D,undef_sampletime(TE))
 end
 
 # TODO We still dont have TransferFunction{..}(number/array) constructors
-Base.convert(::Type{TransferFunction{TE,SisoRational{T}}}, b::Number) where {TE, T} =
-    tf(T(b), undef_sampletime(TE))
-Base.convert(::Type{TransferFunction{TE,SisoZpk{T,TR}}}, b::Number) where {TE, T, TR} =
-    zpk(T(b), undef_sampletime(TE))
-Base.convert(::Type{StateSpace{TE,T,MT}}, b::Number) where {TE, T, MT} =
-    convert(StateSpace{TE,T,MT}, fill(b, (1,1)))
+Base.convert(::Type{TransferFunction{TE,SisoRational{T}}}, d::Number) where {TE, T} =
+    tf(T(d), undef_sampletime(TE))
+
+Base.convert(::Type{TransferFunction{TE,SisoZpk{T,TR}}}, d::Number) where {TE, T, TR} =
+    zpk(T(d), undef_sampletime(TE))
+
+Base.convert(::Type{StateSpace{TE,T}}, d::Number) where {TE, T} =
+    convert(StateSpace{TE,T}, fill(d, (1,1)))
 #
 # Base.convert(::Type{TransferFunction{Continuous,<:SisoRational{T}}}, b::Number) where {T} = tf(T(b), Continuous())
 # Base.convert(::Type{TransferFunction{Continuous,<:SisoZpk{T, TR}}}, b::Number) where {T, TR} = zpk(T(b), Continuous())
 # Base.convert(::Type{StateSpace{Continuous,T, MT}}, b::Number) where {T, MT} = convert(StateSpace{Continuous,T, MT}, fill(b, (1,1)))
 
-#
-# Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, s::TransferFunction) = zpk(s)
-# Base.convert(::Type{<:TransferFunction{<:SisoRational}}, s::TransferFunction) = tf(s)
-
-#
 # function Base.convert{T<:Real,S<:TransferFunction}(::Type{S}, b::VecOrMat{T})
 #     r = Matrix{S}(size(b,2),1)
 #     for j=1:size(b,2)
@@ -64,16 +61,16 @@ function convert(::Type{TransferFunction{TE,S}}, G::TransferFunction) where {TE,
     return TransferFunction{TE,eltype(Gnew_matrix)}(Gnew_matrix, TE(G.timeevol))
 end
 
-function convert(::Type{S}, sys::AbstractStateSpace) where {T, MT, TE, S <:StateSpace{TE,T,MT}}
-    if sys isa S
+function convert(::Type{StateSpace{TE,T}}, sys::AbstractStateSpace) where {TE, T}
+    if sys isa StateSpace{TE, T}
         return sys
     else
-        return StateSpace{TE, T,MT}(convert(MT, sys.A), convert(MT, sys.B), convert(MT, sys.C), convert(MT, sys.D), TE(sys.timeevol))
+        return StateSpace{TE, T}(convert(Matrix{T}, sys.A), convert(Matrix{T}, sys.B), convert(Matrix{T}, sys.C), convert(Matrix{T}, sys.D), TE(sys.timeevol))
     end
 end
 
 # TODO Maybe add convert on matrices
-Base.convert(::Type{HeteroStateSpace{TE1,AT,BT,CT,DT}}, s::StateSpace{TE2,T,MT}) where {TE1,TE2,T,MT,AT,BT,CT,DT} =
+Base.convert(::Type{HeteroStateSpace{TE1,AT,BT,CT,DT}}, s::StateSpace{TE2,T}) where {TE1,TE2,T,AT,BT,CT,DT} =
     HeteroStateSpace{TE1,AT,BT,CT,DT}(s.A,s.B,s.C,s.D,TE1(s.timeevol))
 
 Base.convert(::Type{HeteroStateSpace}, s::StateSpace) = HeteroStateSpace(s)
@@ -84,45 +81,41 @@ Base.convert(::Type{StateSpace}, s::HeteroStateSpace{Continuous}) = StateSpace(s
 function Base.convert(::Type{StateSpace}, G::TransferFunction{TE,<:SisoTf{T0}}) where {TE,T0<:Number}
 
     T = Base.promote_op(/,T0,T0)
-    convert(StateSpace{TE,T,Matrix{T}}, G)
+    convert(StateSpace{TE,T}, G)
 end
 
-
-function Base.convert(::Type{StateSpace{TE,T,MT}}, G::TransferFunction) where {TE,T<:Number, MT<:AbstractArray{T}}
+function Base.convert(::Type{StateSpace{TE,T}}, G::TransferFunction) where {TE,T<:Number}
     if !isproper(G)
         error("System is improper, a state-space representation is impossible")
     end
 
-    # TODO : These are added due to scoped for blocks, but is a hack. This
-    # could be much cleaner.
-    #T = Base.promote_op(/, T0, T0)
-
-    Ac = Bc = Cc = Dc = A = B = C = D = Array{T}(undef, 0, 0)
-    for i=1:ninputs(G)
-        for j=1:noutputs(G)
-            a, b, c, d = siso_tf_to_ss(T, G.matrix[j, i])
-            if j > 1
-                # vcat
-                Ac = blockdiag(Ac, a)
-                Bc = vcat(Bc, b)
-                Cc = blockdiag(Cc, c)
-                Dc = vcat(Dc, d)
-            else
-                Ac, Bc, Cc, Dc = a, b, c, d
-            end
-        end
-        if i > 1
-            # hcat
-            A = blockdiag(A, Ac)
-            B = blockdiag(B, Bc)
-            C = hcat(C, Cc)
-            D = hcat(D, Dc)
-        else
-            A, B, C, D = Ac, Bc, Cc, Dc
+    ny, nu = size(G)
+
+    # A, B, C, D matrices for each element of the transfer function matrix
+    abcd_vec = [siso_tf_to_ss(T, g) for g in G.matrix[:]]
+
+    # Number of states for each transfer function element realization
+    nvec = [size(abcd[1], 1) for abcd in abcd_vec]
+    ntot = sum(nvec)
+
+    A = zeros(T, (ntot, ntot))
+    B = zeros(T, (ntot, nu))
+    C = zeros(T, (ny, ntot))
+    D = zeros(T, (ny, nu))
+
+    inds = -1:0
+    for j=1:nu
+        for i=1:ny
+            k = (j-1)*ny + i
+
+            # states correpsonding to the transfer function element (i,j)
+            inds = (inds.stop+1):(inds.stop+nvec[k])
+
+            A[inds,inds], B[inds,j:j], C[i:i,inds], D[i:i,j:j] = abcd_vec[k]
         end
     end
     # A, B, C = balance_statespace(A, B, C)[1:3] NOTE: Use balance?
-    return StateSpace{TE,T,MT}(A, B, C, D, TE(G.timeevol))
+    return StateSpace{TE,T}(A, B, C, D, TE(G.timeevol))
 end
 
 siso_tf_to_ss(T::Type, f::SisoTf) = siso_tf_to_ss(T, convert(SisoRational, f))
@@ -131,8 +124,8 @@ siso_tf_to_ss(T::Type, f::SisoTf) = siso_tf_to_ss(T, convert(SisoRational, f))
 function siso_tf_to_ss(T::Type, f::SisoRational)
 
     num0, den0 = numvec(f), denvec(f)
-    # Normalize the numerator and denominator,
-    # To allow realization of transfer functions that are proper, but now strictly proper
+    # Normalize the numerator and denominator to allow realization of transfer functions
+    # that are proper, but not strictly proper
     num = num0 / den0[1]
     den = den0 / den0[1]
 
@@ -141,22 +134,22 @@ function siso_tf_to_ss(T::Type, f::SisoRational)
     # Get numerator coefficient of the same order as the denominator
     bN = length(num) == N+1 ? num[1] : 0
 
-    if N==0 #|| num == zero(Polynomial{T})
-        A = fill(zero(T), 0, 0)
-        B = fill(zero(T), 0, 1)
-        C = fill(zero(T), 1, 0)
+    if N == 0 #|| num == zero(Polynomial{T})
+        A = zeros(T, (0, 0))
+        B = zeros(T, (0, 1))
+        C = zeros(T, (1, 0))
     else
-        A = diagm(1 => fill(one(T), N-1))
+        A = diagm(1 => ones(T, N-1))
         A[end, :] .= -reverse(den)[1:end-1]
 
-        B = fill(zero(T), N, 1)
+        B = zeros(T, (N, 1))
         B[end] = one(T)
 
-        C = fill(zero(T), 1, N)
+        C = zeros(T, (1, N))
         C[1:min(N, length(num))] = reverse(num)[1:min(N, length(num))]
         C[:] -= bN * reverse(den)[1:end-1] # Can index into polynomials at greater inddices than their length
     end
-    D = fill(bN, 1, 1)
+    D = fill(bN, (1, 1))
 
     return A, B, C, D
 end