Squash olof3:sampletimes2 into master. Implements type stable sample times

Author: mfalt

.travis.yml

@@ -2,9 +2,11 @@ language: julia
 julia:
   - 1.0
   - 1.3
+  - 1.5
   - nightly
 matrix:
   allow_failures:
+  - julia: 1.3
   - julia: nightly
   fast_finish: true
 notifications:

Project.toml

@@ -27,5 +27,5 @@ IterTools = "1.0"
 LaTeXStrings = "1.0"
 OrdinaryDiffEq = "5.2"
 Plots = "0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 1.0"
-Polynomials = "0.7, 0.8, 1.0"
+Polynomials = "1.0"
 julia = "1.0"

README.md

@@ -1,6 +1,6 @@
 # ControlSystems.jl
 
-[![Build Status](https://travis-ci.org/JuliaControl/ControlSystems.jl.svg?branch=master)](https://travis-ci.org/JuliaControl/ControlSystems.jl)
+[![Build Status](https://travis-ci.com/JuliaControl/ControlSystems.jl.svg?branch=master)](https://travis-ci.com/JuliaControl/ControlSystems.jl)
 [![PkgEval](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/C/ControlSystems.svg)](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/report.html)
 [![Gitter](https://badges.gitter.im/JuliaControl/ControlSystems.jl.svg)](https://gitter.im/JuliaControl/ControlSystems.jl?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge)
 

docs/src/man/creating_systems.md

@@ -10,7 +10,7 @@ The syntax for creating a transfer function is
 ```julia
 tf(num, den, Ts=0)
 ```
-where `num` and `den` are the polinomial coefficients of the numerator and denominator of the polynomial and `Ts` is the sample time.
+where `num` and `den` are the polynomial coefficients of the numerator and denominator of the polynomial and `Ts` is the sample time.
 ### Example:
 ```julia
 tf([1.0],[1,2,1])
@@ -29,7 +29,7 @@ Continuous-time transfer function model
 The transfer functions created using this method will be of type `TransferFunction{SisoRational}`.
 
 ## zpk - Pole-Zero-Gain Representation
-Sometimes it's better to represent the transferfunction by its poles, zeros and gain, this can be done using
+Sometimes it's better to represent the transfer function by its poles, zeros and gain, this can be done using
 ```julia
 zpk(zeros, poles, gain, Ts=0)
 ```
@@ -72,7 +72,7 @@ A state-space system is created using
 ```julia
 ss(A,B,C,D,Ts=0)
 ```
-and they behave similarily to transfer functions. State-space systems with heterogeneous matrix types are also available, which can be used to create systems with static or sized matrices, e.g.,
+and they behave similarly to transfer functions. State-space systems with heterogeneous matrix types are also available, which can be used to create systems with static or sized matrices, e.g.,
 ```jldoctest  HSS
 using StaticArrays
 import ControlSystems.HeteroStateSpace
@@ -83,7 +83,7 @@ import ControlSystems.HeteroStateSpace
 function HeteroStateSpace(A,B,C,D,Ts=0,f::F=to_static) where F
     HeteroStateSpace(f(A),f(B),f(C),f(D),Ts)
 end
-@inline HeteroStateSpace(s,f) = HeteroStateSpace(s.A,s.B,s.C,s.D,s.Ts,f)
+@inline HeteroStateSpace(s,f) = HeteroStateSpace(s.A,s.B,s.C,s.D,s.timeevol,f)
 ControlSystems._string_mat_with_headers(a::SizedArray) = ControlSystems._string_mat_with_headers(Matrix(a)); # Overload for printing purposes
 ```
 Notice the different matrix types used

src/ControlSystems.jl

@@ -6,6 +6,8 @@ export  LTISystem,
         HeteroStateSpace,
         TransferFunction,
         DelayLtiSystem,
+        Continuous,
+        Discrete,
         ss,
         tf,
         zpk,
@@ -84,7 +86,9 @@ export  LTISystem,
         numvec,
         denvec,
         numpoly,
-        denpoly
+        denpoly,
+        iscontinuous,
+        isdiscrete
 
 
 # QUESTION: are these used? LaTeXStrings, Requires, IterTools
@@ -103,6 +107,11 @@ import DSP: conv
 
 abstract type AbstractSystem end
 
+include("types/TimeEvolution.jl")
+## Added interface:
+#   timeevol(Lti) -> TimeEvolution (not exported)
+
+
 include("types/Lti.jl")
 
 include("types/SisoTf.jl")
@@ -116,6 +125,7 @@ include("types/SisoTfTypes/conversion.jl")
 
 include("types/StateSpace.jl")
 
+# TODO Sample time
 include("types/PartionedStateSpace.jl")
 include("types/DelayLtiSystem.jl")
 
@@ -154,13 +164,13 @@ include("plotting.jl")
 @deprecate num numvec
 @deprecate den denvec
 @deprecate norminf hinfnorm
+@deprecate diagonalize(s::AbstractStateSpace, digits) diagonalize(s::AbstractStateSpace)
 
 function covar(D::Union{AbstractMatrix,UniformScaling}, R)
     @warn "This call is deprecated due to ambiguity, use covar(ss(D), R) or covar(ss(D, Ts), R) instead"
     D*R*D'
 end
 
-
 # The path has to be evaluated upon initial import
 const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path())
 

src/analysis.jl

@@ -8,7 +8,7 @@ pole(sys::SisoTf) = error("pole is not implemented for type $(typeof(sys))")
 # converting to zpk before works better in the cases I have tested.
 pole(sys::TransferFunction) = pole(zpk(sys))
 
-function pole(sys::TransferFunction{SisoZpk{T,TR}}) where {T, TR}
+function pole(sys::TransferFunction{<:TimeEvolution,SisoZpk{T,TR}}) where {T, TR}
     # With right TR, this code works for any SisoTf
 
     # Calculate least common denominator of the minors,
@@ -118,7 +118,7 @@ Compute the zeros, poles, and gains of system `sys`.
 
 `k` : Matrix{Float64}, (ny x nu)"""
 function zpkdata(sys::LTISystem)
-    G = convert(TransferFunction{SisoZpk}, sys)
+    G = convert(TransferFunction{typeof(timeevol(sys)),SisoZpk}, sys)
 
     zs = map(x -> x.z, G.matrix)
     ps = map(x -> x.p, G.matrix)
@@ -133,9 +133,8 @@ Compute the natural frequencies, `Wn`, and damping ratios, `zeta`, of the
 poles, `ps`, of `sys`"""
 function damp(sys::LTISystem)
     ps = pole(sys)
-    if !iscontinuous(sys)
-        Ts = sys.Ts == -1 ? 1 : sys.Ts
-        ps = log(ps)/Ts
+    if isdiscrete(sys)
+        ps = log(ps)/sys.Ts
     end
     Wn = abs.(ps)
     order = sortperm(Wn; by=z->(abs(z), real(z), imag(z)))
@@ -446,7 +445,7 @@ function delaymargin(G::LTISystem)
     ϕₘ   *= π/180
     ωϕₘ   = m[3][i]
     dₘ    = ϕₘ/ωϕₘ
-    if G.Ts > 0
+    if isdiscrete(G)
         dₘ /= G.Ts # Give delay margin in number of sample times, as matlab does
     end
     dₘ

src/connections.jl

@@ -22,24 +22,18 @@ Append systems in block diagonal form
 """
 function append(systems::(ST where ST<:AbstractStateSpace)...)
     ST = promote_type(typeof.(systems)...)
-    Ts = systems[1].Ts
-    if !all(s.Ts == Ts for s in systems)
-        error("Sampling time mismatch")
-    end
+    timeevol = common_timeevol(systems...)
     A = blockdiag([s.A for s in systems]...)
     B = blockdiag([s.B for s in systems]...)
     C = blockdiag([s.C for s in systems]...)
     D = blockdiag([s.D for s in systems]...)
-    return ST(A, B, C, D, Ts)
+    return ST(A, B, C, D, timeevol)
 end
 
 function append(systems::TransferFunction...)
-    Ts = systems[1].Ts
-    if !all(s.Ts == Ts for s in systems)
-        error("Sampling time mismatch")
-    end
+    timeevol = common_timeevol(systems...)
     mat = blockdiag([s.matrix for s in systems]...)
-    return TransferFunction(mat, Ts)
+    return TransferFunction(mat, timeevol)
 end
 
 append(systems::LTISystem...) = append(promote(systems...)...)
@@ -64,16 +58,12 @@ function Base.vcat(systems::ST...) where ST <: AbstractStateSpace
     if !all(s.nu == nu for s in systems)
         error("All systems must have same input dimension")
     end
-    Ts = systems[1].Ts
-    if !all(s.Ts == Ts for s in systems)
-        error("Sampling time mismatch")
-    end
     A = blockdiag([s.A for s in systems]...)
     B = vcat([s.B for s in systems]...)
     C = blockdiag([s.C for s in systems]...)
     D = vcat([s.D for s in systems]...)
-
-    return ST(A, B, C, D, Ts)
+    timeevol = common_timeevol(systems...)
+    return ST(A, B, C, D, timeevol)
 end
 
 function Base.vcat(systems::TransferFunction...)
@@ -82,12 +72,10 @@ function Base.vcat(systems::TransferFunction...)
     if !all(s.nu == nu for s in systems)
         error("All systems must have same input dimension")
     end
-    Ts = systems[1].Ts
-    if !all(s.Ts == Ts for s in systems)
-        error("Sampling time mismatch")
-    end
+    timeevol = common_timeevol(systems...)
     mat = vcat([s.matrix for s in systems]...)
-    return TransferFunction(mat, Ts)
+
+    return TransferFunction(mat, timeevol)
 end
 
 Base.vcat(systems::LTISystem...) = vcat(promote(systems...)...)
@@ -98,16 +86,13 @@ function Base.hcat(systems::ST...) where ST <: AbstractStateSpace
     if !all(s.ny == ny for s in systems)
         error("All systems must have same output dimension")
     end
-    Ts = systems[1].Ts
-    if !all(s.Ts == Ts for s in systems)
-        error("Sampling time mismatch")
-    end
+    timeevol = common_timeevol(systems...)
     A = blockdiag([s.A for s in systems]...)
     B = blockdiag([s.B for s in systems]...)
     C = hcat([s.C for s in systems]...)
     D = hcat([s.D for s in systems]...)
 
-    return ST(A, B, C, D, Ts)
+    return ST(A, B, C, D, timeevol)
 end
 
 function Base.hcat(systems::TransferFunction...)
@@ -116,12 +101,10 @@ function Base.hcat(systems::TransferFunction...)
     if !all(s.ny == ny for s in systems)
         error("All systems must have same output dimension")
     end
-    Ts = systems[1].Ts
-    if !all(s.Ts == Ts for s in systems)
-        error("Sampling time mismatch")
-    end
+    timeevol = common_timeevol(systems...)
     mat = hcat([s.matrix for s in systems]...)
-    return TransferFunction(mat, Ts)
+
+    return TransferFunction(mat, timeevol)
 end
 
 Base.hcat(systems::LTISystem...) = hcat(promote(systems...)...)
@@ -139,52 +122,15 @@ function Base.typed_hcat(::Type{T}, X...) where {T<:LTISystem}
     hcat(convert.(T, X)...)
 end
 # Ambiguity
-Base.typed_hcat(::Type{T}, X::Number...) where {T<:LTISystem, N} = hcat(convert.(T, X)...)
+Base.typed_hcat(::Type{S}, X::Number...) where {S<:LTISystem} = hcat(convert.(S, X)...)
+Base.typed_hcat(::Type{S}, X::Union{AbstractArray{<:Number,1}, AbstractArray{<:Number,2}}...) where {S<:LTISystem} = hcat(convert.(S, X)...)
 
 # Catch special cases where inv(sys) might not be possible after promotion, like improper tf
 function /(sys1::Union{StateSpace,AbstractStateSpace}, sys2::LTISystem)
     sys1new, sys2new = promote(sys1, 1/sys2)
     return sys1new*sys2new
 end
 
-# function hvcat(rows::Tuple{Vararg{Int}}, systems::Union{Number,AbstractVecOrMat{<:Number},LTISystem}...)
-#     T = Base.promote_typeof(systems...)
-#     nbr = length(rows)  # number of block rows
-#     rs = Array{T,1}(nbr)
-#     a = 1
-#     for i = 1:nbr
-#         rs[i] = hcat(convert.(T,systems[a:a-1+rows[i]])...)
-#         a += rows[i]
-#     end
-#     vcat(rs...)
-# end
-
-# function _get_common_sampling_time(sys_vec::Union{AbstractVector{LTISystem},AbstractVecOrMat{<:Number},Number})
-#     Ts = -1.0 # Initalize corresponding to undefined sampling time
-#
-#     for sys in sys_vec
-#         if !all(s.Ts == Ts for s in systems])
-#             error("Sampling time mismatch")
-#         end
-#     end
-#
-# end
-
-
-# function Base.hcat{T<:Number}(systems::Union{T,AbstractVecOrMat{T},TransferFunction}...)
-#     S = promote_type(map(e->typeof(e),systems)...) # TODO: Should be simplified
-#
-#     idx_first_tf = findfirst(e -> isa(e, TransferFunction), systems)
-#     Ts = sys_tuple[idx_first_tf].Ts
-#
-#     if S <: TransferFunction
-#         hcat(map(e->convert(TransferFunction,e),systems)...)
-#     else
-#         cat(2,systems...)
-#     end
-# end
-
-
 blockdiag(mats::AbstractMatrix...) = blockdiag(promote(mats...)...)
 
 function blockdiag(mats::AbstractMatrix{T}...) where T
@@ -214,16 +160,16 @@ feedback(L::TransferFunction) = L/(1+L)
 feedback(P1::TransferFunction, P2::TransferFunction) = P1/(1+P1*P2)
 
 #Efficient implementations
-function feedback(L::TransferFunction{T}) where T<:SisoRational
+function feedback(L::TransferFunction{<:TimeEvolution,T}) where T<:SisoRational
     if size(L) != (1,1)
         error("MIMO TransferFunction feedback isn't implemented, use L/(1+L)")
     end
     P = numpoly(L)
     Q = denpoly(L)
-    tf(P, P+Q, L.Ts)
+    tf(P, P+Q, L.timeevol)
 end
 
-function feedback(L::TransferFunction{T}) where {T<:SisoZpk}
+function feedback(L::TransferFunction{TE, T}) where {TE<:TimeEvolution, T<:SisoZpk}
     if size(L) != (1,1)
         error("MIMO TransferFunction feedback isn't implemented, use L/(1+L)")
     end
@@ -233,7 +179,7 @@ function feedback(L::TransferFunction{T}) where {T<:SisoZpk}
     kden = denpol[end] # Get coeff of s^n
     # Create siso system
     sisozpk = T(L.matrix[1].z, roots(denpol), k/kden)
-    return TransferFunction{T}(fill(sisozpk,1,1), L.Ts)
+    return TransferFunction{TE,T}(fill(sisozpk,1,1), L.timeevol)
 end
 
 """
@@ -255,17 +201,13 @@ function feedback(sys::Union{StateSpace, DelayLtiSystem})
 end
 
 function feedback(sys1::StateSpace,sys2::StateSpace)
-    if sys1.Ts != sys2.Ts # FIXME: replace with common_sample_time
-        error("Sampling time mismatch")
-    end
-
+    timeevol = common_timeevol(sys1,sys2)
     !(iszero(sys1.D) || iszero(sys2.D)) && error("There cannot be a direct term (D) in both sys1 and sys2")
     A = [sys1.A+sys1.B*(-sys2.D)*sys1.C sys1.B*(-sys2.C);
          sys2.B*sys1.C  sys2.A+sys2.B*sys1.D*(-sys2.C)]
     B = [sys1.B; sys2.B*sys1.D]
     C = [sys1.C  sys1.D*(-sys2.C)]
-
-    ss(A, B, C, sys1.D, sys1.Ts)
+    ss(A, B, C, sys1.D, timeevol)
 end
 
 
@@ -288,9 +230,7 @@ See Zhou etc. for similar (somewhat less symmetric) formulas.
     U1=:, Y1=:, U2=:, Y2=:, W1=:, Z1=:, W2=Int[], Z2=Int[],
     Wperm=:, Zperm=:, pos_feedback::Bool=false)
 
-    if sys1.Ts != sys2.Ts # FIXME: replace with common_sample_time
-        error("Sampling time mismatch")
-    end
+    timeevol = common_timeevol(sys1,sys2)
 
     if !(isa(Y1, Colon) || allunique(Y1)); @warn "Connecting single output to multiple inputs Y1=$Y1"; end
     if !(isa(Y2, Colon) || allunique(Y2)); @warn "Connecting single output to multiple inputs Y2=$Y2"; end
@@ -359,7 +299,7 @@ See Zhou etc. for similar (somewhat less symmetric) formulas.
                      s2_D12*R2*s1_D21           s2_D11 + α*s2_D12*R2*s1_D22*s2_D21]
     end
 
-    return StateSpace(A, B[:, Wperm], C[Zperm,:], D[Zperm, Wperm], sys1.Ts)
+    return StateSpace(A, B[:, Wperm], C[Zperm,:], D[Zperm, Wperm], timeevol)
 end
 
 
@@ -369,7 +309,7 @@ end
 """
 function feedback2dof(P::TransferFunction,R,S,T)
     !issiso(P) && error("Feedback not implemented for MIMO systems")
-    tf(conv(poly2vec(numpoly(P)[1]),T),zpconv(poly2vec(denpoly(P)[1]),R,poly2vec(numpoly(P)[1]),S))
+    tf(conv(poly2vec(numpoly(P)[1]),T),zpconv(poly2vec(denpoly(P)[1]),R,poly2vec(numpoly(P)[1]),S), P.timeevol)
 end
 
 feedback2dof(B,A,R,S,T) = tf(conv(B,T),zpconv(A,R,B,S))