WIP: C2d single returnval (#436)

* Basic c2d for systems with delays

* single returnval for c2d of StateSpace

* Fixes to tests.

Co-authored-by: olof3 <[email protected]>

Author: olof3

src/ControlSystems.jl

@@ -62,6 +62,7 @@ export  LTISystem,
         lft,
         # Discrete
         c2d,
+        c2d_x0map,
         d2c,
         # Time Response
         step,

src/delay_systems.jl

@@ -60,7 +60,7 @@ function c2d(G::DelayLtiSystem, h::Real, 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]
+    Pd = c2d(G.P.P, h)
     return lft(Pd, X)
 end
 

src/discrete.jl

@@ -2,16 +2,27 @@ export rstd, rstc, dab, c2d_roots2poly, c2d_poly2poly, zpconv#, lsima, indirect_
 
 
 """
-    sysd, x0map = c2d(sys::StateSpace, Ts, method=:zoh)
-    sysd = c2d(sys::TransferFunction, Ts, method=:zoh)
+    sysd= c2d(sys::StateSpace, Ts, method=:zoh)
+    Gd = c2d(G::TransferFunction, Ts, method=:zoh)
 
-Convert the continuous system `sys` into a discrete system with sample time
-`Ts`, using the provided method. Currently only `:zoh`, `:foh` and `:fwdeuler` are provided. Note that the forward-Euler method generally requires the sample time to be very small in relation to the time-constants of the system.
+Convert the continuous-time system `sys` into a discrete-time system with sample time
+`Ts`, using the specified `method` (:zoh`, `:foh` or `:fwdeuler`).
+Note that the forward-Euler method generally requires the sample time to be very small
+relative to the time constants of the system.
 
-Returns the discrete system `sysd`, and for StateSpace systems a matrix `x0map` that transforms the
-initial conditions to the discrete domain by
-`x0_discrete = x0map*[x0; u0]`"""
-function c2d(sys::StateSpace{Continuous}, Ts::Real, method::Symbol=:zoh)
+See also `c2d_x0map`
+"""
+c2d(sys::StateSpace, Ts::Real, method::Symbol=:zoh) = c2d_x0map(sys, Ts, method)[1]
+
+
+"""
+    sysd, x0map = c2d_x0map(sys::StateSpace, Ts, method=:zoh)
+
+Returns the discretization `sysd` of the system `sys` and a matrix `x0map` that
+transforms the initial conditions to the discrete domain by `x0_discrete = x0map*[x0; u0]`
+
+See `c2d` for further details."""
+function c2d_x0map(sys::StateSpace{Continuous}, Ts::Real, method::Symbol=:zoh)
     A, B, C, D = ssdata(sys)
     T = promote_type(eltype.((A,B,C,D))...)
     ny, nu = size(sys)
@@ -220,13 +231,13 @@ function c2d(G::TransferFunction{<:Continuous}, h, args...; kwargs...)
     ny, nu = size(G)
     @assert (ny + nu == 2) "c2d(G::TransferFunction, h) not implemented for MIMO systems"
     sys = ss(G)
-    sysd = c2d(sys, h, args...; kwargs...)[1]
+    sysd = c2d(sys, h, args...; kwargs...)
     return convert(TransferFunction, sysd)
 end
 
 """
     zpc(a,r,b,s)
-    
+
 form conv(a,r) + conv(b,s) where the lengths of the polynomials are equalized by zero-padding such that the addition can be carried out
 """
 function zpconv(a,r,b,s)

src/timeresp.jl

@@ -132,9 +132,9 @@ function lsim(sys::AbstractStateSpace, u::AbstractVecOrMat, t::AbstractVector;
         end
 
         if method === :zoh
-            dsys = c2d(sys, dt, :zoh)[1]
+            dsys = c2d(sys, dt, :zoh)
         elseif method === :foh
-            dsys, x0map = c2d(sys, dt, :foh)
+            dsys, x0map = c2d_x0map(sys, dt, :foh)
             x0 = x0map*[x0; transpose(u)[:,1]]
         else
             error("Unsupported discretization method")
@@ -180,7 +180,7 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
 
     if !iscontinuous(sys) || method === :zoh
         if iscontinuous(sys)
-            dsys = c2d(sys, dt, :zoh)[1]
+            dsys = c2d(sys, dt, :zoh)
         else
             if sys.Ts != dt
                 error("Time vector must match sample time for discrete system")
@@ -224,7 +224,7 @@ If `u` is a function, then `u(x,i)` is called to calculate the control signal ev
                       LinearAlgebra.promote_op(LinearAlgebra.matprod, eltype(B), eltype(u)))
 
     n = size(u, 1)
-  
+
     # Transposing u allows column-wise operations, which apparently is faster.
     ut = transpose(u)
 

test/framework.jl

@@ -70,20 +70,16 @@ macro test_array_vecs_eps(a, b, tol)
     end
 end
 
-macro test_c2d(ex, sys_sol, mat_sol, op)
-    if op == true
-        quote
-            sys, mat = $(esc(ex))
-            @test sys ≈ $(esc(sys_sol)) && mat ≈ $(esc(mat_sol))
-        end
-    else
-        quote
-            sys, mat = $(esc(ex))
-            @test !(sys ≈ $(esc(sys_sol))) || !(mat ≈ $(esc(mat_sol)))
-        end
+# Currently only works for two return values
+macro test_tupleapprox(ex, y1true, y2true)
+    quote
+        y1, y2 = $(esc(ex))
+        @test y1 ≈ $(esc(y1true)) && y2 ≈ $(esc(y2true))
     end
 end
 
 
 approxin(el,col;kwargs...) = any(colel -> isapprox(el, colel; kwargs...), col)
 approxsetequal(s1,s2;kwargs...) = all(approxin(p,s1;kwargs...) for p in s2) && all(approxin(p,s2;kwargs...) for p in s1)
+
+dropwhitespace(str) = filter(!isspace, str)

test/test_discrete.jl

@@ -5,60 +5,49 @@ C_212 = ss([-5 -3; 2 -9], [1; 2], [1 0; 0 1], [0; 0])
 C_221 = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0], [0 0])
 C_222_d = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0; 0 1], [1 0; 0 1])
 
-@test c2d(ss(4*[1 0; 0 1]), 0.5, :zoh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
-@test c2d(ss(4*[1 0; 0 1]), 0.5, :foh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
-@test_c2d(c2d(C_111, 0.01, :zoh),
-    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0], true)
-@test_c2d(c2d(C_111, 0.01, :foh),
+@test c2d_x0map(ss(4*[1 0; 0 1]), 0.5, :zoh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
+@test c2d_x0map(ss(4*[1 0; 0 1]), 0.5, :foh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
+@test_tupleapprox(c2d_x0map(C_111, 0.01, :zoh),
+    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0])
+@test_tupleapprox(c2d_x0map(C_111, 0.01, :foh),
     ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01),
-    [1 -0.009835396005712075], true)
-@test_c2d(c2d(C_212, 0.01, :zoh),
+    [1 -0.009835396005712075])
+@test_tupleapprox(c2d_x0map(C_212, 0.01, :zoh),
     ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
     [0.009466805409666932; 0.019219966830212765], [1 0; 0 1], [0; 0], 0.01),
-    [1 0 0; 0 1 0], true)
-@test_c2d(c2d(C_212, 0.01, :foh),
+    [1 0 0; 0 1 0])
+@test_tupleapprox(c2d_x0map(C_212, 0.01, :foh),
     ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
     [0.008957940478201584; 0.018468989584974498], [1 0; 0 1], [0.004820885889482196;
-    0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675], true)
-@test_c2d(c2d(C_221, 0.01, :zoh),
+    0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675])
+@test_tupleapprox(c2d_x0map(C_221, 0.01, :zoh),
     ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
     [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
-    [1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0], true)
-@test_c2d(c2d(C_221, 0.01, :foh),
+    [1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0])
+@test_tupleapprox(c2d_x0map(C_221, 0.01, :foh),
     ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
     [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
     [1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01),
     [1 0 -0.004917457305816479 9.657141633428213e-5;
-    0 1 -3.219047211142736e-5 -0.009706152723187249], true)
-@test_c2d(c2d(C_222_d, 0.01, :zoh),
+    0 1 -3.219047211142736e-5 -0.009706152723187249])
+@test_tupleapprox(c2d_x0map(C_222_d, 0.01, :zoh),
     ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
     [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
-    [1 0; 0 1], [1 0; 0 1], 0.01), [1 0 0 0; 0 1 0 0], true)
-@test_c2d(c2d(C_222_d, 0.01, :foh),
+    [1 0; 0 1], [1 0; 0 1], 0.01), [1 0 0 0; 0 1 0 0])
+@test_tupleapprox(c2d_x0map(C_222_d, 0.01, :foh),
     ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
     [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
     [1 0; 0 1], [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01),
-    [1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249], true)
-
-# Test some false
-@test_c2d(c2d(2C_111, 0.01, :zoh),  # Factor 2
-    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0], false)
-@test_c2d(c2d(C_111, 0.01, :foh),   #Wrong C
-    ss([0.951229424500714], [0.01902855227625244], [3*3], [0.029506188017136226], 0.01),
-    [1 -0.009835396005712075], false)
-@test_c2d(c2d(C_212, 0.1, :zoh),    # Wrong Ts
-    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
-    [0.009466805409666932; 0.019219966830212765], [1 0; 0 1], [0; 0], 0.01),
-    [1 0 0; 0 1 0], false)
+    [1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249])
 
 # Test c2d for transfer functions
 G = tf([1, 1], [1, 3, 1])
 Gd = c2d(G, 0.2)
 @test Gd ≈ tf([0, 0.165883310712090, -0.135903621603238], [1.0, -1.518831946985175, 0.548811636094027], 0.2) rtol=1e-14
 
 # Issue #391
-@test c2d(tf(C_111), 0.01, :zoh) ≈ tf(c2d(C_111, 0.01, :zoh)[1])
-@test c2d(tf(C_111), 0.01, :foh) ≈ tf(c2d(C_111, 0.01, :foh)[1])
+@test c2d(tf(C_111), 0.01, :zoh) ≈ tf(c2d(C_111, 0.01, :zoh))
+@test c2d(tf(C_111), 0.01, :foh) ≈ tf(c2d(C_111, 0.01, :foh))
 
 # c2d on a zpk model should arguably return a zpk model
 @test_broken typeof(c2d(zpk(G), 1)) <: TransferFunction{<:ControlSystems.SisoZpk}
@@ -72,27 +61,27 @@ Gd = c2d(G, 0.2)
 
 # d2c
 @static if VERSION > v"1.4" # log(matrix) is buggy on previous versions, should be fixed in 1.4 and back-ported to 1.0.6
-    @test d2c(c2d(C_111, 0.01)[1]) ≈ C_111
-    @test d2c(c2d(C_212, 0.01)[1]) ≈ C_212
-    @test d2c(c2d(C_221, 0.01)[1]) ≈ C_221
-    @test d2c(c2d(C_222_d, 0.01)[1]) ≈ C_222_d
+    @test d2c(c2d(C_111, 0.01)) ≈ C_111
+    @test d2c(c2d(C_212, 0.01)) ≈ C_212
+    @test d2c(c2d(C_221, 0.01)) ≈ C_221
+    @test d2c(c2d(C_222_d, 0.01)) ≈ C_222_d
     @test d2c(Gd) ≈ G
 
     sys = ss([0 1; 0 0], [0;1], [1 0], 0)
-    sysd = c2d(sys, 1)[1]
+    sysd = c2d(sys, 1)
     @test d2c(sysd) ≈ sys
 end
 
 # forward euler
-@test c2d(C_111, 1, :fwdeuler)[1].A == I + C_111.A
-@test d2c(c2d(C_111, 0.01, :fwdeuler)[1], :fwdeuler) ≈ C_111
-@test d2c(c2d(C_212, 0.01, :fwdeuler)[1], :fwdeuler) ≈ C_212
-@test d2c(c2d(C_221, 0.01, :fwdeuler)[1], :fwdeuler) ≈ C_221
-@test d2c(c2d(C_222_d, 0.01, :fwdeuler)[1], :fwdeuler) ≈ C_222_d
+@test c2d(C_111, 1, :fwdeuler).A == I + C_111.A
+@test d2c(c2d(C_111, 0.01, :fwdeuler), :fwdeuler) ≈ C_111
+@test d2c(c2d(C_212, 0.01, :fwdeuler), :fwdeuler) ≈ C_212
+@test d2c(c2d(C_221, 0.01, :fwdeuler), :fwdeuler) ≈ C_221
+@test d2c(c2d(C_222_d, 0.01, :fwdeuler), :fwdeuler) ≈ C_222_d
 @test d2c(c2d(G, 0.01, :fwdeuler), :fwdeuler) ≈ G
 
 
-Cd = c2d(C_111, 0.001, :fwdeuler)[1]
+Cd = c2d(C_111, 0.001, :fwdeuler)
 t = 0:0.001:2
 y,_ = step(C_111, t)
 yd,_ = step(Cd, t)
@@ -130,4 +119,3 @@ p = pole(Gcl)
 
 
 end
-

test/test_timeresp.jl

@@ -18,7 +18,7 @@ y, t, x, uout = lsim(sys,u,t,x0=x0) # Continuous time
 th = 1e-6
 @test sum(abs.(x[end,:])) < th
 
-y, t, x, uout = lsim(c2d(sys,0.1)[1],u,t,x0=x0) # Discrete time
+y, t, x, uout = lsim(c2d(sys,0.1),u,t,x0=x0) # Discrete time
 @test sum(abs.(x[end,:])) < th
 
 #Do a manual simulation with uout
@@ -28,7 +28,7 @@ ym, tm, xm = lsim(sys, uout, t, x0=x0)
 
 # Now compare to closed loop
 # Discretization is needed before feedback
-sysd = c2d(sys, 0.1)[1]
+sysd = c2d(sys, 0.1)
 # Create the closed loop system
 sysdfb = ss(sysd.A-sysd.B*L, sysd.B, sysd.C, sysd.D, 0.1)
 #Simulate without input
@@ -90,8 +90,8 @@ yd, td, xd = lsim(sysd, u, t, x0=x0)
 
 # Test lsim with default time vector
 uv = randn(length(t))
-y,t = lsim(c2d(sys,0.1)[1],uv,t,x0=x0)
-yd,td = lsim(c2d(sys,0.1)[1],uv,x0=x0)
+y,t = lsim(c2d(sys,0.1),uv,t,x0=x0)
+yd,td = lsim(c2d(sys,0.1),uv,x0=x0)
 @test yd == y
 @test td == t
 

test/test_transferfunction.jl

@@ -105,9 +105,9 @@ tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005)
 
 # Printing
 res = ("TransferFunction{Continuous,ControlSystems.SisoRational{Int64}}\nInput 1 to output 1\ns^2 + 2s + 3\n-------------\ns^2 + 8s + 15\n\nInput 1 to output 2\ns^2 + 2s + 3\n-------------\ns^2 + 8s + 15\n\nInput 2 to output 1\n    s + 2\n-------------\ns^2 + 8s + 15\n\nInput 2 to output 2\n    s + 2\n-------------\ns^2 + 8s + 15\n\nContinuous-time transfer function model")
-@test sprint(show, C_222) == res
+@test dropwhitespace(sprint(show, C_222)) == dropwhitespace(res)
 res = ("TransferFunction{Discrete{Float64},ControlSystems.SisoRational{Float64}}\nInput 1 to output 1\n1.0z^2 + 2.0z + 3.0\n--------------------\n1.0z^2 - 0.2z - 0.15\n\nInput 1 to output 2\n1.0z^2 + 2.0z + 3.0\n--------------------\n1.0z^2 - 0.2z - 0.15\n\nInput 2 to output 1\n     1.0z + 2.0\n--------------------\n1.0z^2 - 0.2z - 0.15\n\nInput 2 to output 2\n     1.0z + 2.0\n--------------------\n1.0z^2 - 0.2z - 0.15\n\nSample Time: 0.005 (seconds)\nDiscrete-time transfer function model")
-@test sprint(show, D_222) == res
+@test dropwhitespace(sprint(show, D_222)) == dropwhitespace(res)
 
 
 @test tf(zpk([1.0 2; 3 4])) == tf([1 2; 3 4])