add d2c_exact

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,28 @@ end
 
 d2c(sys::TransferFunction{<:Discrete}, args...) = tf(d2c(ss(sys), args...))
 
+
+"""
+    d2c_exact(sys::AbstractStateSpace{<:Discrete})
+
+Translate a discrete-time system to a continuous-time system by the substitution ``z = e^{sT_s}``.
+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 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 negative delays,
+i.e., this system cannot be simulated in the time domain.
+"""
+function d2c_exact(sys::AbstractStateSpace{<:Discrete})
+    T = sys.Ts
+    A,B,C,D = ssdata(sys)
+    z = delay(-T)
+    LR = append([z for _ in 1:sys.nx]...) - ss(A + I)
+    C*feedback(I(sys.nx), LR)*B + D
+end
+
 # c2d and d2c for covariance and cost matrices =================================
 """
     Qd     = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Ts;             opt=:o)

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) = basetype(sys)(zero(sys.D), sys.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_discrete.jl

@@ -146,6 +146,16 @@ 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 = d2c_exact(sys)
+
+    # bodeplot(sys, w, lab="Discrete SS")
+    # bodeplot!(sysc, w, lab="Continuous SS with delays", l=:dash, size=(800, 800))
+    w = exp10.(LinRange(-2, 2, 200))
+    @test freqresp(sys, w) ≈ freqresp(sysc, w) atol = 1e-8
 end
 
 
@@ -233,3 +243,4 @@ Qd_van_load = [
 @test norm(Qd - Qd_van_load) < 1e-6
 
 
+end