Basic c2d for systems with delays (#343)

Co-authored-by: olof3 <olof3@control.lth.se>

src/delay_systems.jl

@@ -36,6 +36,35 @@ function evalfr(sys::DelayLtiSystem, s)
     return P11_fr + P12_fr/(delay_matrix_inv_fr - P22_fr)*P21_fr
 end
 
+
+"""
+    delayd_ss(τ, Ts)
+Discrete-time statespace realization of a delay τ sampled with period Ts,
+i.e. of z^-N where N = τ / Ts.
+
+τ must be a multiple of Ts.
+"""
+function delayd_ss(τ::Number, h::Number)
+    n = Int(round(τ / h))
+    if !(τ - n*h ≈ 0)
+        error("The delay τ must be a multiple of the sample time Ts")
+    end
+    ss(diagm(1 => ones(n-1)), [zeros(n-1,1); 1], [1 zeros(1,n-1)], 0, h)
+end
+
+"""
+    c2d(G::DelayLtiSystem, Ts, method=:zoh)
+"""
+function c2d(G::DelayLtiSystem, h::Real, method=:zoh)
+    if !(method === :zoh)
+        error("c2d for DelayLtiSystems only supports zero-order hold")
+    end
+    X = append([delayd_ss(τ, h) for τ in G.Tau]...)
+    Pd = c2d(G.P.P, h)[1]
+    return lft(Pd, X)
+end
+
+
 """
     `y, t, x = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...)`
 

test/test_delayed_systems.jl

@@ -97,6 +97,19 @@ d = exp(-2*s)
 @test_throws ErrorException exp([-2*s; -s])
 @test_throws ErrorException exp(2*s) # Non-causal
 
+
+Ω = [0.5, 1.0, 1.5, 2.0]
+# Test for internal function delayd_ss
+@test freqresp(ControlSystems.delayd_ss(1.0, 0.2), Ω)[:] ≈ exp.(-im*Ω) atol=1e-14
+@test freqresp(ControlSystems.delayd_ss(3.2, 0.4), Ω)[:] ≈ exp.(-3.2*im*Ω) atol=1e-14
+@test_throws ErrorException freqresp(ControlSystems.delayd_ss(3.2, 0.5), Ω)
+
+# Simple tests for c2d of DelayLtiSystems
+@test freqresp(c2d(feedback(ss(0,1,1,0), delay(1.5)), 0.5), Ω) ≈ [0.5/((z - 1) + 0.5*z^-3) for z in exp.(im*Ω*0.5)]
+@test freqresp(c2d(feedback(delay(1.5), delay(1.0)), 0.5), Ω) ≈ [z^-3/(1 + z^-5) for z in exp.(im*Ω*0.5)]
+@test freqresp(c2d(feedback(0.5, delay(1.5)), 0.5), Ω) ≈ [0.5/(1 + 0.5*z^-3) for z in exp.(im*Ω*0.5)]
+
+
 # Random conversions
 sys1 = DelayLtiSystem(1.0/s)
 @test sys1.P.A == sys1.P.D == fill(0,1,1)