diff --git a/lib/ControlSystemsBase/src/discrete.jl b/lib/ControlSystemsBase/src/discrete.jl index acd938853..680842784 100644 --- a/lib/ControlSystemsBase/src/discrete.jl +++ b/lib/ControlSystemsBase/src/discrete.jl @@ -67,6 +67,13 @@ function c2d_x0map(sys::AbstractStateSpace{<:Continuous}, Ts::Real, method::Symb elseif method === :fwdeuler Ad, Bd, Cd, Dd = (I+Ts*A), Ts*B, C, D x0map = I(nx) + elseif method === :bwdeuler + Ad = inv(I - Ts*A) + TsB = Ts*B + Bd = (Ad * TsB) + Cd = C*Ad + Dd = Cd * TsB + D + x0map = I(nx) elseif method === :tustin a = w_prewarp == 0 ? Ts/2 : tan(w_prewarp*Ts/2)/w_prewarp a > 0 || throw(DomainError("A positive w_prewarp must be provided for method Tustin")) @@ -77,7 +84,7 @@ function c2d_x0map(sys::AbstractStateSpace{<:Continuous}, Ts::Real, method::Symb Dd = a*Cd*B + D x0map = Matrix{T}(I, nx, nx) elseif method === :matched - error("NotImplemented: Only `:zoh`, `:foh`, :tustin and `:fwdeuler` implemented so far") + error("NotImplemented: Only `:zoh`, `:foh`, :tustin, `:fwdeuler` and `bwdeuler` implemented so far") else error("Unsupported method: ", method) end @@ -111,6 +118,11 @@ function d2c(sys::AbstractStateSpace{<:Discrete}, method::Symbol=:zoh; w_prewarp Ac = (A-I)./sys.Ts Bc = B./sys.Ts Cc, Dc = C, D + elseif method === :bwdeuler + Ac = (I - inv(A))/sys.Ts + Bc = (A\B) ./ sys.Ts + Cc = C/A + Dc = D - sys.Ts*C*Bc elseif method === :tustin a = w_prewarp == 0 ? sys.Ts/2 : tan(w_prewarp*sys.Ts/2)/w_prewarp a > 0 || throw(DomainError("A positive w_prewarp must be provided for method Tustin")) diff --git a/lib/ControlSystemsBase/test/test_discrete.jl b/lib/ControlSystemsBase/test/test_discrete.jl index c02a2ade0..4e304b839 100644 --- a/lib/ControlSystemsBase/test/test_discrete.jl +++ b/lib/ControlSystemsBase/test/test_discrete.jl @@ -77,9 +77,9 @@ sysd = c2d(sys, 1) @test d2c(sysd) ≈ sys -# forward euler / tustin +# forward euler / tustin / backward euler @test c2d(C_111, 1, :fwdeuler).A == I + C_111.A -for method in (:fwdeuler, :tustin) +for method in (:fwdeuler, :tustin, :bwdeuler) @test d2c(c2d(C_111, 0.01, method), method) ≈ C_111 atol = sqrt(eps()) @test d2c(c2d(C_212, 0.01, method), method) ≈ C_212 atol = sqrt(eps()) @test d2c(c2d(C_221, 0.01, method), method) ≈ C_221 atol = sqrt(eps())