Add support for pure D systems in lsim function and corresponding tests

lib/ControlSystemsBase/src/timeresp.jl

@@ -196,6 +196,14 @@ function lsim(sys::AbstractStateSpace, u::AbstractVecOrMat, t::AbstractVector;
         error("Time vector t must be uniformly spaced")
     end
 
+    # Handle pure D system (no states)
+    if nx == 0
+        x = Matrix{eltype(u)}(undef, 0, length(t)) # no states
+        y = sys.D * u
+        dsys = sys
+        return SimResult(y, t, x, u, dsys)
+    end
+
     if iscontinuous(sys)
         if method === :zoh
             dsys = c2d(sys, dt, :zoh)
@@ -280,6 +288,18 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
 
     dt = t[2] - t[1]
 
+    # Handle pure D system (no states)
+    if nx == 0
+        uout = Matrix{T}(undef, nu, length(t))
+        for i = 1:length(t)
+            uout[:, i] = u(T[], t[i])
+        end
+        x = Matrix{T}(undef, 0, length(t)) # no states
+        y = sys.D * uout
+        simsys = sys
+        return SimResult(y, t, x, uout, simsys)
+    end
+
     if !iscontinuous(sys) || method ∈ (:zoh, :tustin, :foh, :fwdeuler)
         if iscontinuous(sys)
             simsys = c2d(sys, dt, method)

src/timeresp.jl

@@ -37,6 +37,18 @@ function lsim(sys::AbstractStateSpace{Continuous}, u::Function, t::AbstractVecto
 
     dt = t[2] - t[1]
 
+    # Handle pure D system (no state)
+    if nx == 0
+        uout = Matrix{T}(undef, nu, length(t))
+        for i = 1:length(t)
+            uout[:, i] = u(T[], t[i])
+        end
+        x = Matrix{T}(undef, 0, length(t)) # no states
+        y = sys.D * uout
+        simsys = sys
+        return SimResult(y, t, x, uout, simsys)
+    end
+
     if method === :zoh
         simsys = c2d(sys, dt, :zoh)
         x,uout = ControlSystemsBase.ltitr(simsys.A, simsys.B, u, t, T.(x0))

test/test_timeresp.jl

@@ -1,3 +1,9 @@
+# Test lsim for pure D system (no states) with matrix input
+sysDm = ss([2.0])
+tDm = 0:0.1:1
+uDm = fill(3.0, 1, length(tDm))
+resDm = lsim(sysDm, uDm, tDm)
+@test resDm.y[:] == fill(2.0 * 3.0, length(tDm))
 import OrdinaryDiffEq
 using OrdinaryDiffEq: solve
 using ControlSystems
@@ -35,6 +41,14 @@ th = 1e-6
 # Error for nonuniformly spaced vector
 @test_throws ErrorException lsim(sys, [1 2 3 4], [0, 1, 1, 2])
 
+
+# Test lsim for pure D system (no states)
+sysD = ss([2.0])
+tD = 0:0.1:1
+uD(x, t) = [3.0]
+resD = lsim(sysD, uD, tD)
+@test resD.y[:] == fill(2.0 * 3.0, length(tD))
+
 # Test for problem with broadcast dimensions
 @test lsim(sys, zeros(1, 5), 0:0.2:0.8)[1][:] == zeros(5)