Return a result structure from lsim

This is by design a non-breaking change to the lsim interface since the structure allows both getindex and destructuring to behave just like if lsim returned a tuple of arrays like before. Indeed, the tests for lsim were not touched in this commit. This commit also adds a plot recipe to the result structure. All plot recipes for lsimplot, stepplot, impulseplot have been replaced by the new recipe. This is a breaking change since the names of the previous plots no longer exist. A slight change is that the plots for a step response no longer show the text "step response", but I think that's an acceptable change, the user can supply any title they prefer themselves.
baggepinnen · Oct 20, 2021 · May 29, 2021 · May 30, 2021 · May 30, 2021 · May 30, 2021
commit a9cf9a75d2bf9b70053de4bc539378f67c0b9448
diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl
@@ -116,6 +116,8 @@ using MacroTools
 
 abstract type AbstractSystem end
 
+
+include("types/result_types.jl")
 include("types/TimeEvolution.jl")
 ## Added interface:
 #   timeevol(Lti) -> TimeEvolution (not exported)

diff --git a/src/plotting.jl b/src/plotting.jl
@@ -113,110 +113,39 @@ function getLogTicks(x, minmax)
 end
 
 
-@userplot Lsimplot
-
-"""
-    fig = lsimplot(sys::LTISystem, u, t; x0=0)
-    lsimplot(LTISystem[sys1, sys2...], u, t; x0)
-
-Calculate the time response of the `LTISystem`(s) to input `u`. If `x0` is
-not specified, a zero vector is used.
-"""
-lsimplot
-
-@recipe function lsimplot(p::Lsimplot; x0=0)
-
-    systems = p.args[1]
-
-    if !isa(systems,AbstractArray)
-        systems = [systems]
-    end
-    if !_same_io_dims(systems...)
-        error("All systems must have the same input/output dimensions")
-    end
-    ny, nu = size(systems[1])
-    layout --> (ny,1)
-    s2i(i,j) = LinearIndices((ny,1))[j,i]
-    for (si,s) in enumerate(systems)
-        s = systems[si]
-        y, t = lsim(s, p.args[2:end]...; x0 = x0 == 0 ? zeros(nstates(s)) : x0)
-        seriestype := iscontinuous(s) ? :path : :steppost
+# This will be called on plot(lsim(sys, args...))
+@recipe function simresultplot(r::SimResult; plotu=false)
+    ny, nu = r.ny, r.nu
+    t = r.t
+    multiseries = size(r.y, 3) # step and impulse produce multiple results
+    layout --> ((plotu ? ny + nu : ny), 1)
+    # s2i(i,j) = LinearIndices((ny,1))[j,i]
+    seriestype := iscontinuous(r.sys) ? :path : :steppost
+    for ms in 1:multiseries
         for i=1:ny
-            ytext = (ny > 1) ? "Amplitude to: y($i)" : "Amplitude"
+            ytext = (ny > 1) ? "y($i)" : "y"
             @series begin
                 xguide  --> "Time (s)"
                 yguide  --> ytext
                 title   --> "System Response"
-                subplot --> s2i(1,i)
-                label     --> "\$G_{$(si)}\$"
-                t,  y[i, :]
+                label   --> (multiseries > 1 ? "From u($(ms))" : "")
+                subplot --> i#s2i(1,i)
+                t,  r.y[i, :, ms]
             end
         end
