Merge pull request #902 from JuliaControl/d2c_exact

Add exact frequency-domain translation of discrete systems to continuous time

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -84,6 +84,7 @@ export  LTISystem,
         c2d,
         c2d_x0map,
         d2c,
+        d2c_exact,
         # Time Response
         step,
         impulse,

lib/ControlSystemsBase/src/discrete.jl

@@ -91,6 +91,8 @@ end
 Convert discrete-time system to a continuous time system, assuming that the discrete-time system was discretized using `method`. Available methods are `:zoh, :fwdeuler´.
 
 - `w_prewarp`: Frequency for pre-warping when using the Tustin method, has no effect for other methods.
+
+See also [`d2c_exact`](@ref).
 """
 function d2c(sys::AbstractStateSpace{<:Discrete}, method::Symbol=:zoh; w_prewarp=0)
     A, B, C, D = ssdata(sys)
@@ -125,6 +127,38 @@ end
 
 d2c(sys::TransferFunction{<:Discrete}, args...) = tf(d2c(ss(sys), args...))
 
+
+"""
+    d2c_exact(sys::AbstractStateSpace{<:Discrete}, method = :causal)
+
+Translate a discrete-time system to a continuous-time system by one of the substitutions
+- ``z^{-1} = e^{-sT_s}`` if `method = :causal` (default)
+- ``z = e^{sT_s}``  if `method = :acausal`
+The translation is exact in the frequency domain, i.e.,
+the frequency response of the resulting continuous-time system is identical to
+the frequency response of the discrete-time system.
+
+This method of translation is useful when analyzing hybrid continuous/discrete systems in the frequency domain and high accuracy is required.
+
+The resulting system will be be a static system in feedback with pure delays. When `method = :causal`, the delays will be positive, resulting in a causal system that can be simulated in the time domain. When `method = :acausal`, the delays will be negative, resulting in an acausal system that **can not** be simulated in the time domain. The acausal translation results in a smaller system with half as many delay elements in the feedback path.
+"""
+function d2c_exact(sys::AbstractStateSpace{<:Discrete}, method=:causal)
+    T = sys.Ts
+    A,B,C,D = ssdata(sys)
+    if method === :acausal
+        z = delay(-T)
+        LR = append([z for _ in 1:sys.nx]...) - ss(A + I)
+        C*feedback(I(sys.nx), LR)*B + D
+    elseif method === :causal
+        z1 = delay(T)
+        ZI1 = append([z1 for _ in 1:sys.nx]...)
+        LR = ZI1 * ss(-A)
+        C*ZI1*feedback(I(sys.nx), LR)*B + D
+    else
+        error("Unknown method: $method. Choose method = :acausal or :causal")
+    end
+end
+
 # c2d and d2c for covariance and cost matrices =================================
 """
     Qd     = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Ts;             opt=:o)

lib/ControlSystemsBase/src/types/DelayLtiSystem.jl

@@ -94,9 +94,16 @@ function +(sys::DelayLtiSystem{T1,S}, n::T2) where {T1,T2<:Number,S}
 end
 
 
+"""
+    /(G1, G2::DelayLtiSystem)
+
+Compute ``G_1 * G_2^{-1} where ``G_2`` is a DelayLtiSystem.
+Throws a SingularException if ``G_2`` is not invertible.
+"""
 function /(anything, sys::DelayLtiSystem)
-    all(iszero, sys.Tau) || error("A delayed system can not be inverted. Consider use of the function `feedback`.")
-    /(anything, sys.P.P) # If all delays are zero, invert the inner system
+    ny,nu = size(sys)
+    ny == nu || error("The denominator system must be square")
+    return anything * feedback(I(nu), sys - I(nu))
 end
 
 for other_type in [:Number, :AbstractMatrix, :LTISystem]

lib/ControlSystemsBase/src/types/LFTSystem.jl

@@ -1,6 +1,10 @@
 abstract type LFTSystem{TE, T} <: LTISystem{TE} end
 timeevol(sys::LFTSystem) = timeevol(sys.P)
 
+Base.zero(sys::LFTSystem) = ss(zero(sys.P.D), sys.P.timeevol)
+Base.zero(::Type{<:LFTSystem{Continuous, F}}) where {F} = ss([zero(F)], Continuous())
+
+
 function *(n::Number, sys::LFTSystem)
     new_C = [sys.P.C1*n; sys.P.C2]
     new_D = [sys.P.D11*n sys.P.D12*n; sys.P.D21 sys.P.D22]

lib/ControlSystemsBase/test/test_delayed_systems.jl

@@ -56,6 +56,13 @@ P2_fr = (im*ω .+ 1) ./ (im*ω .+ 2)
 @test freqresp(P2 * delay(1), ω)[:] ≈ P2_fr .* exp.(-im*ω) rtol=1e-15
 @test freqresp(delay(1) * P2, ω)[:] ≈ P2_fr .* exp.(-im*ω) rtol=1e-15
 
+## Division / feedback
+@test freqresp(1/(1+P1), ω) ≈ freqresp(feedback(I(size(P1, 1)), P1), ω) rtol=1e-15
+
+## Zero
+@test zero(typeof(P1)) == ss(0.0)
+@test zero(P1) == ss(0.0)
+
 ## append
 P12 = append(P1, P2)
 G12 = [P1 tf(0); tf(0) P2]
@@ -179,9 +186,10 @@ w = 10 .^ (-2:0.1:2)
 
 @test propertynames(delay(1.0)) == (:P, :Tau, :nu, :ny)
 
-@test_throws ErrorException 1/f2
-@test_throws ErrorException randn(2,2)/f2
-@test_throws ErrorException f2/f2
+@test_throws SingularException 1/f2
+@test_throws SingularException randn(2,2)/f2
+@test_throws SingularException f2/f2
+@test 1/(I(2)+f2) == feedback(I(2), f2)
 
 #FIXME: A lot more tests, including MIMO systems in particular
 

lib/ControlSystemsBase/test/test_discrete.jl

@@ -146,6 +146,24 @@ p = poles(Gcl)
 pd = c2d_poly2poly(A, 0.1)
 @test pd ≈ denvec(c2d(P, 0.1))[]
 
+
+@testset "d2c_exact" begin
+    @info "Testing d2c_exact"
+    sys = ssrand(2, 3, 2, Ts = 1)
+    sysc_causal = d2c_exact(sys, :causal)
+    sysc_acausal = d2c_exact(sys, :acausal)
+    @test_throws ErrorException d2c_exact(sys, :KentMorgan)
+    # bodeplot(sys, w, lab="Discrete SS")
+    # bodeplot!(sysc_causal, w, lab="Continuous SS with causal delays", l=:dash, size=(800, 800))
+    # bodeplot!(sysc_acausal, w, lab="Continuous SS with acausal (negative) delays", l=:dash, size=(800, 800))
+    w = exp10.(LinRange(-2, 2, 200))
+    @test freqresp(sys, w) ≈ freqresp(sysc_causal, w) atol = 1e-8
+    @test freqresp(sys, w) ≈ freqresp(sysc_acausal, w) atol = 1e-8
+
+    # Requires full ControlSystems but does pass
+    # rd = step(sys, 0:10)
+    # rc = step(sysc_causal, 0:10)
+    # @test rd.y ≈ rc.y atol = 1e-8
 end
 
 
@@ -233,3 +251,4 @@ Qd_van_load = [
 @test norm(Qd - Qd_van_load) < 1e-6
 
 
+end