Merge pull request #220 from JuliaControl/moredelays

Complete version of delays

README.md

@@ -17,6 +17,19 @@ Pkg.add("ControlSystems")
 ```
 
 ## News
+### 2019-05-28
+#### Delay systems
+- We now support systems with time delays. Example:
+```julia
+sys = tf(1, [1,1])*delay(1)
+stepplot(sys, 5) # Compilation time might be long for first simulation
+nyquistplot(sys)
+```
+#### New examples
+- [Delayed systems (frequency domain)](https://github.com/JuliaControl/ControlSystems.jl/blob/master/example/delayed_lti_system.jl)
+- [Delayed systems (time domain)](https://github.com/JuliaControl/ControlSystems.jl/blob/master/example/delayed_lti_timeresp.jl)
+- [Systems with uncertainty](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/examples/controlsystems.jl)
+- [Robust PID optimization](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/examples/robust_controller_opt.jl)
 ### 2019-05-22
 New state-space type `HeteroStateSpace` that accepts matrices of heterogeneous types: [example using `StaticArrays`](https://juliacontrol.github.io/ControlSystems.jl/latest/man/creating_systems/#Creating-State-Space-Systems-1).
 ### 2019-01-31

src/connections.jl

@@ -135,7 +135,9 @@ function Base._cat_t(::Val{2}, T::Type{<:LTISystem}, X...)
 end
 
 # Used in typed_hvcat
-Base.typed_hcat(::Type{T}, X...) where {T<:LTISystem} = hcat(convert.(T, X)...)
+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)...)
 

src/delay_systems.jl

@@ -23,7 +23,7 @@ end
 
 
 """
-    `y, t, x = lsim(sys::DelayLtiSystem, t::AbstractArray{<:Real}; u=(out, t) -> (out .= 0), x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...)`
+    `y, t, x = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...)`
 
     Simulate system `sys`, over time `t`, using input signal `u`, with initial state `x0`, using method `alg` .
 
@@ -38,7 +38,7 @@ end
 
     Returns: times `t`, and `y` and `x` at those times.
 """
-function lsim(sys::DelayLtiSystem{T}, u, t::AbstractArray{<:T}; x0=fill(zero(T), nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...) where T
+function lsim(sys::DelayLtiSystem{T,S}, u, t::AbstractArray{<:Real}; x0=fill(zero(T), nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...) where {T,S}
     # Make u! in-place function of u
     u! = if isa(u, Number) || isa(u,AbstractVector) # Allow for u to be a constant number or vector
         (uout, t) -> uout .= u
@@ -70,7 +70,7 @@ function dde_param(dx, x, h!, p, t)
 end
 
 # TODO Discontinuities in u are not handled well yet.
-function _lsim(sys::DelayLtiSystem{T}, Base.@nospecialize(u!), t::AbstractArray{<:T}, x0::Vector{T}, alg; kwargs...) where T
+function _lsim(sys::DelayLtiSystem{T,S}, Base.@nospecialize(u!), t::AbstractArray{<:Real}, x0::Vector{T}, alg; kwargs...) where {T,S}
     P = sys.P
 
     if ~iszero(P.D22)
@@ -94,22 +94,22 @@ function _lsim(sys::DelayLtiSystem{T}, Base.@nospecialize(u!), t::AbstractArray{
     h!(out, p, t) = (out .= 0)      # History function
 
     p = (A, B1, B2, C2, D21, Tau, u!, uout, hout, tmp)
-    prob = DDEProblem{true}(dde_param, x0, h!, (t[1], t[end]), p, constant_lags=sys.Tau)
+    prob = DDEProblem{true}(dde_param, x0, h!, (T(t[1]), T(t[end])), p, constant_lags=sys.Tau)
 
     sol = DelayDiffEq.solve(prob, alg; saveat=t, kwargs...)
 
-    x = sol.u::Array{Array{T,1},1} # the states are labeled u in DelayDiffEq
+    x = sol.u::Vector{Vector{T}} # the states are labeled u in DelayDiffEq
 
-    y = Array{T,2}(undef, noutputs(sys), length(t))
-    d = Array{T,2}(undef, size(C2,1), length(t))
+    y = Matrix{T}(undef, noutputs(sys), length(t))
+    d = Matrix{T}(undef, size(C2,1), length(t))
     # Build y signal (without d term)
     for k = 1:length(t)
         u!(uout, t[k])
         y[:,k] = C1*x[k] + D11*uout
         #d[:,k] = C2*x[k] + D21*uout
     end
 
-    dtmp = Array{T,1}(undef, size(C2,1))
+    dtmp = Vector{T}(undef, size(C2,1))
     # Function to evaluate d(t)_i at an arbitrary time
     # X is continuous, so interpoate, u is not
     function dfunc!(tmp::Array{T1}, t, i) where T1
@@ -135,7 +135,7 @@ end
 
 
 # We have to default to something, look at the sys.P.P and delays
-function _bounds_and_features(sys::DelayLtiSystem, plot::Symbol)
+function _bounds_and_features(sys::DelayLtiSystem, plot::Val)
     ws, pz =  _bounds_and_features(sys.P.P, plot)
     logtau = log10.(abs.(sys.Tau))
     logtau = filter(x->x>4, logtau) # Ignore low frequency

src/synthesis.jl

@@ -192,9 +192,9 @@ Forms the negative feedback interconnection
 ```
 If no second system is given, negative identity feedback is assumed
 """
-function feedback(sys::StateSpace)
-    sys.ny != sys.nu && error("Use feedback(sys1::StateSpace,sys2::StateSpace) if sys.ny != sys.nu")
-    feedback(sys,ss(Matrix{numeric_type(sys)}(I,sys.ny,sys.ny)))
+function feedback(sys::Union{StateSpace, DelayLtiSystem})
+    ninputs(sys) != noutputs(sys) && error("Use feedback(sys1, sys2) if number of inputs != outputs")
+    feedback(sys,ss(Matrix{numeric_type(sys)}(I,size(sys)...)))
 end
 
 function feedback(sys1::StateSpace,sys2::StateSpace)

src/types/DelayLtiSystem.jl

@@ -1,3 +1,8 @@
+"""
+    struct DelayLtiSystem{T, S <: Real} <: LTISystem
+
+Represents an LTISystem with internal time-delay. See `?delay` for a convenience constructor.
+"""
 struct DelayLtiSystem{T,S<:Real} <: LTISystem
     P::PartionedStateSpace{StateSpace{T,Matrix{T}}}
     Tau::Vector{S} # The length of the vector tau implicitly defines the partitionging of P
@@ -13,6 +18,11 @@ end
 
 # QUESTION: would psys be a good standard variable name for a PartionedStateSpace
 #           and perhaps dsys for a delayed system, (ambigous with discrete system though)
+"""
+    DelayLtiSystem{T, S}(sys::StateSpace, Tau::AbstractVector{S}=Float64[]) where {T <: Number, S <: Real}
+
+Create a delayed system by speciying both the system and time-delay vector. NOTE: if you want to create a system with simple input or output delays, use the Function `delay(τ)`.
+"""
 function DelayLtiSystem{T,S}(sys::StateSpace, Tau::AbstractVector{S} = Float64[]) where {T<:Number,S<:Real}
     nu = ninputs(sys) - length(Tau)
     ny = noutputs(sys) - length(Tau)
@@ -49,6 +59,10 @@ end
 function Base.convert(::Type{DelayLtiSystem{T1,S}}, d::T2) where {T1,T2 <: Number,S}
     DelayLtiSystem{T1,S}(StateSpace(T1(d)))
 end
+function Base.convert(::Type{DelayLtiSystem{T1,S}}, d::T2) where {T1,T2 <: AbstractArray,S}
+    DelayLtiSystem{T1,S}(StateSpace(T1.(d)))
+end
+
 Base.convert(::Type{<:DelayLtiSystem}, sys::TransferFunction)  = DelayLtiSystem(sys)
 # Catch convertsion between T
 Base.convert(::Type{V}, sys::DelayLtiSystem)  where {T, V<:DelayLtiSystem{T}} =
@@ -156,10 +170,19 @@ function feedback(sys1::DelayLtiSystem, sys2::DelayLtiSystem)
     DelayLtiSystem(psys_new.P, Tau_new)
 end
 
-function delay(tau::S, T::Type{<:Number}=Float64) where S
+"""
+    delay(T::Type{<:Number}, tau)
+
+Create a pure time delay. If `T` is not specified, the default is to choose `promote_type(T, typeof(tau))`
+"""
+function delay(T::Type{<:Number}, tau)
     return DelayLtiSystem(ControlSystems.ss([zero(T) one(T); one(T) zero(T)]), [T(tau)])
 end
 
+function delay(tau::S) where S
+    delay(promote_type(Float64,S), tau)
+end
+
 # function exp(G::TransferFunction)
 #     if (size(G.matrix) != [1, 1]) || ~isone(G.matrix[1].den) || length(G.matrix[1].num) >= 2
 #         error("exp only accepts TransferFunction arguemns of the form (a*s + b)")

src/types/PartionedStateSpace.jl

@@ -211,3 +211,7 @@ function hcat_1(systems::PartionedStateSpace...)
     sysnew = StateSpace(A, [B1 B2], [C1; C2], [D11 D12; D21 D22], Ts)
     return PartionedStateSpace(sysnew, sum(s -> s.nu1, systems), ny1)
 end
+
+function Base.convert(::Type{<:PartionedStateSpace}, sys::T) where T<: StateSpace
+    PartionedStateSpace(sys,sys.nu,sys.ny)
+end

src/types/promotion.jl

@@ -21,7 +21,7 @@
 
 # Abstract type pyramid =============================================================
 
-# NOTE: mycket som inte verkar användas...
+# NOTE: a lot that doesnt seem to be used...
 #Base.promote_rule(::Type{StateSpace{T1}}, ::Type{StateSpace{T2}}) where {T1,T2} = StateSpace{promote_type(T1, T2)}
 
 # NOTE: Is the below thing correct always?
@@ -84,8 +84,8 @@ function Base.promote_rule(::Type{StateSpace{T1, MT1}}, ::Type{MT2}) where {T1,
     StateSpace{eltype(MT), MT}
 end
 
-Base.promote_rule(::Type{DelayLtiSystem{T1}}, ::Type{MT1}) where {T1, MT1<:AbstractMatrix} =
-    DelayLtiSystem{promote_type(T1, eltype(MT1))}
+Base.promote_rule(::Type{DelayLtiSystem{T1,S}}, ::Type{MT1}) where {T1, S, MT1<:AbstractMatrix} =
+    DelayLtiSystem{promote_type(T1, eltype(MT1)),S}
 
 #Base.promote_rule{S<:TransferFunction{<:SisoTf}}(::Type{S}, ::Type{<:Real}) = S
 

src/utilities.jl

@@ -6,6 +6,7 @@ numeric_type(sys::SisoTf) = numeric_type(typeof(sys))
 
 numeric_type(::Type{TransferFunction{S}}) where S = numeric_type(S)
 numeric_type(::Type{<:StateSpace{T}}) where T = T
+numeric_type(::Type{<:DelayLtiSystem{T}}) where {T} = T
 numeric_type(sys::LTISystem) = numeric_type(typeof(sys))
 
 

test/test_delayed_systems.jl

@@ -116,6 +116,9 @@ w = 10 .^ (-2:0.1:2)
 
 # Test step
 println("Simulating first delay system:")
+@time step(delay(1)*tf(1,[1.,1]))
+@time step(delay(1)*tf(1,[1,1]))
+
 @time y1, t1, x1 = step([s11;s12], 10)
 @time @test y1[:,2] ≈ step(s12, t1)[1] rtol = 1e-14