-    end
-end
-
-@userplot Stepplot
-@userplot Impulseplot
-"""
-    stepplot(sys[, tfinal]; kwargs...) or stepplot(sys[, t]; kwargs...)
-
-Plot step response of  `sys` until final time `tfinal` or at time points in the vector `t`.
-If not defined, suitable values are chosen based on `sys`.
-See also [`step`](@ref)
-
-`kwargs` is sent as argument to Plots.plot.
-"""
-stepplot
-
-"""
-    impulseplot(sys[, tfinal]; kwargs...) or impulseplot(sys[, t]; kwargs...)
-
-Plot impulse response of `sys` until final time `tfinal` or at time points in the vector `t`.
-If not defined, suitable values are chosen based on `sys`.
-See also [`impulse`](@ref)
-
-`kwargs` is sent as argument to Plots.plot.
-"""
-impulseplot
-
-for (func, title, typ) = ((step, "Step Response", Stepplot), (impulse, "Impulse Response", Impulseplot))
-    funcname = Symbol(func,"plot")
-
-    @recipe function f(p::typ)
-        systems = p.args[1]
-        if !isa(systems, AbstractArray)
-            systems = [systems]
-        end
-        if !_same_io_dims(systems...)
-            error("All systems must have the same input/output dimensions")
-        end
-        ny, nu = size(systems[1])
-        layout --> (ny,nu)
-        titles = fill("", 1, ny*nu)
-        title --> titles
-        s2i(i,j) = LinearIndices((ny,nu))[i,j]
-        for (si,s) in enumerate(systems)
-            y,t = func(s, p.args[2:end]...)
-            for i=1:ny
-                for j=1:nu
-                    ydata = reshape(y[i, :, j], size(t, 1))
-                    style = iscontinuous(s) ? :path : :steppost
-                    ttext = (nu > 1 && i==1) ? title*" from: u($j) " : title
-                    titles[s2i(i,j)] = ttext
-                    ytext = (ny > 1 && j==1) ? "Amplitude to: y($i)" : "Amplitude"
-                    @series begin
-                        seriestype --> style
-                        xguide --> "Time (s)"
-                        yguide --> ytext
-                        subplot --> s2i(i,j)
-                        label --> "\$G_{$(si)}\$"
-                        t, ydata
-                    end
-                end
+    end 
+    if plotu # bug in recipe syste, can't use `plotu || return`
+        for i=1:nu
+            utext = (nu > 1) ? "u($i)" : "u"
+            @series begin
+                xguide  --> "Time (s)"
+                yguide  --> utext
+                subplot --> ny+i
+                label --> ""
+                t,  r.u[i, :]
             end
         end
     end
-
 end
 
 """

diff --git a/src/timeresp.jl b/src/timeresp.jl
@@ -15,18 +15,19 @@ function Base.step(sys::AbstractStateSpace, t::AbstractVector; method=:cont, kwa
     T = promote_type(eltype(sys.A), Float64)
     ny, nu = size(sys)
     nx = nstates(sys)
-    u = (x,t)->[one(eltype(t))]
+    u_element = [one(eltype(t))] # to avoid allocating this multiple times
+    u = (x,t)->u_element
     x0 = zeros(nx)
     if nu == 1
-        y, tout, x, _ = lsim(sys, u, t; x0=x0, method=method, kwargs...)
+        y, tout, x, uout = lsim(sys, u, t; x0, method, kwargs...)
     else
         x = Array{T}(undef, nx, length(t), nu)
         y = Array{T}(undef, ny, length(t), nu)
         for i=1:nu
-            y[:,:,i], tout, x[:,:,i],_ = lsim(sys[:,i], u, t; x0=x0, method=method, kwargs...)
+            y[:,:,i], tout, x[:,:,i], uout = lsim(sys[:,i], u, t; x0, method, kwargs...)
         end
     end
-    return y, t, x
+    return SimResult(y, t, x, uout, sys)
 end
 
 Base.step(sys::LTISystem, tfinal::Real; kwargs...) = step(sys, _default_time_vector(sys, tfinal); kwargs...)
@@ -42,7 +43,7 @@ vector `t` is not provided, one is calculated based on the system pole
 locations.
 
 `y` has size `(ny, length(t), nu)`, `x` has size `(nx, length(t), nu)`"""
-function impulse(sys::AbstractStateSpace, t::AbstractVector; method=:cont, kwargs...)
+function impulse(sys::AbstractStateSpace, t::AbstractVector; kwargs...)
     T = promote_type(eltype(sys.A), Float64)
     ny, nu = size(sys)
     nx = nstates(sys)
