Merge pull request #945 from JuliaControl/hcatres

add concatenation for SimResult

lib/ControlSystemsBase/src/timeresp.jl

@@ -144,7 +144,7 @@ Continuous-time systems are simulated using an ODE solver if `u` is a function (
 If `u` is a function, then `u(x,i)` (for discrete systems) or `u(x,t)` (for continuous ones) is called to calculate the control signal at every iteration (time instance used by solver). This can be used to provide a control law such as state feedback `u(x,t) = -L*x` calculated by `lqr`.
 To simulate a unit step at `t=t₀`, use `(x,t)-> t ≥ t₀`, for a ramp, use `(x,t)-> t`, for a step at `t=5`, use `(x,t)-> (t >= 5)` etc.
 
-*Note:* The function `u` will be called once before simulating to verify that it returns an array of the correct dimensions. This can cause problems if `u` is stateful. You can disable this check by passing `check_u = false`.
+*Note:* The function `u` will be called once before simulating to verify that it returns an array of the correct dimensions. This can cause problems if `u` is stateful or has other side effects. You can disable this check by passing `check_u = false`.
 
 For maximum performance, see function [`lsim!`](@ref), available for discrete-time systems only.
 
@@ -206,7 +206,7 @@ function lsim(sys::AbstractStateSpace, u::AbstractVecOrMat, t::AbstractVector;
             error("Unsupported discretization method: $method")
         end
     else
-        if sys.Ts != dt
+        if !(sys.Ts ≈ dt)
             error("Time vector must match the sample time of the discrete-time system, $(sys.Ts): got $dt")
         end
         dsys = sys
@@ -284,8 +284,8 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
         if iscontinuous(sys)
             simsys = c2d(sys, dt, :zoh)
         else
-            if sys.Ts != dt
-                error("Time vector must match sample time for discrete system")
+            if !(sys.Ts ≈ dt)
+                error("Time vector interval ($dt) must match sample time for discrete system ($(sys.Ts))")
             end
             simsys = sys
         end

lib/ControlSystemsBase/src/types/result_types.jl

@@ -24,6 +24,14 @@ y, t, x, u = result
 - `x::Tx`
 - `u::Tu`
 - `sys::Ts`
+
+## Concatenation of SimResults
+
+Two SimResults can be concatenated in time using `hcat`, or `[res1 res2]`, the rules for this are as follows:
+- If the start time of the second result is one sample interval after the end time of the first result, the results are concatenated and the length of the result is the sum of the lengths of the two results.
+- If the start time of the second result is equal to the end time of the first result, _and_ the initial state of the second result is equal to the final state of the first result, the results are concatenated omitting the initial point from the second result, which would otherwise have been repeated. The length of the result is the sum of the lengths of the two results minus one.
+If none of the above holds, a warning is issued and the result has the length of the sum of the lengths of the two results.
+- If the sample intervals of the two results are different, an error is thrown.
 """
 struct SimResult{Ty, Tt, Tx, Tu, Ts} <: AbstractSimResult # Result of lsim
     y::Ty
@@ -42,6 +50,18 @@ function Base.getindex(r::SimResult, v::AbstractVector)
     return getfield.((r,), v) 
 end
 
+function Base.hcat(r1::SimResult, r2::SimResult)
+    r1.sys.Ts == r2.sys.Ts || throw(ArgumentError("Sampling-time mismatch"))
+    if r1.x[:, end] == r2.x[:, 1] && r1.t[end] == r2.t[1]
+        r1.u[:, end] == r2.u[:, 1] || @warn "Concatenated SimResults have different inputs at the join"
+        # This is a common case when r2 starts with initial conditions from r1
+        return SimResult(hcat(r1.y, r2.y[:, 2:end]), vcat(r1.t, r2.t[2:end]), hcat(r1.x, r2.x[:, 2:end]), hcat(r1.u, r2.u[:, 2:end]), r1.sys)
+    elseif !(r1.t[end] + r1.sys.Ts ≈ r2.t[1])
+        @warn "Concatenated SimResults do not appear to be continuous in time, the first ends at t=$(r1.t[end]) and the second starts at t=$(r2.t[1]). With sample interval Ts=$(r1.sys.Ts), the second simulation was expected to start at t=$(r1.t[end] + r1.sys.Ts) To start a simulation at a non-zero time, pass a time vector to lsim."
+    end
+    SimResult(hcat(r1.y, r2.y), vcat(r1.t, r2.t), hcat(r1.x, r2.x), hcat(r1.u, r2.u), r1.sys)
+end
+
 issiso(r::SimResult) = issiso(r.sys)
 
 # to allow destructuring, e.g., y,t,x = lsim(sys, u)

lib/ControlSystemsBase/test/test_timeresp.jl

@@ -255,4 +255,24 @@ si = stepinfo(res)
 @test si.undershoot ≈ 27.98 rtol=0.01
 plot(si)
 
+# Test concatenation of SimResults
+u = ones(1, 100)
+sysd = c2d(sys,0.1)
+res1 = lsim(sysd,u)
+res2 = lsim(sysd,u; x0 = res1.x[:, end])
+@test_logs (:warn, r"Concatenated SimResults do not appear to be continuous in time") [res1 res2]
+
+res2 = lsim(sysd,u,res1.t[end]:0.1:res1.t[end]+9.9; x0 = res1.x[:, end])
+res12 = [res1 res2]
+@test length(res12.t) == length(res1.t) + length(res2.t) - 1 # -1 since we do not include the initial time point from the second result which overlaps with the first
+
+res2 = lsim(sysd,u)
+@test_logs (:warn, r"Concatenated SimResults do not appear to be continuous in time") [res1 res2]
+res12 = [res1 res2]
+@test length(res12.t) == length(res1.t) + length(res2.t) # not -1 since we do include the initial time point from the second result if they do not appear to be continuous in time
+
+res2 = lsim(sysd,u, res1.t[end]+0.1:0.1:res1.t[end]+10)
+@test_nowarn res12 = [res1 res2]
+@test length(res12.t) == length(res1.t) + length(res2.t) # not -1 since we do do include the initial time point from the second result if they do not appear to be continuous in time
+
 end