@@ -58,29 +59,33 @@ function impulse(sys::AbstractStateSpace, t::AbstractVector; method=:cont, kwarg
         x0s = zeros(T, nx, nu)
     end
     if nu == 1 # Why two cases # QUESTION: Not type stable?
-        y, t, x,_ = lsim(sys, u, t; x0=x0s[:], method=method, kwargs...)
+        y, t, x,uout = lsim(sys, u, t; x0=x0s[:], kwargs...)
-        y, t, x,uout = lsim(sys, u, t; x0=x0s[:], kwargs...)
+        y, t, x, uout = lsim(sys, u, t; x0=x0s[:], kwargs...)
-        y, t, x,uout = lsim(sys, u, t; x0=x0s[:], kwargs...)
+        y, t, x, uout = lsim(sys, u, t; x0=x0s[:], kwargs...)
     else
         x = Array{T}(undef, nx, length(t), nu)
         y = Array{T}(undef, ny, length(t), nu)
         for i=1:nu
-            y[:,:,i], t, x[:,:,i],_ = lsim(sys[:,i], u, t; x0=x0s[:,i], method=method, kwargs...)
+            y[:,:,i], t, x[:,:,i],uout = lsim(sys[:,i], u, t; x0=x0s[:,i], kwargs...)
-            y[:,:,i], t, x[:,:,i],uout = lsim(sys[:,i], u, t; x0=x0s[:,i], kwargs...)
+            y[:,:,i], t, x[:,:,i], uout = lsim(sys[:,i], u, t; x0=x0s[:,i], kwargs...)
-            y[:,:,i], t, x[:,:,i],uout = lsim(sys[:,i], u, t; x0=x0s[:,i], kwargs...)
+            y[:,:,i], t, x[:,:,i], uout = lsim(sys[:,i], u, t; x0=x0s[:,i], kwargs...)
         end
     end
-    return y, t, x
+    return SimResult(y, t, x, uout, sys)
 end
 
 impulse(sys::LTISystem, tfinal::Real; kwargs...) = impulse(sys, _default_time_vector(sys, tfinal); kwargs...)
 impulse(sys::LTISystem; kwargs...) = impulse(sys, _default_time_vector(sys); kwargs...)
 impulse(sys::TransferFunction, t::AbstractVector; kwargs...) = impulse(ss(sys), t; kwargs...)
 
 """
-    y, t, x = lsim(sys, u[, t]; x0, method])
-    y, t, x, uout = lsim(sys, u::Function, t; x0, method)
+    result = lsim(sys, u[, t]; x0, method])
+    result = lsim(sys, u::Function, t; x0, method)
 
 Calculate the time response of system `sys` to input `u`. If `x0` is ommitted,
 a zero vector is used.
 
-`y`, `x`, `uout` has time in the second dimension. Initial state `x0` defaults to zero.
+The result structure contains the fields `y, t, x, u` and can be destructured automatically by iteration, e.g.,
+```julia
+y, t, x, u = result
+```
+`y, `x`, `u` have time in the second dimension. Initial state `x0` defaults to zero.
 
 Continuous time systems are simulated using an ODE solver if `u` is a function. If `u` is an array, the system is discretized before simulation. For a lower level inteface, see `?Simulator` and `?solve`
 
@@ -136,7 +141,7 @@ function lsim(sys::AbstractStateSpace, u::AbstractVecOrMat, t::AbstractVector;
             dsys, x0map = c2d_x0map(sys, dt, :foh)
             x0 = x0map*[x0; u[:,1]]
         else
-            error("Unsupported discretization method")
+            error("Unsupported discretization method: $method")
         end
     else
         if sys.Ts != dt
@@ -147,16 +152,16 @@ function lsim(sys::AbstractStateSpace, u::AbstractVecOrMat, t::AbstractVector;
 
     x = ltitr(dsys.A, dsys.B, u, x0)
     y = sys.C*x + sys.D*u
-    return y, t, x
+    return SimResult(y, t, x, u, dsys) # saves the system that actually produced the simulation
 end
 
 function lsim(sys::AbstractStateSpace{<:Discrete}, u::AbstractVecOrMat; kwargs...)
     t = range(0, length=size(u, 2), step=sys.Ts)
     lsim(sys, u, t; kwargs...)
 end
 
-@deprecate lsim(sys, u, t, x0) lsim(sys, u, t; x0=x0)
-@deprecate lsim(sys, u, t, x0, method) lsim(sys, u, t; x0=x0, method=method)
+@deprecate lsim(sys, u, t, x0) lsim(sys, u, t; x0)
+@deprecate lsim(sys, u, t, x0, method) lsim(sys, u, t; x0, method)
 
 function lsim(sys::AbstractStateSpace, u::Function, tfinal::Real; kwargs...)
     t = _default_time_vector(sys, tfinal)
@@ -185,22 +190,23 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
 
     if !iscontinuous(sys) || method === :zoh
         if iscontinuous(sys)
-            dsys = c2d(sys, dt, :zoh)
+            simsys = c2d(sys, dt, :zoh)
         else
             if sys.Ts != dt
                 error("Time vector must match sample time for discrete system")
             end
-            dsys = sys
+            simsys = sys
         end
-        x,uout = ltitr(dsys.A, dsys.B, u, t, T.(x0))
+        x,uout = ltitr(simsys.A, simsys.B, u, t, T.(x0))
     else
         p = (sys.A, sys.B, u)
         sol = solve(ODEProblem(f_lsim, x0, (t[1], t[end]), p), alg; saveat=t, kwargs...)
         x = reduce(hcat, sol.u)
         uout = reduce(hcat, u(x[:, i], t[i]) for i in eachindex(t))
+        simsys = sys
     end
     y = sys.C*x + sys.D*uout
-    return y, t, x, uout
+    return SimResult(y, t, x, uout, simsys) # saves the system that actually produced the simulation
 end
 
 

diff --git a/src/types/result_types.jl b/src/types/result_types.jl
@@ -0,0 +1,43 @@
+abstract type AbstractResult end
+
+## SimResult: the output of lsim etc. ==========================================
+
+struct SimResult{Ty, Tt, Tx, Tu, Ts} <: AbstractResult
+    y::Ty
+    t::Tt
+    x::Tx
+    u::Tu
+    sys::Ts
+end
+
+# To emulate, e.g., lsim(sys, u)[1] -> y
+function Base.getindex(r::SimResult, i::Int)
+    return getfield(r, i) 
+end
+
+function Base.getindex(r::SimResult, v::AbstractVector)
+    return getfield.((r,), v) 
+end
+
+# to allow destructuring, e.g., y,t,x = lsim(sys, u)
+# This performs explicit iteration in the type domain to ensure inferability
+Base.iterate(r::SimResult)              = (r.y, Val(:t))
+Base.iterate(r::SimResult, ::Val{:t})   = (r.t, Val(:x))
+Base.iterate(r::SimResult, ::Val{:x})   = (r.x, Val(:u))
+Base.iterate(r::SimResult, ::Val{:u})   = (r.u, Val(:sys))
+Base.iterate(r::SimResult, ::Val{:sys}) = (r.sys, Val(:done))
+Base.iterate(r::SimResult, ::Val{:done}) = nothing
+
+
+function Base.getproperty(r::SimResult, s::Symbol)
+    s ∈ fieldnames(SimResult) && (return getfield(r, s))
+    if s === :nx
+        return size(r.x, 1)
+    elseif s === :nu
+        return size(r.u, 1)
+    elseif s === :ny
+        return size(r.y, 1)
+    else
+        throw(ArgumentError("Unsupported property $s"))
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -6,6 +6,7 @@ import DSP: conv            # In test_conversion and test_synthesis
 include("framework.jl")
 
 my_tests = [
+            "test_result_types",
             "test_timeevol",
             "test_statespace",
             "test_transferfunction",            

diff --git a/test/test_plots.jl b/test/test_plots.jl
@@ -28,10 +28,10 @@ function getexamples()
     #Only siso for now
     nicholsgen() = nicholsplot(tf1,ws)
 
-    stepgen() = stepplot(sys, ts[end], l=(:dash, 4))
-    impulsegen() = impulseplot(sys, ts[end], l=:blue)
+    stepgen() = plot(step(sys, ts[end]), l=(:dash, 4))
+    impulsegen() = plot(impulse(sys, ts[end]), l=:blue)
     L = lqr(sysss.A, sysss.B, [1 0; 0 1], [1 0; 0 1])
-    lsimgen() = lsimplot(sysss, (x,i)->-L*x, ts; x0=[1;2])
+    lsimgen() = plot(lsim(sysss, (x,i)->-L*x, ts; x0=[1;2]), plotu=true)
 
     margingen() = marginplot([tf1, tf2], ws)
     gangoffourgen() = begin

diff --git a/test/test_result_types.jl b/test/test_result_types.jl
@@ -0,0 +1,40 @@
+@testset "SimResult" begin
+    @info "Testing SimResult"
+
+    import ControlSystems: SimResult, ssrand
+
+    y = randn(1,10)
+    u = randn(1,10)
+    x = randn(1,10)
+    t = 1:10
+    sys = ssrand(1, 1, 1, Ts=1)
+
+    r = SimResult(y,t,x,u,sys)
+
+    # test getindex
+    @test r[1] == y
+    @test r[2] == t
+    @test r[3] == x
+    @test r[4] == u
+    @test r[5] == sys
+
+    # test destructuring
+    y2,t2,x2,u2 = r # incomplete destructuring
+    @test y2 === y
+    @test t2 === t
+    @test x2 === x
+    @test u2 === u
+
+    _,_,_,_,sys2 = r # complete destructuring
+    @test sys2 === sys
+
+    # test inferability of destructuring
+    foo(r::SimResult) = (a,b,c,d,e) = r
+    @inferred foo(r)
+
+    # test properties
+    @test r.nx == 1
+    @test r.nu == 1
+    @test r.ny == 1
+end
+